Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Start-to-Finish Example: Applying Boundary Conditions in Curvilinear Coordinates in Three Dimensions\n",
    "\n",
    "## Authors: Zach Etienne & Terrence Pierre Jacques\n",
    "\n",
    "## This paper validates boundary condition algorithms for curvilinear coordinates like Spherical and Cylindrical, as mentioned in the [SENR/NRPy+ paper](https://arxiv.org/abs/1712.07658). Using an Eigen-Coordinate system, it differentiates between inner and outer boundaries, applying parity conditions to arbitrary-rank tensors. The approach is validated by successfully navigating coordinate singularities through a combination of cell-centered grids, tensor rescaling, and a stable Method of Lines time stepping algorithm.\n",
    "\n",
    "<a id='intro'></a>\n",
    "\n",
    "## Introduction: Applying boundary conditions in curvilinear coordinates\n",
    "$$\\label{intro}$$\n",
    "\n",
    "We face four challenges when solving PDEs in curvilinear coordinates:\n",
    "\n",
    "1. [Challenge 1](#innerouterbcs): Unlike ordinary Cartesian coordinate grids, not all boundary points on uniform curvilinear coordinate grids are outer boundary points.\n",
    "1. [Challenge 2](#coordinversion): Figuring out the locations to which boundary points map requires a coordinate inversion, but some coordinate systems are not easily invertible.\n",
    "1. [Challenge 3](#parity): Tensors and vectors in curvilinear coordinates can change directions across \"inner\" boundaries, changing sign as a result. \n",
    "1. [Challenge 4](#coordsingularities): Coordinate singularities can appear, causing a stiffening or divergence of terms in PDEs."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='innerouterbcs'></a>\n",
    "\n",
    "### Challenge 1: Inner versus outer boundary points\n",
    "$$\\label{innerouterbcs}$$\n",
    "\n",
    "Unlike Cartesian coordinates, the boundaries of our grid in generic curvilinear coordinates are not all outer boundaries. \n",
    "\n",
    "Consider first a computational grid in Cartesian coordinates, with $(x_0,x_1,x_2)=(x,y,z)$, that is *uniform* (i.e., with $\\Delta x$, $\\Delta y$, and $\\Delta z$ all set to constant values). This grid may extend over arbitrary coordinate ranges $x_i \\in [x_{i, \\rm min},x_{i, \\rm max}]$.\n",
    "\n",
    "By contrast, consider now a uniform grid in spherical coordinates $(x_0,x_1,x_2)=(r,\\theta,\\phi)$ with constant spacing $\\Delta r$, $\\Delta \\theta$, and $\\Delta \\phi$ between grid points in $r$, $\\theta$, and $\\phi$, respectively. Further, let's assume that these grids span all possible values of $\\theta$ and $\\phi$, with $r=0$ included in the domain. Then our numerical domain must satisfy the following relations\n",
    "\n",
    "+ $x_0 = r \\in [0,{\\rm RMAX}]$,\n",
    "+ $x_1 = \\theta \\in [0,\\pi]$, and\n",
    "+ $x_2 = \\phi \\in [-\\pi,\\pi]$. (Notice how we do not choose $x_2= \\phi \\in [0,2\\pi]$ so that our conversion from Cartesian to spherical coordinates is compatible with the output range from the ${\\rm atan2}(y,x)$ function: $\\phi={\\rm atan2}(y,x)\\in[-\\pi,\\pi]$.)\n",
    "\n",
    "Notice that, unlike Cartesian coordinates, the boundaries of this numerical grid in spherical coordinates are not all outer boundaries. For example, data stored at $\\phi=-\\pi$ will be identical to data at $\\phi=\\pi$, *regardless of $r$ or $\\theta$*. I.e., $\\phi$ satisfies *periodic* boundary conditions only. Further, $\\theta=0$ presents a more complicated boundary condition, in which points with negative $\\theta$ map to points with $|\\theta|$ but at an angle of $\\phi\\to \\phi+\\pi$. Finally, negative $r$ points will map to positive $r$ points on the other side of the origin. We call these boundaries *inner* boundaries, as they generally map to other points in the interior (as opposed to the outer boundaries) of the grid. \n",
    "\n",
    "As we generally cannot apply an outer boundary condition to the inner boundaries, these boundaries will need to be treated differently.\n",
    "\n",
    "On our numerical grids, this poses some difficulty, as finite difference derivatives we compute within the numerical domain require that the grid be extended beyond the domain boundaries. In spherical coordinates, this means that we need additional grid points at, e.g.,  $r<0$, $\\theta<0$, and $\\phi>\\pi$, just to name a few. Whether they be on outer or inner boundaries, we call grid points in the extended region *ghost zones*.\n",
    "\n",
    "Numerical grids of $N$th order accuracy generally possess $N/2$ ghost zone points in the boundary regions (i.e., $x_i < x_{i,\\rm min}$ and $x_i > x_{i, \\rm max}$). While in Cartesian coordinates, these ghost zone points map to regions outside the grid domain $x_i \\in [x_{i, \\rm min},x_{i, \\rm max}]$, in spherical coordinates, *most* ghost zone points map to regions *inside* the grid domain. For example, for some $\\tilde{r}\\in [0,{\\rm RMAX}]$ and $\\tilde{\\theta}\\in[0,\\pi]$, the ghost zone point $(\\tilde{r},\\tilde{\\theta},2\\pi+\\Delta \\phi/2)$ would map to the interior point $(\\tilde{r},\\tilde{\\theta},\\Delta \\phi/2)$ because the $\\phi$ coordinate is periodic. Thus when given a ghost zone point in some arbitrary curvilinear coordinate system, we are faced with the problem of addressing the following two questions:\n",
    "1. Does a given ghost point map to an interior point, or is it an outer boundary point (i.e., a point exterior to the domain)?\n",
    "1. If the ghost zone point maps to an interior point, to which interior point does it map?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='coordinversion'></a>\n",
    "\n",
    "### Challenge 2: Inverting coordinates\n",
    "$$\\label{coordinversion}$$\n",
    "\n",
    "Coordinate systems within NRPy+ are generally Spherical-like, Cylindrical-like, SymTP-like (where SymTP is a prolate-spheroidal coordinate system), or Cartesian-like. For example, SinhSphericalv2 coordinates are exactly the same as Spherical coordinates, except we choose an odd function for the radial coordinate $r$ as a function of $x_0$:\n",
    "\n",
    "$$\n",
    "r(x_0) = {\\rm AMPL} \\left[ {\\rm const\\_dr} x_0 + \\sinh\\left(\\frac{x_0}{\\rm SINHW}\\right) / \\sinh\\left(\\frac{1}{\\rm SINHW}\\right) \\right].\n",
    "$$\n",
    "\n",
    "While this coordinate choice exhibits nice properties for certain cases, the function $x_0(r)$ is not a closed-form expression. Thus finding the mapping of ghost zone points in the radial direction would require a root finder. \n",
    "\n",
    "*Is there an easier way of dealing with this problem than with a root finder?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='parity'></a>\n",
    "\n",
    "### Challenge 3: Parity: changes of direction in vectors and tensors across inner boundaries\n",
    "$$\\label{parity}$$\n",
    "\n",
    "When applying inner boundary conditions to vectors and tensors, we must consider how the direction or *parity* of vector and tensor components change across the inner boundary.\n",
    "\n",
    "Suppose we have a vector $v^\\rho$ defined at ghost zone $(-\\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\n",
    "\\begin{align}\n",
    "v^\\rho(-\\rho,\\phi,z)&=-v^\\rho(\\rho,\\phi+\\pi,z) \\\\\n",
    "&= \\mathbf{e}^\\rho\\left(-\\rho,\\phi,z\\right) \\cdot \\mathbf{e}^\\rho\\left(\\rho,\\phi+\\pi,z\\right)v^\\rho(\\rho,\\phi+\\pi,z),\n",
    "\\end{align}\n",
    "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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='coordsingularities'></a>\n",
    "\n",
    "### Challenge 4: Coordinate singularities\n",
    "$$\\label{coordsingularities}$$\n",
    "\n",
    "Most non-Cartesian, orthogonal coordinate systems (like spherical coordinates) possess *coordinate singularities*. \n",
    "\n",
    "For example, coordinate singularities in spherical coordinates lie along $\\theta=0$ and $\\theta=\\pi$; these are points where the coordinate system focuses to a single point. For example, the coordinate singularity at the North Pole is the reason why all directions are south there. Critically, these singularities manifest as points where the reference metric or its inverse crosses through zero or diverges to $\\infty$. As we derived in a [previous module](Tutorial-ScalarWaveCurvilinear.ipynb), the Laplacian in spherical polar coordinates takes the form\n",
    "$$\n",
    "\\nabla^2 u = \\partial_r^2 u + \\frac{1}{r^2} \\partial_\\theta^2 u + \\frac{1}{r^2 \\sin^2 \\theta} \\partial_\\phi^2 u +  \\frac{2}{r} \\partial_r u + \\frac{\\cos\\theta}{r^2 \\sin\\theta} \\partial_\\theta u,\n",
    "$$\n",
    "which diverges at $r=0$ and $\\sin\\theta=0-$precisesly at the $\\theta=0$ and $\\theta=\\pi$ coordinate singularity. \n",
    "\n",
    "To avoid this divergence, we simply choose that our numerical grids be **cell-centered**. \n",
    "\n",
    "I.e., given the desired bounds of the grid interior to be \n",
    "\n",
    "\\begin{align}\n",
    "x_0 &\\in [x_{0,\\ \\rm min},x_{0,\\ \\rm max}]\\\\\n",
    "x_1 &\\in [x_{1,\\ \\rm min},x_{1,\\ \\rm max}]\\\\\n",
    "x_2 &\\in [x_{2,\\ \\rm min},x_{2,\\ \\rm max}],\n",
    "\\end{align}\n",
    "\n",
    "${\\rm NGHOSTS}$ to be the number of ghost zones (assumed the same in all directions), and $\\{N_0,N_1,N_2\\}$ to be the desired number of points in the grid interior in the $\\{x_0,x_1,x_2\\}$ directions, respectively, then the numerical grid spacing in each respective direction will be given by\n",
    "\n",
    "\\begin{align}\n",
    "dx_0 &= \\frac{x_{0,\\ \\rm max} - x_{0,\\ \\rm min}}{N_0} \\\\\n",
    "dx_1 &= \\frac{x_{1,\\ \\rm max} - x_{1,\\ \\rm min}}{N_1} \\\\\n",
    "dx_2 &= \\frac{x_{2,\\ \\rm max} - x_{2,\\ \\rm min}}{N_2}.\n",
    "\\end{align}\n",
    "\n",
    "Given the above definitions, the complete set of indices $\\{{\\rm i0},{\\rm i1},{\\rm i2}\\}$ located at $\\{x_{0,{\\rm i0}},x_{1,{\\rm i1}},x_{2,{\\rm i2}}\\}$ as follows:\n",
    "\n",
    "\\begin{align}\n",
    "x_{0,{\\rm i0}} &= x_{0,\\ \\rm min} + \\left[({\\rm i0}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_0 \\\\\n",
    "x_{1,{\\rm i1}} &= x_{1,\\ \\rm min} + \\left[({\\rm i1}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_1 \\\\\n",
    "x_{2,{\\rm i2}} &= x_{2,\\ \\rm min} + \\left[({\\rm i2}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_2 \\\\\n",
    "\\end{align}\n",
    "where\n",
    "\n",
    "* ${\\rm i0}\\in [0,N_0+2\\cdot{\\rm NGHOSTS})$\n",
    "* ${\\rm i1}\\in [0,N_1+2\\cdot{\\rm NGHOSTS})$\n",
    "* ${\\rm i2}\\in [0,N_2+2\\cdot{\\rm NGHOSTS})$,\n",
    "\n",
    "which guarantees the interior is covered by exactly $\\{N_0,N_1,N_2\\}$ grid points, the boundaries are covered by ${\\rm NGHOSTS}$ ghost zones, and we maintain cell centering."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So for example, if we choose a numerical grid in *spherical* coordinates $\\{r,\\theta,\\phi\\}$, with 3 ghost zone points (needed for e.g., 6th-order-accurate centered finite differencing), and we want the grid interior to be sampled with $\\{N_r,N_\\theta,N_\\phi\\}$ grid points, then we have\n",
    "\n",
    "\\begin{align}\n",
    "{\\rm NGHOSTS} &= 3 \\\\\n",
    "dr &= \\frac{r_{\\rm max} - 0}{N_r} \\\\\n",
    "d\\theta &= \\frac{\\pi - 0}{N_\\theta} \\\\\n",
    "d\\phi &= \\frac{\\pi - (-\\pi)}{N_\\phi} \\\\\n",
    "r_{{\\rm i0}} &= 0 + \\left[({\\rm i0}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_0 \\\\\n",
    "&= \\left[({\\rm i0}-3) + \\frac{1}{2}\\right] dx_0 \\\\\n",
    "\\theta_{{\\rm i1}} &= 0 + \\left[({\\rm i1}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_1 \\\\\n",
    "&= \\left[({\\rm i1}-3) + \\frac{1}{2}\\right] dx_1 \\\\\n",
    "\\phi_{{\\rm i2}} &= -\\pi + \\left[({\\rm i2}-{\\rm NGHOSTS}) + \\frac{1}{2}\\right] dx_2 \\\\\n",
    "&= -\\pi + \\left[({\\rm i2}-3) + \\frac{1}{2}\\right] dx_2, \\\\\n",
    "\\end{align}\n",
    "\n",
    "where again\n",
    "* ${\\rm i0}\\in [0,N_r+2\\cdot{\\rm NGHOSTS})$\n",
    "* ${\\rm i1}\\in [0,N_\\theta+2\\cdot{\\rm NGHOSTS})$\n",
    "* ${\\rm i2}\\in [0,N_\\phi+2\\cdot{\\rm NGHOSTS})$,\n",
    "\n",
    "which guarantees the interior is covered by exactly $\\{N_r,N_\\theta,N_\\phi\\}$ grid points, the boundaries are covered by ${\\rm NGHOSTS}$ ghost zones, and we maintain cell centering.\n",
    "\n",
    "Notice that in NRPy+, we use the [physics](https://en.wikipedia.org/wiki/Spherical_coordinate_system) notation for spherical coordinates, where $\\theta$ is the polar and $\\phi$ is the azimuthal angle. Also, we choose $\\phi$ to range from $-\\pi$ to $+\\pi$, which is most useful since it is compatible with output from [`atan2`](https://en.wikipedia.org/wiki/Atan2).\n",
    "\n",
    "**Exercise to student**: Given the prescription above, why do the integers $N_\\theta$ and $N_\\phi$ need to be even?\n",
    "\n",
    "As Laplacians like these appear on the right-hand sides of, e.g., the scalar wave equation in curvilinear coordinates, we still have a problem of some terms becoming quite large as the coordinate singularity is approached. This issue manifests as a stiffening of the PDE, requiring that we be very careful about the precise [Method of Lines](Tutorial-Method_of_Lines-C_Code_Generation.ipynb) timestepping algorithm used. See [Cordero-Carrión & Cerdá-Durán](https://arxiv.org/abs/1211.5930) for information on dealing with this subtlety in a second-order Runge-Kutta Method of Lines context; it was later found that the standard RK4 method maintains stable solutions to PDEs affected by this sort of stiffening.\n",
    "\n",
    "The above discussion focuses primarily on scalar fields. However, when solving PDEs involving vectors and tensors, the vectors and tensors themselves can exhibit divergent behavior at coordinate singularities. The good news is, this singular behavior is well-understood in terms of the scale factors of the reference metric, enabling us to define rescaled versions of these quantities that are well-behaved (so that, e.g., they can be finite differenced).\n",
    "\n",
    "For example, given a smooth vector *in a 3D Cartesian basis* $\\bar{\\Lambda}^{i}$, all components $\\bar{\\Lambda}^{x}$, $\\bar{\\Lambda}^{y}$, and $\\bar{\\Lambda}^{z}$ will be smooth (by assumption). When changing the basis to spherical coordinates (applying the appropriate Jacobian matrix transformation), we will find that since $\\phi = \\arctan(y/x)$, $\\bar{\\Lambda}^{\\phi}$ is given by\n",
    "\n",
    "\\begin{align}\n",
    "\\bar{\\Lambda}^{\\phi} &= \\frac{\\partial \\phi}{\\partial x} \\bar{\\Lambda}^{x} + \n",
    "\\frac{\\partial \\phi}{\\partial y} \\bar{\\Lambda}^{y} + \n",
    "\\frac{\\partial \\phi}{\\partial z} \\bar{\\Lambda}^{z} \\\\\n",
    "&= -\\frac{y}{x^2+y^2} \\bar{\\Lambda}^{x} + \n",
    "\\frac{x}{x^2+y^2} \\bar{\\Lambda}^{y} \\\\\n",
    "&= -\\frac{y}{(r \\sin\\theta)^2} \\bar{\\Lambda}^{x} + \n",
    "\\frac{x}{(r \\sin\\theta)^2} \\bar{\\Lambda}^{y} \\\\\n",
    "&= -\\frac{r \\sin\\theta \\sin\\phi}{(r \\sin\\theta)^2} \\bar{\\Lambda}^{x} + \n",
    "\\frac{r \\sin\\theta \\cos\\phi}{(r \\sin\\theta)^2} \\bar{\\Lambda}^{y}\\\\\n",
    "&= -\\frac{\\sin\\phi}{r \\sin\\theta} \\bar{\\Lambda}^{x} + \n",
    "\\frac{\\cos\\phi}{r \\sin\\theta} \\bar{\\Lambda}^{y}\\\\\n",
    "\\end{align}\n",
    "\n",
    "Thus $\\bar{\\Lambda}^{\\phi}$ diverges at all points where $r\\sin\\theta=0$ (or equivalently where $x=y=0$; i.e., the $z$-axis) due to the $\\frac{1}{r\\sin\\theta}$ that appear in the Jacobian transformation.\n",
    "\n",
    "This divergence might pose no problem on cell-centered grids that avoid $r \\sin\\theta=0$, except that the BSSN equations require that *first and second derivatives* of quantities like $\\bar{\\Lambda}^{\\phi}$ be computed. Usual strategies for numerical approximation of these derivatives (e.g., finite difference methods) will \"see\" these divergences and errors generally will not drop to zero with increased numerical sampling of the functions at points near where the functions diverge.\n",
    "\n",
    "However, notice that if we define $\\lambda^{\\phi}$ such that\n",
    "\n",
    "$$\\bar{\\Lambda}^{\\phi} = \\frac{1}{r\\sin\\theta} \\lambda^{\\phi},$$\n",
    "\n",
    "then $\\lambda^{\\phi}$ will be smooth and non-divergent as well. The strategy when computing derivatives of $\\bar{\\Lambda}^{\\phi}$, therefore, is to perform the product rule on the above expression, computing derivatives of the scale factors *analytically* (i.e., exactly using a computer algebra system like SymPy), and smooth terms like $\\lambda^{\\phi}$ with finite-difference derivatives.\n",
    "\n",
    "Avoiding such singularities can be generalized to arbitrary coordinate systems, so long as $\\lambda^i$ is defined as:\n",
    "\n",
    "$$\\bar{\\Lambda}^{i} = \\frac{\\lambda^i}{\\text{scalefactor[i]}} ,$$\n",
    "\n",
    "where scalefactor\\[i\\] is the $i$th scale factor in the given coordinate system. This idea can be extended to covariant (lowered-index) vectors and arbitrary tensors, as described in [the BSSN quantities tutorial notebook](Tutorial-BSSN_quantities.ipynb#rescaling_tensors).\n",
    "\n",
    "**In summary, Challenge 4 is addressed via a combination of cell-centered grids, tensor rescaling, and a stable Method of Lines time-stepping algorithm. This tutorial notebook will therefore focus on addressing Challenges 1 through 3, which, coincidentally, are addressed via an appropriate boundary condition algorithm.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#basic_algorithm): Overview of boundary condition algorithm in curvilinear coordinates\n",
    "    1. [Step 1.a](#challenge1): Addressing Challenge 1: Distinguishing inner from outer boundary points\n",
    "    1. [Step 1.b](#challenge2): Addressing Challenge 2: Eigen-Coordinate systems\n",
    "    1. [Step 1.c](#challenge3): Addressing Challenge 3: Applying parity conditions to arbitrary-rank tensors\n",
    "1. [Step 2](#ccode_bc_struct): `bc_struct` & friends: Data structures for storing boundary condition information\n",
    "1. [Step 3](#nrpycodegen): NRPy+-based C code generation for parity conditions\n",
    "    1. [Step 3.a](#dotproducts): Python function `parity_conditions_symbolic_dot_products()`: Set up C code for computing unit-vector dot products (=parity) for each of the 10 parity condition types\n",
    "    1. [Step 3.b](#set_parity_type): Populate `NRPy_basic_defines.h` with gridfunction parity types, which are set based on gridfunction name\n",
    "1. [Step 4](#bcstruct_set_up): `bcstruct_set_up()`: C function for setting up `bcstruct`\n",
    "    1. [Step 4.a](#bcstruct_set_up_eigencoord_mapping): `EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt()`: Perform $(x,y,z)\\to (x_0,x_1,x_2), \\left(i_0,i_1,i_2\\right)$ mappings using Eigen-Coordinate approach\n",
    "    1. [Step 4.b](#bcstruct_set_up_inner_bc_parity): `set_parity_for_inner_boundary_single_pt()`: Perform $(x,y,z)\\to \\{(x_0,x_1,x_2), \\left(i_0,i_1,i_2\\right)\\}$ mappings using Eigen-Coordinate approach\n",
    "    1. [Step 4.c](#bcstruct_set_up_add_to_cfunc): Add `bcstruct_set_up()` to NRPy+ C function dictionary\n",
    "1. [Step 5](#apply_bcs_inner_only): `apply_bcs_inner_only()`: Inner boundary condition C function\n",
    "1. [Step 6](#outer_bcs): Outer boundary condition C functions\n",
    "    1. [Step 6.a](#extrapolation_bcs) `\"extrapolation\"` outer boundary conditions: apply quadratic polynomial extrapolation\n",
    "    1. [Step 6.b](#radiation_bcs): `\"radiation\"` outer boundary conditions\n",
    "        1. [Step 6.b.i](#radiation_bcs_theory): Core *ansatz*, and implications\n",
    "        1. [Step 6.b.ii](#radiation_bcs_numerical): Numerical implementation\n",
    "        1. [Step 6.b.iii](#radiation_bcs_numerical_inv_jacobian): The $\\partial_r f$ term: Computing $\\frac{\\partial x^i}{\\partial r}$\n",
    "        1. [Step 6.b.iv](#radiation_bcs_numerical_partial_i_f): The $\\partial_r f$ term: Computing $\\partial_i f$ with arbitrary-offset finite-difference derivatives\n",
    "        1. [Step 6.b.v](#radiation_bcs_numerical_compute_partial_r_f): The $\\partial_r f$ term: Putting it all together in `compute_partial_r_f()`\n",
    "        1. [Step 6.b.vi](#radiation_bcs_evaluating_k): Evaluating the $k$ in the $k/r^3$ term\n",
    "        1. [Step 6.b.vii](#radiation_bcs_apply_bcs_radiation): `radiation_bcs_single_pt()`: Apply radiation BCs to a single point\n",
    "        1. [Step 6.b.viii](#apply_bcs_outerradiation_and_inner): `apply_bcs_outerradiation_and_inner()`: Apply radiation BCs at all outer boundary points and inner BCs at all inner boundary points\n",
    "1. [Step 7](#start2finish): `CurviBC_Playground.c`: Start-to-Finish C code module for testing & validating curvilinear boundary conditions\n",
    "    1. [Step 7.a](#register_gfs): Register gridfunctions of all 10 parity types in NRPy+\n",
    "    1. [Step 7.b](#validate): Set up test data for Curvilinear Boundary Conditions code validation\n",
    "    1. [Step 7.c](#mainc): `CurviBC_Playground`'s `main.c` Code\n",
    "    1. [Step 7.d](#curvibc_setupall): Add all CurviBC C codes to C function dictionary, and add CurviBC definitions to `NRPy_basic_defines.h`\n",
    "1. [Step 8](#add_cfunction_dicts): Add all C functions to `Cfunction_dict`, also output `NRPy_basic_defines.h` and `NRPy_function_prototypes.h`\n",
    "1. [Step 9](#senr_compare): Validation: Compare with original SENR results\n",
    "1. [Step 10](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='basic_algorithm'></a>\n",
    "\n",
    "# Step 1: Overview of boundary condition algorithm in curvilinear coordinates \\[Back to [top](#toc)\\]\n",
    "$$\\label{basic_algorithm}$$\n",
    "\n",
    "Here we **review** the basic algorithm for addressing Challenges [1](#innerouterbcs) [2](#coordinversion), and [3](#parity) discussed in the [Introduction](#intro) above. \n",
    "\n",
    "The algorithm itself is **implemented** as C code in Steps [**2**](#bc_struct) (data structures), [**3**](#set_up__bc_gz_map_and_parity_condns) (searching entire grid for inner and outer boundary points and setting parities), [**4**](#set_bcstruct) (setting data structures for quick and efficient implementation of outer boundaries), and [**5**](#apply_bcs_curvilinear) (function to apply inner & outer boundary conditions).\n",
    "\n",
    "<a id='challenge1'></a>\n",
    "\n",
    "## Step 1.a: Addressing Challenge 1: Distinguishing inner from outer boundary points \\[Back to [top](#toc)\\]\n",
    "$$\\label{challenge1}$$\n",
    "\n",
    "At each ghost zone grid point $\\mathbf{d}_{\\rm gz}=(x_0,x_1,x_2)$, we will do the following:\n",
    "\n",
    "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)$. \n",
    "    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.\n",
    "    1. Ghost zone points for which $\\mathbf{d}_{\\rm new} \\equiv \\mathbf{d}_{\\rm gz}$ are on the outer boundary of the grid, and standard outer boundary conditions should be applied.\n",
    "\n",
    "In detail, the algorithm is as follows:\n",
    "\n",
    "1. Convert the coordinate $(x_0,x_1,x_2)$ for the ghost zone point to Cartesian coordinates $\\left(x(x_0,x_1,x_2),y(x_0,x_1,x_2),z(x_0,x_1,x_2)\\right)$. For example, if we choose ordinary spherical coordinates $(x_0,x_1,x_2)=(r,\\theta,\\phi)$, then\n",
    "    + $x(r,\\theta,\\phi) = r \\sin(\\theta) \\cos(\\phi) = x_0 \\sin(x_1) \\cos(x_2)$\n",
    "    + $y(r,\\theta,\\phi) = r \\sin(\\theta) \\sin(\\phi) = x_0 \\sin(x_1) \\sin(x_2)$\n",
    "    + $z(r,\\theta,\\phi) = r \\cos(\\theta)            = x_0 \\cos(x_1)$\n",
    "1. Once we have $(x,y,z)$, we then find the corresponding value $(x_0,x_1,x_2)_{\\rm in/OB}=(r,\\theta,\\phi)_{\\rm in/OB}$ *in the grid interior or outer boundary*, via the simple inverse formula:\n",
    "    + $r_{\\rm in/OB}      = x_{0, \\rm in/OB} = \\sqrt{x^2+y^2+z^2} \\in [0,\\infty)$\n",
    "    + $\\theta_{\\rm in/OB} = x_{1, \\rm in/OB} = {\\rm acos}\\left(\\frac{z}{\\sqrt{x^2+y^2+z^2}}\\right) \\in [0,\\pi]$\n",
    "    + $\\phi_{\\rm in/OB}   = x_{2, \\rm in/OB} = {\\rm atan2}(y,x) \\in [-\\pi,\\pi]$ [Wikipedia article on atan2()](https://en.wikipedia.org/w/index.php?title=Atan2&oldid=859313982)\n",
    "\n",
    "1. If $(x_0,x_1,x_2)_{\\rm in/OB}$ is the same as the original $(x_0,x_1,x_2)$, then we know $(x_0,x_1,x_2)$ is an outer boundary point (in spherical coordinates, at $r>{\\rm RMAX}$), and we store `(i0,i1,i2)`$_{\\rm in/OB} = (-1,-1,-1)$. Otherwise, we know that $(x_0,x_1,x_2)$ maps to some interior point at index `(i0,i1,i2)`, which we store:\n",
    "    + $\\rm{i0}_{\\rm in/OB}=\\frac{r_{\\rm in/OB}      -      r_{\\rm min}}{\\Delta r}      - \\frac{1}{2}$\n",
    "    + $\\rm{i1}_{\\rm in/OB}=\\frac{\\theta_{\\rm in/OB} - \\theta_{\\rm min}}{\\Delta \\theta} - \\frac{1}{2}$\n",
    "    + $\\rm{i2}_{\\rm in/OB}=\\frac{\\phi_{\\rm in/OB}   -   \\phi_{\\rm min}}{\\Delta \\phi}   - \\frac{1}{2}$\n",
    "\n",
    "1. When updating a ghost zone point `(i0,i1,i2)` in the domain exterior, if the corresponding `(i0,i1,i2)`$_{\\rm in/OB}$ was set to $(-1,-1,-1)$, then we apply outer boundary conditions. Otherwise, we simply copy the data from the interior point at `(i0,i1,i2)`$_{\\rm in/OB}$ to `(i0,i1,i2)`.\n",
    "\n",
    "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.\n",
    "\n",
    "<a id='challenge2'></a>\n",
    "\n",
    "## Step 1.b: Addressing Challenge 2: Eigen-Coordinate Systems \\[Back to [top](#toc)\\]\n",
    "$$\\label{challenge2}$$\n",
    "\n",
    "Suppose we were to rewrite the spherical coordinate $r$ as an arbitrary odd function of $x_0$ instead of $x_0$ itself. In that case, $r(-x_0)=-r(x_0)$, and all parity conditions remain unchanged. However, the inverse function, $x_0(r)$, may not be writable as a closed-form expression, requiring a Newton-Raphson root finder to find the appropriate boundary mappings. \n",
    "\n",
    "To greatly simplify the algorithm in the case of arbitrary $r(x_0)$ in Spherical-like coordinates, or $\\rho(x_0)$ or $z(x_2)$ in Cylindrical-like coordinates, we note that the coordinate mappings *and* parities for all Spherical-like coordinate systems are identical to the mappings and parities for ordinary Spherical coordinates. The same holds true for Cylindrical-like and SymTP-like coordinate systems. Thus so long as we know the correct \"Eigen-Coordinate system\" (i.e., Spherical in the case of SinhSpherical or SinhSphericalv2; Cylindrical in the case of SinhCylindrical; SymTP in the case of SinhSymTP; etc.) there is no need for a Newton-Raphson root finder to set up the boundary conditions.\n",
    "\n",
    "<a id='challenge3'></a>\n",
    "\n",
    "## Step 1.c: Addressing Challenge 3: Applying parity conditions to arbitrary-rank tensors \\[Back to [top](#toc)\\]\n",
    "$$\\label{challenge3}$$\n",
    "\n",
    "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.\n",
    "\n",
    "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$).\n",
    "\n",
    "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\"):\n",
    "$$\n",
    "P_1 = \\mathbf{e}^i\\left(\\mathbf{d}_{\\rm gz}\\right) \\cdot \\mathbf{e}^i\\left(\\mathbf{d}_{\\rm new}\\right).\n",
    "$$\n",
    "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}$\n",
    "1. Repeat the above for the remaining components of the rank-$N$ tensor $j\\in \\{2,3,...,N\\}$, storing each $P_j$.\n",
    "1. The tensor mapping from $\\mathbf{d}_{\\rm gz}$ to $\\mathbf{d}_{\\rm new}$ for this tensor $T^{ijk...}_{mnp...}$ will be given by\n",
    "$$\n",
    "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}.\n",
    "$$\n",
    "\n",
    "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 \n",
    "+ 1 parity condition (the trivial +1) for scalars (rank-0 tensors)\n",
    "+ 3 parity conditions for all rank-1 tensors (corresponding to each direction)\n",
    "+ 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)\n",
    "\n",
    "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. The following Table outlines the correct conditions for each:\n",
    "\n",
    "The appropriate dot products determining parity condition are assigned to each gridfunction based on the following numbering:\n",
    "\n",
    "Tensor type | Parity type | Dot product(s) determining parity condition (see equation above)\n",
    "--- | --- | ---\n",
    "Scalar (Rank-0 tensor) | 0 | (*none*)\n",
    "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)$\n",
    "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)$\n",
    "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)$\n",
    "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$\n",
    "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]$\n",
    "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]$\n",
    "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$\n",
    "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]$\n",
    "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$\n",
    "\n",
    "\n",
    "In the following few steps, we will document the data structures for and implementation of this boundary condition algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ccode_bc_struct'></a>\n",
    "\n",
    "# Step 2:  `bc_struct` & friends: Data structures for storing boundary condition information \\[Back to [top](#toc)\\]\n",
    "$$\\label{ccode_bc_struct}$$\n",
    "\n",
    "Here we define `bc_struct`, a C data structure that stores all information needed to apply boundary conditions at all boundary points. \n",
    "\n",
    "Information needed to fill each ghost zone is based on whether the ghost zone is an inner or an outer boundary point.\n",
    "\n",
    "**`innerpt_bc_struct`: C struct for inner boundary points**\n",
    "\n",
    "If the ghost zone point is an inner boundary point, it maps to *a different point* either on the grid outer boundary or the grid interior. Further, vectors or tensors will flip sign based on their parity across this boundary (see section above). To fill in each inner boundary point we store data to the C struct `inner_bc_struct`:\n",
    "\n",
    "```C\n",
    "typedef struct __innerpt_bc_struct__ {\n",
    "  int dstpt;  // dstpt is the 3D grid index IDX3S(i0,i1,i2) of the inner boundary point (i0,i1,i2)\n",
    "  int srcpt;  // srcpt is the 3D grid index (a la IDX3S) to which the inner boundary point maps\n",
    "  int8_t parity[10];  // parity[10] is a calculation of dot products for the 10 independent parity types\n",
    "} innerpt_bc_struct;\n",
    "```\n",
    "\n",
    "**`outerpt_bc_struct`: C struct for outer boundary points**\n",
    "\n",
    "Unlike inner boundary points, ghost zone points that are outer boundary points map to *themselves*, and are filled using outer boundary conditions. Further, outer boundary ghost zones must be filled in from the inside outward, as e.g., the outermost ghost zones may depend on ghost zones closer to the grid interior being set. Thus filling in outer boundary points generally requires data at neighboring gridpoints *in the direction of the grid interior*. For example in Cartesian coordinates, filling in the $x=x_{\\rm max}$ face of the grid depends on data at $x<x_{\\rm max}$. Similarly, filling in the $y=y_{\\rm min}$ face will require data at $y>y_{\\rm min}$.\n",
    "\n",
    "So some indication of not only the gridpoint `(i0,i1,i2)` but also the face on which it exists, must be stored. We set the latter to the 1-byte integers `FACEX0`, `FACEX1`, and `FACEX2`, such that if we are on the `x0=x0_max` face, `FACEX0=-1`, `FACEX1=FACEX2=0`, to indicate that `x0=x0_max` points will be filled in using data at `x0<x0_max`. Thus to fill in each outer boundary point we store data to the C struct `outer_bc_struct`:\n",
    "\n",
    "```C\n",
    "typedef struct __outerpt_bc_struct__ {\n",
    "  short i0,i1,i2;  // the outer boundary point grid index (i0,i1,i2), on the 3D grid\n",
    "  int8_t FACEX0,FACEX1,FACEX2;  // 1-byte integers that store\n",
    "  //                               FACEX0,FACEX1,FACEX2 = +1, 0, 0 if on the i0=i0min face,\n",
    "  //                               FACEX0,FACEX1,FACEX2 = -1, 0, 0 if on the i0=i0max face,\n",
    "  //                               FACEX0,FACEX1,FACEX2 =  0,+1, 0 if on the i1=i2min face,\n",
    "  //                               FACEX0,FACEX1,FACEX2 =  0,-1, 0 if on the i1=i1max face,\n",
    "  //                               FACEX0,FACEX1,FACEX2 =  0, 0,+1 if on the i2=i2min face, or\n",
    "  //                               FACEX0,FACEX1,FACEX2 =  0, 0,-1 if on the i2=i2max face,\n",
    "} outerpt_bc_struct;\n",
    "```\n",
    "\n",
    "**`bc_info_struct`: C struct for other boundary information**\n",
    "\n",
    "The `bc_info_struct` simply stores the total number of inner and outer boundary points. To be scalable when updating boundary conditions on multiple cores, as many boundary points are updated simultaneously as possible. Inner boundary points can all be updated simultaneously, while only one ghostzone level of min and max faces of outer boundary points can be updated simultaneously.\n",
    "\n",
    "```C\n",
    "typedef struct __bc_info_struct__ {\n",
    "  int num_inner_boundary_points;  // stores total number of inner boundary points\n",
    "  int num_pure_outer_boundary_points[NGHOSTS][3];  // stores number of outer boundary points on each\n",
    "  //                                                  ghostzone level and direction (update min and\n",
    "  //                                                  max faces simultaneously on multiple cores)\n",
    "  int bc_loop_bounds[NGHOSTS][6][6];  // stores outer boundary loop bounds. Unused after bcstruct_set_up()\n",
    "} bc_info_struct;\n",
    "```\n",
    "\n",
    "**`bc_struct`: A struct of structs for storing all boundary condition information**\n",
    "\n",
    "As described above, \n",
    "\n",
    "```C\n",
    "typedef struct __bc_struct__ {\n",
    "  innerpt_bc_struct *restrict inner_bc_array;  // information needed for updating each inner boundary point\n",
    "  outerpt_bc_struct *restrict pure_outer_bc_array[NGHOSTS*3]; // information needed for updating each outer\n",
    "  //                                                             boundary point\n",
    "  bc_info_struct bc_info;  // stores number of inner and outer boundary points, needed for setting loop\n",
    "  //                          bounds and parallelizing over as many boundary points as possible.\n",
    "} bc_struct;\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:14.965489Z",
     "iopub.status.busy": "2021-09-08T15:47:14.964975Z",
     "iopub.status.idle": "2021-09-08T15:47:15.183581Z",
     "shell.execute_reply": "2021-09-08T15:47:15.183117Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 0: Add NRPy's directory to the path\n",
    "# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n",
    "import os,sys\n",
    "\n",
    "# Step P1: Import needed NRPy+ core modules:\n",
    "from outputC import outputC      # NRPy+: Core C code output module\n",
    "from outputC import outC_NRPy_basic_defines_h_dict  # NRPy+: Core C code output module\n",
    "from outputC import add_to_Cfunction_dict  # NRPy+: Core C code output module\n",
    "from outputC import Cfunction    # NRPy+: Core C code output module\n",
    "from outputC import indent_Ccode # NRPy+: Core C code output module\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import grid as gri               # NRPy+: Functions having to do with numerical grids\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface\n",
    "import finite_difference as fin  # NRPy+: Finite-difference module\n",
    "import shutil, os, sys           # Standard Python modules for multiplatform OS-level functions\n",
    "from UnitTesting.assert_equal import check_zero  # NRPy+: Checks whether an expression evaluates to zero.\n",
    "\n",
    "# Step P2: Create C code output directory:\n",
    "Ccodesrootdir = os.path.join(\"CurviBoundaryConditions_Ccodes/\")\n",
    "# P2.a: First remove C code output directory if it exists\n",
    "# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
    "shutil.rmtree(Ccodesrootdir, ignore_errors=True)\n",
    "# P2.b: Then create a fresh directory\n",
    "# Then create a fresh directory\n",
    "cmd.mkdir(Ccodesrootdir)\n",
    "\n",
    "# CoordSystem = \"Cylindrical\"\n",
    "CoordSystem = \"SinhSpherical\"\n",
    "outer_bcs_type = \"extrapolation\"\n",
    "# outer_bcs_type = \"radiation\"\n",
    "RADIATION_BC_FD_ORDER = 4  # Set to -1 to adopt the same value as finite_difference::FD_CENTDERIVS_ORDER\n",
    "\n",
    "# Set the finite differencing order to 4; although this module doesn't compute\n",
    "#   finite difference derivatives, it does need to set the number of ghost zone\n",
    "#   cells, which is generally based on NGHOSTS, which itself depends\n",
    "#   on finite_difference::FD_CENTDERIVS_ORDER.\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", 4)\n",
    "\n",
    "_unused = par.Cparameters(\"int\", __name__, \"outer_bc_type\", \"EXTRAPOLATION_OUTER_BCS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define/declare core data structures for curvilinear BCs within `NRPy_basic_defines.h`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:15.187904Z",
     "iopub.status.busy": "2021-09-08T15:47:15.187482Z",
     "iopub.status.idle": "2021-09-08T15:47:15.189166Z",
     "shell.execute_reply": "2021-09-08T15:47:15.188788Z"
    }
   },
   "outputs": [],
   "source": [
    "# First register basic C data structures/macros inside NRPy_basic_defines.h\n",
    "def NRPy_basic_defines_CurviBC_data_structures():\n",
    "    return r\"\"\"\n",
    "// NRPy+ Curvilinear Boundary Conditions: Core data structures\n",
    "// Documented in: Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb\n",
    "\n",
    "#define EXTRAPOLATION_OUTER_BCS 0  // used to identify/specify params.outer_bc_type\n",
    "#define RADIATION_OUTER_BCS     1  // used to identify/specify params.outer_bc_type\n",
    "\n",
    "typedef struct __innerpt_bc_struct__ {\n",
    "  int dstpt;  // dstpt is the 3D grid index IDX3S(i0,i1,i2) of the inner boundary point (i0,i1,i2)\n",
    "  int srcpt;  // srcpt is the 3D grid index (a la IDX3S) to which the inner boundary point maps\n",
    "  int8_t parity[10];  // parity[10] is a calculation of dot products for the 10 independent parity types\n",
    "} innerpt_bc_struct;\n",
    "\n",
    "typedef struct __outerpt_bc_struct__ {\n",
    "  short i0,i1,i2;  // the outer boundary point grid index (i0,i1,i2), on the 3D grid\n",
    "  int8_t FACEX0,FACEX1,FACEX2;  // 1-byte integers that store\n",
    "  //                               FACEX0,FACEX1,FACEX2 = +1, 0, 0 if on the i0=i0min face,\n",
    "  //                               FACEX0,FACEX1,FACEX2 = -1, 0, 0 if on the i0=i0max face,\n",
    "  //                               FACEX0,FACEX1,FACEX2 =  0,+1, 0 if on the i1=i2min face,\n",
    "  //                               FACEX0,FACEX1,FACEX2 =  0,-1, 0 if on the i1=i1max face,\n",
    "  //                               FACEX0,FACEX1,FACEX2 =  0, 0,+1 if on the i2=i2min face, or\n",
    "  //                               FACEX0,FACEX1,FACEX2 =  0, 0,-1 if on the i2=i2max face,\n",
    "} outerpt_bc_struct;\n",
    "\n",
    "typedef struct __bc_info_struct__ {\n",
    "  int num_inner_boundary_points;  // stores total number of inner boundary points\n",
    "  int num_pure_outer_boundary_points[NGHOSTS][3];  // stores number of outer boundary points on each\n",
    "  //                                                  ghostzone level and direction (update min and\n",
    "  //                                                  max faces simultaneously on multiple cores)\n",
    "  int bc_loop_bounds[NGHOSTS][6][6];  // stores outer boundary loop bounds. Unused after bcstruct_set_up()\n",
    "} bc_info_struct;\n",
    "\n",
    "typedef struct __bc_struct__ {\n",
    "  innerpt_bc_struct *restrict inner_bc_array;  // information needed for updating each inner boundary point\n",
    "  outerpt_bc_struct *restrict pure_outer_bc_array[NGHOSTS*3]; // information needed for updating each outer\n",
    "  //                                                             boundary point\n",
    "  bc_info_struct bc_info;  // stores number of inner and outer boundary points, needed for setting loop\n",
    "  //                          bounds and parallelizing over as many boundary points as possible.\n",
    "} bc_struct;\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpycodegen'></a>\n",
    "\n",
    "#  Step 3: NRPy+-based C code generation for parity conditions \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpycodegen}$$\n",
    "\n",
    "Much of the algorithm needed for setting up `bcstruct` requires a loop over all gridpoints on the numerical grid. As the precise numerical grids are chosen at C runtime, that part of the algorithm must be run entirely within a static C code.\n",
    "\n",
    "However, there are two parts to the overall algorithm that must be generated by NRPy+, namely\n",
    "\n",
    "1. [Step 3.a](#dotproducts): `parity_conditions_symbolic_dot_products()`: Based on the chosen reference metric, sets up the needed unit-vector dot products for each of the 10 parity condition types.\n",
    "1. [Step 3.b](#set_parity_type): Set parity type for each gridfunction registered within NRPy+, based on the digits at the end of each gridfunction name, append result to `dirname+gridfunction_defines.h`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='dotproducts'></a>\n",
    "\n",
    "## Step 3.a: Python function `parity_conditions_symbolic_dot_products()`: Set up C code for computing unit-vector dot products (=parity) for each of the 10 parity condition types \\[Back to [top](#toc)\\]\n",
    "$$\\label{dotproducts}$$\n",
    "\n",
    "Next, we generate the C code necessary to perform needed dot products for filling in the parity condition arrays inside `bcstruct`. \n",
    "\n",
    "Using the unit vectors defined in `rfm.UnitVectors[][]` (in `reference_metric.py`), each unit vector takes as input either $\\mathbf{d}_{\\rm gz} = (x_0,x_1,x_2)_{\\rm IB}$=`(xx0,xx1,xx2)` or $\\mathbf{d}_{\\rm new} = (x_0,x_1,x_2)_{\\rm in}$=`(xx0_inbounds,xx1_inbounds,xx2_inbounds)` as summarized in the table above in [Step 1](#challenge2). We paste the table here again, for quick reference. Pay special attention to the parity type numbering, as this is the convention adopted in the code:\n",
    "\n",
    "Tensor type | Parity type | Dot product(s) determining parity condition (\n",
    "--- | --- | ---\n",
    "Scalar (Rank-0 tensor) | 0 | (*none*)\n",
    "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)$\n",
    "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)$\n",
    "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)$\n",
    "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$\n",
    "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]$\n",
    "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]$\n",
    "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$\n",
    "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]$\n",
    "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$\n",
    "\n",
    "Looping over all 10 parity types, the corresponding symbolic expressions for dot product(s) is output to C code. For example, in Spherical coordinates, parity type 1's dot product output as C code is given by: \n",
    "\n",
    "```\n",
    "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);\n",
    "```\n",
    "\n",
    "To wit (as described above), there are 10 parity types for BSSN evolved variables, which include scalars, vectors, and symmetric rank-2 tensors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:15.196318Z",
     "iopub.status.busy": "2021-09-08T15:47:15.195711Z",
     "iopub.status.idle": "2021-09-08T15:47:15.198183Z",
     "shell.execute_reply": "2021-09-08T15:47:15.197651Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set unit-vector dot products (=parity) for each of the 10 parity condition types\n",
    "def parity_conditions_symbolic_dot_products():\n",
    "    parity = ixp.zerorank1(DIM=10)\n",
    "    UnitVectors_inner = ixp.zerorank2()\n",
    "    xx0_inbounds,xx1_inbounds,xx2_inbounds = sp.symbols(\"xx0_inbounds xx1_inbounds xx2_inbounds\", real=True)\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            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)\n",
    "    # Type 0: scalar\n",
    "    parity[0] = sp.sympify(1)\n",
    "    # Type 1: i0-direction vector or one-form\n",
    "    # Type 2: i1-direction vector or one-form\n",
    "    # Type 3: i2-direction vector or one-form\n",
    "    for i in range(3):\n",
    "        for Type in range(1,4):\n",
    "            parity[Type] += rfm.UnitVectors[Type-1][i]*UnitVectors_inner[Type-1][i]\n",
    "    # Type 4: i0i0-direction rank-2 tensor\n",
    "    # parity[4] = parity[1]*parity[1]\n",
    "    # Type 5: i0i1-direction rank-2 tensor\n",
    "    # Type 6: i0i2-direction rank-2 tensor\n",
    "    # Type 7: i1i1-direction rank-2 tensor\n",
    "    # Type 8: i1i2-direction rank-2 tensor\n",
    "    # Type 9: i2i2-direction rank-2 tensor\n",
    "    count = 4\n",
    "    for i in range(3):\n",
    "        for j in range(i,3):\n",
    "            parity[count] = parity[i+1]*parity[j+1]\n",
    "            count = count + 1\n",
    "\n",
    "    lhs_strings = []\n",
    "    for i in range(10):\n",
    "        lhs_strings.append(\"REAL_parity_array[\"+str(i)+\"]\")\n",
    "    outstr = \"\"\"\n",
    "    // NRPy+ Curvilinear Boundary Conditions: Unit vector dot products for all\n",
    "    //      ten parity conditions, in given coordinate system.\n",
    "    //      Needed for automatically determining sign of tensor across coordinate boundary.\n",
    "    // Documented in: Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb\n",
    "\"\"\"\n",
    "    return outstr + outputC(parity, lhs_strings, filename=\"returnstring\", params=\"preindent=2\")\n",
    "\n",
    "#     print(\"\\n\\nExample: parity type 1's dot product is given by: \\n\"+lhs_strings[1]+\" = \"+str(parity[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='set_parity_type'></a>\n",
    "\n",
    "## Step 3.b: Populate `NRPy_basic_defines.h` with gridfunction parity types, which are set based on gridfunction name \\[Back to [top](#toc)\\]\n",
    "$$\\label{set_parity_type}$$\n",
    "\n",
    "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 rather robust algorithm, because `gri.register_gridfunctions()` in [grid.py](../edit/grid.py) will throw an error if a gridfunction's **base** name ends in an integer. This strict syntax was added with the express purpose of making it easier to set parity types based solely on the gridfunction name.\n",
    "\n",
    "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 `NRPy_basic_defines.h`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:15.214520Z",
     "iopub.status.busy": "2021-09-08T15:47:15.211805Z",
     "iopub.status.idle": "2021-09-08T15:47:15.215728Z",
     "shell.execute_reply": "2021-09-08T15:47:15.216134Z"
    }
   },
   "outputs": [],
   "source": [
    "def NRPy_basic_defines_set_gridfunction_defines_with_parity_types(verbose=True):\n",
    "    # First add human-readable gridfunction aliases (grid.py) to NRPy_basic_defines dictionary,\n",
    "    evolved_variables_list, auxiliary_variables_list, auxevol_variables_list = gri.gridfunction_lists()[0:3]\n",
    "\n",
    "    # Step 3.b: set the parity conditions on all gridfunctions in gf_list,\n",
    "    #       based on how many digits are at the end of their names\n",
    "    def set_parity_types(list_of_gf_names):\n",
    "        parity_type = []\n",
    "        for name in list_of_gf_names:\n",
    "            for gf in gri.glb_gridfcs_list:\n",
    "                if gf.name == name:\n",
    "                    parity_type__orig_len = len(parity_type)\n",
    "                    if gf.DIM < 3 or gf.DIM > 4:\n",
    "                        print(\"Error: Cannot currently specify parity conditions on gridfunctions with DIM<3 or >4.\")\n",
    "                        sys.exit(1)\n",
    "                    if gf.rank == 0:\n",
    "                        parity_type.append(0)\n",
    "                    elif gf.rank == 1:\n",
    "                        if gf.DIM == 3:\n",
    "                            parity_type.append(int(gf.name[-1]) + 1)  # = 1 for e.g., beta^0; = 2 for e.g., beta^1, etc.\n",
    "                        elif gf.DIM == 4:\n",
    "                            parity_type.append(int(gf.name[-1]))  # = 0 for e.g., b4^0; = 1 for e.g., beta^1, etc.\n",
    "                    elif gf.rank == 2:\n",
    "                        if gf.DIM == 3:\n",
    "                            # element of a list; a[-2] the\n",
    "                            # second-to-last element, etc.\n",
    "                            idx0 = gf.name[-2]\n",
    "                            idx1 = gf.name[-1]\n",
    "                            if idx0 == \"0\" and idx1 == \"0\":\n",
    "                                parity_type.append(4)\n",
    "                            elif (idx0 == \"0\" and idx1 == \"1\") or (idx0 == \"1\" and idx1 == \"0\"):\n",
    "                                parity_type.append(5)\n",
    "                            elif (idx0 == \"0\" and idx1 == \"2\") or (idx0 == \"2\" and idx1 == \"0\"):\n",
    "                                parity_type.append(6)\n",
    "                            elif idx0 == \"1\" and idx1 == \"1\":\n",
    "                                parity_type.append(7)\n",
    "                            elif (idx0 == \"1\" and idx1 == \"2\") or (idx0 == \"2\" and idx1 == \"1\"):\n",
    "                                parity_type.append(8)\n",
    "                            elif idx0 == \"2\" and idx1 == \"2\":\n",
    "                                parity_type.append(9)\n",
    "                        elif gf.DIM == 4:\n",
    "                            idx0 = gf.name[-2]\n",
    "                            idx1 = gf.name[-1]\n",
    "                            # g4DD00 = g_{tt} : parity type = 0\n",
    "                            # g4DD01 = g_{tx} : parity type = 1\n",
    "                            # g4DD02 = g_{ty} : parity type = 2\n",
    "                            # g4DD0a = g_{ta} : parity type = a\n",
    "                            if idx0 == \"0\":\n",
    "                                parity_type.append(int(idx1))\n",
    "                            elif idx1 == \"0\":\n",
    "                                parity_type.append(int(idx0))\n",
    "                            if idx0 == \"1\" and idx1 == \"1\":\n",
    "                                parity_type.append(4)\n",
    "                            elif (idx0 == \"1\" and idx1 == \"2\") or (idx0 == \"2\" and idx1 == \"1\"):\n",
    "                                parity_type.append(5)\n",
    "                            elif (idx0 == \"1\" and idx1 == \"3\") or (idx0 == \"3\" and idx1 == \"1\"):\n",
    "                                parity_type.append(6)\n",
    "                            elif idx0 == \"2\" and idx1 == \"2\":\n",
    "                                parity_type.append(7)\n",
    "                            elif (idx0 == \"2\" and idx1 == \"3\") or (idx0 == \"3\" and idx1 == \"2\"):\n",
    "                                parity_type.append(8)\n",
    "                            elif idx0 == \"3\" and idx1 == \"3\":\n",
    "                                parity_type.append(9)\n",
    "                    if len(parity_type) == parity_type__orig_len:\n",
    "                        print(\"Error: Could not figure out parity type for \"+gf.gftype+\" gridfunction: \" + gf.name,gf.DIM,gf.name[-2],gf.name[-1],gf.rank)\n",
    "                        sys.exit(1)\n",
    "        if len(parity_type) != len(list_of_gf_names):\n",
    "            print(\"Error: For some reason the length of the parity types list did not match the length of the gf list.\")\n",
    "            sys.exit(1)\n",
    "        return parity_type\n",
    "\n",
    "    evol_parity_type = set_parity_types(evolved_variables_list)\n",
    "    aux_parity_type = set_parity_types(auxiliary_variables_list)\n",
    "    auxevol_parity_type = set_parity_types(auxevol_variables_list)\n",
    "\n",
    "    # Output all gridfunctions to Ccodesrootdir/gridfunction_defines.h\n",
    "    # ... then append to the file the parity type for each gridfunction.\n",
    "    outstr = \"\"\"\n",
    "/* PARITY TYPES FOR ALL GRIDFUNCTIONS.\n",
    " * SEE \\\"Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb\\\" FOR DEFINITIONS. */\n",
    "\"\"\"\n",
    "    if len(evolved_variables_list) > 0:\n",
    "        outstr += \"static const int8_t evol_gf_parity[\" + str(len(evolved_variables_list)) + \"] = { \"\n",
    "        for i in range(len(evolved_variables_list) - 1):\n",
    "            outstr += str(evol_parity_type[i]) + \", \"\n",
    "        outstr += str(evol_parity_type[len(evolved_variables_list) - 1]) + \" };\\n\"\n",
    "\n",
    "    if len(auxiliary_variables_list) > 0:\n",
    "        outstr += \"static const int8_t aux_gf_parity[\" + str(len(auxiliary_variables_list)) + \"] = { \"\n",
    "        for i in range(len(auxiliary_variables_list) - 1):\n",
    "            outstr += str(aux_parity_type[i]) + \", \"\n",
    "        outstr += str(aux_parity_type[len(auxiliary_variables_list) - 1]) + \" };\\n\"\n",
    "\n",
    "    if len(auxevol_variables_list) > 0:\n",
    "        outstr += \"static const int8_t auxevol_gf_parity[\" + str(len(auxevol_variables_list)) + \"] = { \"\n",
    "        for i in range(len(auxevol_variables_list) - 1):\n",
    "            outstr += str(auxevol_parity_type[i]) + \", \"\n",
    "        outstr += str(auxevol_parity_type[len(auxevol_variables_list) - 1]) + \" };\\n\"\n",
    "\n",
    "    if verbose == True:\n",
    "        for i in range(len(evolved_variables_list)):\n",
    "            print(\"Evolved gridfunction \\\"\" + evolved_variables_list[i] + \"\\\" has parity type \" + str(\n",
    "                evol_parity_type[i]) + \".\")\n",
    "        for i in range(len(auxiliary_variables_list)):\n",
    "            print(\"Auxiliary gridfunction \\\"\" + auxiliary_variables_list[i] + \"\\\" has parity type \" + str(\n",
    "                aux_parity_type[i]) + \".\")\n",
    "        for i in range(len(auxevol_variables_list)):\n",
    "            print(\"AuxEvol gridfunction \\\"\" + auxevol_variables_list[i] + \"\\\" has parity type \" + str(\n",
    "                auxevol_parity_type[i]) + \".\")\n",
    "    return outstr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bcstruct_set_up'></a>\n",
    "\n",
    "# Step 4: `bcstruct_set_up()`: C function for setting up `bcstruct` \\[Back to [top](#toc)\\]\n",
    "$$\\label{bcstruct_set_up}$$\n",
    "\n",
    "`bcstruct_set_up()` sets up `bcstruct`, the C struct that contains all information needed to fill in grid boundaries. It depends on two C helper functions, which are included in `bcstruct_set_up.c`:\n",
    "\n",
    "1. `EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds()` uses the [Eigen-Coordinate](#challenge2) approach described above to find the mapping of boundary point `(x0,x1,x2)` $\\to$ `(x,y,z)` $\\to$ `(x0',x1',x2')`. If `(x0,x1,x2)` is an inner boundary point, then `(x0',x1',x2')` $\\equiv$ `(x0,x1,x2)_inbounds` will map to a different gridpoint-either in the grid interior or at an outer boundary point. I.e., when `(x0,x1,x2)`$\\ne$`(x0',x1',x2')`, `(x0,x1,x2)` is an *inner* boundary gridpoint. Otherwise the boundary point `(x0,x1,x2)` is an *outer* boundary gridpoint.\n",
    "1. `set_parity_for_inner_boundary_single_pt()` sets the parity condition for inner boundary points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bcstruct_set_up_eigencoord_mapping'></a>\n",
    "\n",
    "## Step 4.a:  `EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt()`: Perform $(x,y,z)\\to (x_0,x_1,x_2), \\left(i_0,i_1,i_2\\right)$ mappings using Eigen-Coordinate approach \\[Back to [top](#toc)\\]\n",
    "$$\\label{bcstruct_set_up_eigencoord_mapping}$$\n",
    "\n",
    "\n",
    "First we generate the C code needed for applying boundary conditions in generic coordinate systems, using the [Eigen-Coordinate](#challenge2) approach described above:\n",
    "* $\\left(x(x_0,x_1,x_2),y(x_0,x_1,x_2),z(x_0,x_1,x_2)\\right) \\to \\left(x_0(x,y,z),x_1(x,y,z),x_2(x,y,z)\\right), \\left(i_0,i_1,i_2\\right)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Cfunction__EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt():\n",
    "    desc = \"\"\"EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt():\n",
    "  A coordinate system's \"eigencoordinate\" is the simplest member\n",
    "  of its family; all spherical-like coordinate systems have\n",
    "  Spherical as their eigencoordinate. The same is true for\n",
    "  cylindrical-like (Cylindrical is eigencoordinate),\n",
    "  Cartesian-like (Cartesian is the eigencoordinate), and\n",
    "  SymTP-like (SymTP is the eigencoordinate) coordinates.\n",
    "\n",
    "  For a given gridpoint (i0,i1,i2) and corresponding coordinate\n",
    "  (x0,x1,x2), this function performs the dual mapping\n",
    "  (x0,x1,x2) -> (Cartx,Carty,Cartz) -> (x0,x1,x2)'\n",
    "  Note that (x0,x1,x2) IS NOT ALWAYS equal to (x0,x1,x2)';\n",
    "  For example consider in Spherical coordinates\n",
    "  (x0,x1,x2)=(r,theta,phi)=(-0.1,pi/4,pi/4).\n",
    "  This point will map to (x0,x1,x2)', in which x0>0,\n",
    "  because the inversion r=sqrt(Cartx^2+Carty^2+Cartz^2)\n",
    "  is always positive. In this case, (x0,x1,x2) is considered\n",
    "  an *inner* boundary point, and on a cell-centered grid\n",
    "  is guaranteed to map to a grid point in the grid interior;\n",
    "  filling in this point requires copying data, and possibly\n",
    "  multiplying by a +/- 1 if the data is from a gridfunction\n",
    "  storing tensors/vectors.\n",
    "\"\"\"\n",
    "    c_type = \"static void\"\n",
    "    name = \"EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt\"\n",
    "    params = \"\"\"const paramstruct *restrict params, REAL *restrict xx[3],\n",
    "                                                               const int i0, const int i1, const int i2,\n",
    "                                                               REAL x0x1x2_inbounds[3], int i0i1i2_inbounds[3]\"\"\"\n",
    "    body = r\"\"\"\n",
    "  // This is a 3-step algorithm:\n",
    "  // Step 1: (x0,x1,x2) -> (Cartx,Carty,Cartz)\n",
    "  //         Find the Cartesian coordinate that (x0,x1,x2)\n",
    "  //         maps to, assuming (x0,x1,x2) is the eigen-\n",
    "  //         coordinate. Note that we assume (x0,x1,x2)\n",
    "  //         has the same grid boundaries in both the\n",
    "  //         original coordinate and the eigencoordinate.\n",
    "  // Step 2: (Cartx,Carty,Cartz) -> (x0,x1,x2)'\n",
    "  //         Find the interior eigencoordinate point\n",
    "  //         (x0,x1,x2)' to which (Cartx,Carty,Cartz)\n",
    "  //         maps, as well as the corresponding\n",
    "  //         gridpoint integer index (i0,i1,i2). For\n",
    "  //         cell-centered grids, (x0,x1,x2) will always\n",
    "  //         overlap exactly (to roundoff error) a point\n",
    "  //         on the numerical grid.\n",
    "  // Step 3: Sanity check\n",
    "  //         Convert x0(i0_inbounds),x1(i1_inbounds),x2(i2_inbounds) -> (Cartx,Carty,Cartz),\n",
    "  //         and check that\n",
    "  //         (Cartx,Carty,Cartz) == (Cartx(x0(i0)),Cartx(x1(i1)),Cartx(x2(i2)))\n",
    "  //         If not, error out!\n",
    "\"\"\"\n",
    "    # Load up the EigenCoordinate corresponding to reference_metric::CoordSystem\n",
    "    CoordSystem_orig = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\",rfm.get_EigenCoord())\n",
    "    rfm.reference_metric()\n",
    "\n",
    "    # Step 1: Output C code for the Eigen-Coordinate mapping from xx->Cartesian':\n",
    "    body += r\"\"\"\n",
    "  // Step 1: Convert the (curvilinear) coordinate (x0,x1,x2) to Cartesian coordinates\n",
    "  REAL xCart[3];  // where (x,y,z) is output\n",
    "  {\n",
    "    // xx_to_Cart for EigenCoordinate \"\"\"+rfm.get_EigenCoord()+r\"\"\" (orig coord = \"\"\"+CoordSystem_orig+r\"\"\"):\n",
    "    REAL xx0 = xx[0][i0];\n",
    "    REAL xx1 = xx[1][i1];\n",
    "    REAL xx2 = xx[2][i2];\n",
    "\"\"\"+ \\\n",
    "    outputC([rfm.xx_to_Cart[0],rfm.xx_to_Cart[1],rfm.xx_to_Cart[2]],\n",
    "            [\"xCart[0]\",\"xCart[1]\",\"xCart[2]\"],\n",
    "            \"returnstring\", params=\"preindent=2\")+\"  }\\n\"\n",
    "    body += r\"\"\"\n",
    "  REAL Cartx = xCart[0];\n",
    "  REAL Carty = xCart[1];\n",
    "  REAL Cartz = xCart[2];\"\"\"\n",
    "\n",
    "    # Step 2: Output C code for the Eigen-Coordinate mapping from Cartesian->xx':\n",
    "    body += r\"\"\"\n",
    "  // Step 2: Find the (i0_inbounds,i1_inbounds,i2_inbounds) corresponding to the above Cartesian coordinate.\n",
    "  //   If (i0_inbounds,i1_inbounds,i2_inbounds) is in a ghost zone, then it must equal (i0,i1,i2), and\n",
    "  //      the point is an outer boundary point.\n",
    "  //   Otherwise (i0_inbounds,i1_inbounds,i2_inbounds) is in the grid interior, and data at (i0,i1,i2)\n",
    "  //      must be replaced with data at (i0_inbounds,i1_inbounds,i2_inbounds), but multiplied by the\n",
    "  //      appropriate parity condition (+/- 1).\n",
    "  REAL Cart_to_xx0_inbounds,Cart_to_xx1_inbounds,Cart_to_xx2_inbounds;\n",
    "\"\"\"\n",
    "    # Step 2.a: Sanity check: First make sure that rfm.Cart_to_xx has been set. Error out if not!\n",
    "    if rfm.Cart_to_xx[0] == 0 or rfm.Cart_to_xx[1] == 0 or rfm.Cart_to_xx[2] == 0:\n",
    "        print(\"ERROR: rfm.Cart_to_xx[], which maps Cartesian -> xx, has not been set for\")\n",
    "        print(\"       reference_metric::CoordSystem = \"+par.parval_from_str(\"reference_metric::CoordSystem\"))\n",
    "        print(\"       Boundary conditions in curvilinear coordinates REQUiRE this be set.\")\n",
    "        sys.exit(1)\n",
    "    # Step 2.b: Output C code for the Eigen-Coordinate mapping from Cartesian->xx:\n",
    "    body += \"\"\"  // Cart_to_xx for EigenCoordinate \"\"\"+rfm.get_EigenCoord()+r\"\"\" (orig coord = \"\"\"+CoordSystem_orig+\");\\n\"\n",
    "    body += outputC([rfm.Cart_to_xx[0],rfm.Cart_to_xx[1],rfm.Cart_to_xx[2]],\n",
    "                    [\"Cart_to_xx0_inbounds\",\"Cart_to_xx1_inbounds\",\"Cart_to_xx2_inbounds\"],\n",
    "                    filename=\"returnstring\", params=\"preindent=1\")\n",
    "    body += r\"\"\"\n",
    "  // Next compute xxmin[i]. By definition,\n",
    "  //    xx[i][j] = xxmin[i] + ((REAL)(j-NGHOSTS) + (1.0/2.0))*dxxi;\n",
    "  // -> xxmin[i] = xx[i][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxxi\n",
    "  const REAL xxmin[3] = {\n",
    "    xx[0][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx0,\n",
    "    xx[1][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx1,\n",
    "    xx[2][0] - ((REAL)(0-NGHOSTS) + (1.0/2.0))*dxx2 };\n",
    "\n",
    "  // Finally compute i{0,1,2}_inbounds (add 0.5 to account for rounding down)\n",
    "  const int i0_inbounds = (int)( (Cart_to_xx0_inbounds - xxmin[0] - (1.0/2.0)*dxx0 + ((REAL)NGHOSTS)*dxx0)/dxx0 + 0.5 );\n",
    "  const int i1_inbounds = (int)( (Cart_to_xx1_inbounds - xxmin[1] - (1.0/2.0)*dxx1 + ((REAL)NGHOSTS)*dxx1)/dxx1 + 0.5 );\n",
    "  const int i2_inbounds = (int)( (Cart_to_xx2_inbounds - xxmin[2] - (1.0/2.0)*dxx2 + ((REAL)NGHOSTS)*dxx2)/dxx2 + 0.5 );\n",
    "\"\"\"\n",
    "\n",
    "    # Restore reference_metric::CoordSystem back to the original CoordSystem\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\",CoordSystem_orig)\n",
    "    rfm.reference_metric()\n",
    "\n",
    "    # Step 3:\n",
    "    body += \"\"\"\n",
    "  // Step 3: Convert x0(i0_inbounds),x1(i1_inbounds),x2(i2_inbounds) -> (Cartx,Carty,Cartz),\n",
    "  //         and check that\n",
    "  //         (Cartx,Carty,Cartz) == (Cartx(x0(i0)),Cartx(x1(i1)),Cartx(x2(i2)))\n",
    "  //         If not, error out!\n",
    "\n",
    "  // Step 3.a: Compute x0(i0_inbounds),x1(i1_inbounds),x2(i2_inbounds):\n",
    "  const REAL x0_inbounds = xx[0][i0_inbounds];\n",
    "  const REAL x1_inbounds = xx[1][i1_inbounds];\n",
    "  const REAL x2_inbounds = xx[2][i2_inbounds];\n",
    "\n",
    "  // Step 3.b: Compute {x,y,z}Cart_from_xx, as a\n",
    "  //           function of i0,i1,i2\n",
    "  REAL xCart_from_xx, yCart_from_xx, zCart_from_xx;\n",
    "  {\n",
    "    // xx_to_Cart for Coordinate \"\"\"+CoordSystem_orig+r\"\"\"):\n",
    "    REAL xx0 = xx[0][i0];\n",
    "    REAL xx1 = xx[1][i1];\n",
    "    REAL xx2 = xx[2][i2];\n",
    "\"\"\"+ \\\n",
    "    outputC([rfm.xx_to_Cart[0],rfm.xx_to_Cart[1],rfm.xx_to_Cart[2]],\n",
    "            [\"xCart_from_xx\",\"yCart_from_xx\",\"zCart_from_xx\"],\n",
    "            \"returnstring\", params=\"preindent=2,includebraces=False\")\n",
    "    body += r\"\"\"  }\n",
    "\n",
    "  // Step 3.c: Compute {x,y,z}Cart_from_xx_inbounds, as a\n",
    "  //           function of i0_inbounds,i1_inbounds,i2_inbounds\n",
    "  REAL xCart_from_xx_inbounds, yCart_from_xx_inbounds, zCart_from_xx_inbounds;\n",
    "  {\n",
    "    // xx_to_Cart_inbounds for Coordinate \"\"\"+CoordSystem_orig+r\"\"\"):\n",
    "    REAL xx0 = xx[0][i0_inbounds];\n",
    "    REAL xx1 = xx[1][i1_inbounds];\n",
    "    REAL xx2 = xx[2][i2_inbounds];\n",
    "\"\"\"+ \\\n",
    "    outputC([rfm.xx_to_Cart[0],rfm.xx_to_Cart[1],rfm.xx_to_Cart[2]],\n",
    "            [\"xCart_from_xx_inbounds\",\"yCart_from_xx_inbounds\",\"zCart_from_xx_inbounds\"],\n",
    "            \"returnstring\", params=\"preindent=2,includebraces=False\")\n",
    "    body += r\"\"\"  }\n",
    "\n",
    "  // Step 3.d: Compare xCart_from_xx to xCart_from_xx_inbounds;\n",
    "  //           they should be identical!!!\n",
    "#define EPS_REL 1e-8\n",
    "  const REAL norm_factor = sqrt(xCart_from_xx*xCart_from_xx + yCart_from_xx*yCart_from_xx + zCart_from_xx*zCart_from_xx) + 1e-15;\n",
    "  if(fabs( (double)(xCart_from_xx - xCart_from_xx_inbounds) ) > EPS_REL * norm_factor ||\n",
    "     fabs( (double)(yCart_from_xx - yCart_from_xx_inbounds) ) > EPS_REL * norm_factor ||\n",
    "     fabs( (double)(zCart_from_xx - zCart_from_xx_inbounds) ) > EPS_REL * norm_factor) {\n",
    "    fprintf(stderr,\"Error in \"\"\"+CoordSystem_orig+r\"\"\" coordinate system: Cartesian disagreement: ( %.15e %.15e %.15e ) != ( %.15e %.15e %.15e ) | xx: %e %e %e -> %e %e %e | %d %d %d\\n\",\n",
    "            (double)xCart_from_xx,(double)yCart_from_xx,(double)zCart_from_xx,\n",
    "            (double)xCart_from_xx_inbounds,(double)yCart_from_xx_inbounds,(double)zCart_from_xx_inbounds,\n",
    "            xx[0][i0],xx[1][i1],xx[2][i2],\n",
    "            xx[0][i0_inbounds],xx[1][i1_inbounds],xx[2][i2_inbounds],\n",
    "            Nxx_plus_2NGHOSTS0,Nxx_plus_2NGHOSTS1,Nxx_plus_2NGHOSTS2);\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  // Step 4: Set output arrays.\n",
    "  x0x1x2_inbounds[0] = xx[0][i0_inbounds];\n",
    "  x0x1x2_inbounds[1] = xx[1][i1_inbounds];\n",
    "  x0x1x2_inbounds[2] = xx[2][i2_inbounds];\n",
    "  i0i1i2_inbounds[0] = i0_inbounds;\n",
    "  i0i1i2_inbounds[1] = i1_inbounds;\n",
    "  i0i1i2_inbounds[2] = i2_inbounds;\n",
    "\"\"\"\n",
    "    __function_prototype_ignore, Cfunc = Cfunction(\n",
    "        desc=desc, c_type=c_type, name=name, params=params,\n",
    "        body=body)\n",
    "    return Cfunc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bcstruct_set_up_inner_bc_parity'></a>\n",
    "\n",
    "## Step 4.b:  `set_parity_for_inner_boundary_single_pt()`: Perform $(x,y,z)\\to \\{(x_0,x_1,x_2), \\left(i_0,i_1,i_2\\right)\\}$ mappings using Eigen-Coordinate approach \\[Back to [top](#toc)\\]\n",
    "$$\\label{bcstruct_set_up_inner_bc_parity}$$\n",
    "\n",
    "Next, we define the C function `set_parity_for_inner_boundary_single_pt()`, which implements the algorithm described above in [Challenge 3](#challenge3).\n",
    "\n",
    "Given the inputs: \n",
    "* inner boundary point $(x_0,x_1,x_2)$ situated at grid index `(i0,i1,i2)`, and\n",
    "* the interior (non-boundary) point to which it maps $(x0,x1,x2)'=$ `(x0x1x2_inbounds[0],x0x1x2_inbounds[1],x0x1x2_inbounds[2])`\n",
    "\n",
    "this function evaluates dot products of the unit vectors evaluated at `(i0_inbounds,i1_inbounds,i2_inbounds)` and `(i0,i1,i2)`, for all 10 gridfunction parities currently supported in NRPy+. The C code for computing the needed symbolic dot products is generated above in [Step 3](#dotproducts)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Cfunction__set_parity_for_inner_boundary_single_pt():\n",
    "    desc = \"\"\"set_parity_for_inner_boundary_single_pt():\n",
    "  Given (x0,x1,x2)=(xx0,xx1,xx2) and\n",
    "  (x0,x1,x2)'=(x0x1x2_inbounds[0],x0x1x2_inbounds[1],x0x1x2_inbounds[2])\n",
    "  (see description of\n",
    "  EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt()\n",
    "  above for more details), here we compute the parity conditions\n",
    "  for all 10 tensor types supported by NRPy+.\n",
    "\"\"\"\n",
    "    c_type = \"static void\"\n",
    "    name = \"set_parity_for_inner_boundary_single_pt\"\n",
    "    params = \"\"\"const paramstruct *restrict params, const REAL xx0,const REAL xx1,const REAL xx2,\n",
    "                                             const REAL x0x1x2_inbounds[3], const int idx,\n",
    "                                             innerpt_bc_struct *restrict innerpt_bc_arr\"\"\"\n",
    "    body = r\"\"\"\n",
    "  const REAL xx0_inbounds = x0x1x2_inbounds[0];\n",
    "  const REAL xx1_inbounds = x0x1x2_inbounds[1];\n",
    "  const REAL xx2_inbounds = x0x1x2_inbounds[2];\n",
    "\n",
    "  REAL REAL_parity_array[10];\n",
    "  {\n",
    "    // Evaluate dot products needed for setting parity\n",
    "    //     conditions at a given point (xx0,xx1,xx2),\n",
    "    //     using C code generated by NRPy+\n",
    "\"\"\" + parity_conditions_symbolic_dot_products() + r\"\"\"\n",
    "  }\n",
    "\n",
    "  // Next perform sanity check on parity array output: should be +1 or -1 to within 8 significant digits:\n",
    "  for(int whichparity=0;whichparity<10;whichparity++) {\n",
    "    //printf(\"Good? Parity %d evaluated to %e\\n\",whichparity,(double)REAL_parity_array[whichparity]);\n",
    "    if( fabs(REAL_parity_array[whichparity]) < 1 - 1e-8 || fabs(REAL_parity_array[whichparity]) > 1 + 1e-8 ) {\n",
    "      fprintf(stderr,\"Error at point (%e %e %e), which maps to (%e %e %e).\\n\",\n",
    "              xx0,xx1,xx2, xx0_inbounds,xx1_inbounds,xx2_inbounds);\n",
    "      fprintf(stderr,\"Parity evaluated to %e , which is not within 8 significant digits of +1 or -1.\\n\",\n",
    "              REAL_parity_array[whichparity]);\n",
    "      exit(1);\n",
    "    }\n",
    "    // The typecast (int8_t)REAL_parity_array[whichparity] *does not work*.\n",
    "    //  Thankfully, we've already checked whether REAL_parity_array[whichparity]\n",
    "    //  is within 1e-8 of +/- 1, so here we just check the sign of the\n",
    "    //  REAL_parity_array to find the correct value of innerpt_bc_arr[idx].parity[parity].\n",
    "    for(int parity=0;parity<10;parity++) {\n",
    "      innerpt_bc_arr[idx].parity[parity] = 1;\n",
    "      if(REAL_parity_array[parity] < 0) innerpt_bc_arr[idx].parity[parity] = -1;\n",
    "    }\n",
    "  } // END for(int whichparity=0;whichparity<10;whichparity++)\n",
    "\"\"\"\n",
    "    __function_prototype_ignore, Cfunc = Cfunction(\n",
    "        desc=desc, c_type=c_type, name=name, params=params,\n",
    "        body=body)\n",
    "    return Cfunc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bcstruct_set_up_add_to_cfunc'></a>\n",
    "\n",
    "## Step 4.c: Add `bcstruct_set_up()` to NRPy+ C function dictionary \\[Back to [top](#toc)\\]\n",
    "$$\\label{bcstruct_set_up_add_to_cfunc}$$\n",
    "\n",
    "With the core routines `EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt` and `set_parity_for_inner_boundary_single_pt()` now defined, we can now add `bcstruct_set_up()` to the NRPy+ C function dictionary.\n",
    "\n",
    "At each inner and outer boundary point, this function sets the `innerpt_bc_struct` and `outerpt_bc_struct`, respectively. Note that\n",
    "\n",
    "1. Inner boundary points always map to a different grid point. An inner boundary point falls into one of the following two categories\n",
    "    1. pure inner boundary points (that map to a grid point entirely within the grid interior), or\n",
    "    1. inner-maps-to-outer boundary points (inner boundary points that map to a grid point on the outer boundary).\n",
    "1. Outer boundary points always map to themselves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:15.236052Z",
     "iopub.status.busy": "2021-09-08T15:47:15.235636Z",
     "iopub.status.idle": "2021-09-08T15:47:15.237789Z",
     "shell.execute_reply": "2021-09-08T15:47:15.237434Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_bcstruct_set_up(rel_path_to_Cparams=os.path.join(\".\")):\n",
    "    includes = [os.path.join(rel_path_to_Cparams, \"NRPy_basic_defines.h\"),\n",
    "                os.path.join(rel_path_to_Cparams, \"NRPy_function_prototypes.h\")]\n",
    "    prefunc  = Cfunction__EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt()\n",
    "    prefunc += Cfunction__set_parity_for_inner_boundary_single_pt()\n",
    "    desc = r\"\"\"At each coordinate point (x0,x1,x2) situated at grid index (i0,i1,i2):\n",
    "\n",
    "Step 1: Set up inner boundary structs bcstruct->inner_bc_array[].\n",
    "  Recall that at each inner boundary point we must set innerpt_bc_struct:\n",
    "    typedef struct __innerpt_bc_struct__ {\n",
    "      int dstpt;  // dstpt is the 3D grid index IDX3S(i0,i1,i2) of the inner boundary point (i0,i1,i2)\n",
    "      int srcpt;  // srcpt is the 3D grid index (a la IDX3S) to which the inner boundary point maps\n",
    "      int8_t parity[10];  // parity[10] is a calculation of dot products for the 10 independent parity types\n",
    "    } innerpt_bc_struct;\n",
    "  At each ghostzone (i.e., each point within NGHOSTS points from grid boundary):\n",
    "    Call EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt().\n",
    "        This function converts the curvilinear coordinate (x0,x1,x2) to the corresponding\n",
    "        Cartesian coordinate (x,y,z), then finds the grid point\n",
    "        (i0_inbounds,i1_inbounds,i2_inbounds) in the grid interior or outer boundary\n",
    "        corresponding to this Cartesian coordinate (x,y,z).\n",
    "    If (i0,i1,i2) *is not* the same as (i0_inbounds,i1_inbounds,i2_inbounds),\n",
    "        then we are at an inner boundary point. We must set\n",
    "        Set bcstruct->inner_bc_array for this point, which requires we specify\n",
    "        both (i0_inbounds,i1_inbounds,i2_inbounds) [just found!] and parity\n",
    "        conditions for this gridpoint. The latter is found & specified within the\n",
    "        function set_parity_for_inner_boundary_single_pt().\n",
    "    If (i0,i1,i2) *is* the same as (i0_inbounds,i1_inbounds,i2_inbounds),\n",
    "        then we are at an outer boundary point. Take care of outer BCs in Step 2.\n",
    "Step 2: Set up outer boundary structs bcstruct->outer_bc_array[which_gz][face][idx2d]:\n",
    "  Recall that at each inner boundary point we must set outerpt_bc_struct:\n",
    "    typedef struct __outerpt_bc_struct__ {\n",
    "      short i0,i1,i2;  // the outer boundary point grid index (i0,i1,i2), on the 3D grid\n",
    "      int8_t FACEX0,FACEX1,FACEX2;  // 1-byte integers that store\n",
    "      //                               FACEX0,FACEX1,FACEX2 = +1, 0, 0 if on the i0=i0min face,\n",
    "      //                               FACEX0,FACEX1,FACEX2 = -1, 0, 0 if on the i0=i0max face,\n",
    "      //                               FACEX0,FACEX1,FACEX2 =  0,+1, 0 if on the i1=i2min face,\n",
    "      //                               FACEX0,FACEX1,FACEX2 =  0,-1, 0 if on the i1=i1max face,\n",
    "      //                               FACEX0,FACEX1,FACEX2 =  0, 0,+1 if on the i2=i2min face, or\n",
    "      //                               FACEX0,FACEX1,FACEX2 =  0, 0,-1 if on the i2=i2max face,\n",
    "    } outerpt_bc_struct;\n",
    "  Outer boundary points are filled from the inside out, two faces at a time.\n",
    "    E.g., consider a Cartesian coordinate grid that has 14 points in each direction,\n",
    "    including the ghostzones, with NGHOSTS=2.\n",
    "    We first fill in the lower x0 face with (i0=1,i1={2,11},i2={2,11}). We fill these\n",
    "    points in first, since they will in general (at least in the case of extrapolation\n",
    "    outer BCs) depend on e.g., i0=2 and i0=3 points.\n",
    "    Simultaneously we can fill in the upper x0 face with (i0=12,i1={2,11},i2={2,11}),\n",
    "    since these points depend only on e.g., i0=11 and i0=10 (again assuming extrap. BCs).\n",
    "    Next we can fill in the lower x1 face: (i0={1,12},i1=2,i2={2,11}). Notice these\n",
    "    depend on i0 min and max faces being filled. The remaining pattern goes like this:\n",
    "    Upper x1 face: (i0={1,12},i1=12,i2={2,11})\n",
    "    Lower x2 face: (i0={1,12},i1={1,12},i2=1)\n",
    "    Upper x2 face: (i0={1,12},i1={1,12},i2=12)\n",
    "    Lower x0 face: (i0=0,i1={1,12},i2={1,12})\n",
    "    Upper x0 face: (i0=13,i1={1,12},i2={1,12})\n",
    "    Lower x1 face: (i0={0,13},i1=0,i2={2,11})\n",
    "    Upper x1 face: (i0={0,13},i1=13,i2={2,11})\n",
    "    Lower x2 face: (i0={0,13},i1={0,13},i2=0)\n",
    "    Upper x2 face: (i0={0,13},i1={0,13},i2=13)\n",
    "  Note that we allocate a outerpt_bc_struct at *all* boundary points,\n",
    "    regardless of whether the point is an outer or inner point. However\n",
    "    the struct is set only at outer boundary points. This is slightly\n",
    "    wasteful, but only in memory, not in CPU.\n",
    "\"\"\"\n",
    "    c_type = \"void\"\n",
    "    name = \"bcstruct_set_up\"\n",
    "    params = \"const paramstruct *restrict params, REAL *restrict xx[3], bc_struct *restrict bcstruct\"\n",
    "    body = r\"\"\"\n",
    "  ////////////////////////////////////////\n",
    "  // STEP 1: SET UP INNER BOUNDARY STRUCTS\n",
    "  {\n",
    "    // First count the number of inner points.\n",
    "    int num_inner = 0;\n",
    "    LOOP_OMP(\"omp parallel for reduction(+:num_inner)\",\n",
    "             i0,0,Nxx_plus_2NGHOSTS0,  i1,0,Nxx_plus_2NGHOSTS1,  i2,0,Nxx_plus_2NGHOSTS2) {\n",
    "      const int i0i1i2[3] = { i0,i1,i2 };\n",
    "      if(!IS_IN_GRID_INTERIOR(i0i1i2, Nxx_plus_2NGHOSTS0,Nxx_plus_2NGHOSTS1,Nxx_plus_2NGHOSTS2, NGHOSTS)) {\n",
    "        REAL x0x1x2_inbounds[3];\n",
    "        int i0i1i2_inbounds[3];\n",
    "        EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt(params, xx, i0,i1,i2, x0x1x2_inbounds,i0i1i2_inbounds);\n",
    "        if(i0 == i0i1i2_inbounds[0] && i1==i0i1i2_inbounds[1] && i2==i0i1i2_inbounds[2]) {\n",
    "          // this is a pure outer boundary point.\n",
    "        } else {\n",
    "          // this is an inner boundary point, which maps either\n",
    "          //  to the grid interior or to an outer boundary point\n",
    "          num_inner++;\n",
    "        }\n",
    "      }\n",
    "    }\n",
    "    // Store num_inner to bc_info:\n",
    "    bcstruct->bc_info.num_inner_boundary_points = num_inner;\n",
    "\n",
    "    // Next allocate memory for inner_boundary_points:\n",
    "    bcstruct->inner_bc_array = (innerpt_bc_struct *restrict)malloc( sizeof(innerpt_bc_struct)*num_inner );\n",
    "  }\n",
    "\n",
    "  // Then set inner_bc_array:\n",
    "  {\n",
    "    int which_inner = 0;\n",
    "    LOOP_NOOMP(i0,0,Nxx_plus_2NGHOSTS0,  i1,0,Nxx_plus_2NGHOSTS1,  i2,0,Nxx_plus_2NGHOSTS2) {\n",
    "      const int i0i1i2[3] = { i0,i1,i2 };\n",
    "      if(!IS_IN_GRID_INTERIOR(i0i1i2, Nxx_plus_2NGHOSTS0,Nxx_plus_2NGHOSTS1,Nxx_plus_2NGHOSTS2, NGHOSTS)) {\n",
    "        REAL x0x1x2_inbounds[3];\n",
    "        int i0i1i2_inbounds[3];\n",
    "        EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt(params, xx, i0,i1,i2, x0x1x2_inbounds,i0i1i2_inbounds);\n",
    "        if(i0 == i0i1i2_inbounds[0] && i1==i0i1i2_inbounds[1] && i2==i0i1i2_inbounds[2]) {\n",
    "          // this is a pure outer boundary point.\n",
    "        } else {\n",
    "          bcstruct->inner_bc_array[which_inner].dstpt = IDX3S(i0,i1,i2);\n",
    "          bcstruct->inner_bc_array[which_inner].srcpt = IDX3S(i0i1i2_inbounds[0],i0i1i2_inbounds[1],i0i1i2_inbounds[2]);\n",
    "          //printf(\"%d / %d\\n\",which_inner, bc_info->num_inner_boundary_points);\n",
    "          set_parity_for_inner_boundary_single_pt(params, xx[0][i0],xx[1][i1],xx[2][i2],\n",
    "                                                  x0x1x2_inbounds, which_inner, bcstruct->inner_bc_array);\n",
    "\n",
    "          which_inner++;\n",
    "        }\n",
    "      }\n",
    "    }\n",
    "  }\n",
    "\n",
    "  ////////////////////////////////////////\n",
    "  // STEP 2: SET UP OUTER BOUNDARY STRUCTS\n",
    "  // First set up loop bounds for outer boundary condition updates,\n",
    "  //   store to bc_info->bc_loop_bounds[which_gz][face][]. Also\n",
    "  //   allocate memory for outer_bc_array[which_gz][face][]:\n",
    "  int imin[3] = { NGHOSTS, NGHOSTS, NGHOSTS };\n",
    "  int imax[3] = { Nxx_plus_2NGHOSTS0-NGHOSTS, Nxx_plus_2NGHOSTS1-NGHOSTS, Nxx_plus_2NGHOSTS2-NGHOSTS };\n",
    "  for(int which_gz=0;which_gz<NGHOSTS;which_gz++) {\n",
    "    const int x0min_face_range[6] = { imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2] };  imin[0]--;\n",
    "    const int x0max_face_range[6] = { imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2] };  imax[0]++;\n",
    "    const int x1min_face_range[6] = { imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2] };  imin[1]--;\n",
    "    const int x1max_face_range[6] = { imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2] };  imax[1]++;\n",
    "    const int x2min_face_range[6] = { imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2] };  imin[2]--;\n",
    "    const int x2max_face_range[6] = { imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1 };  imax[2]++;\n",
    "\n",
    "    int face=0;\n",
    "    ////////////////////////\n",
    "    // x0min and x0max faces: Allocate memory for outer_bc_array and set bc_loop_bounds:\n",
    "    //                        Note that x0min and x0max faces have exactly the same size.\n",
    "    //                   Also, note that face/2 --v   offsets this factor of 2 ------------------------------------------v\n",
    "    bcstruct->pure_outer_bc_array[3*which_gz + face/2] = (outerpt_bc_struct *restrict)malloc(sizeof(outerpt_bc_struct) * 2 *\n",
    "                                                                                             ((x0min_face_range[1]-x0min_face_range[0]) *\n",
    "                                                                                              (x0min_face_range[3]-x0min_face_range[2]) *\n",
    "                                                                                              (x0min_face_range[5]-x0min_face_range[4])));\n",
    "    // x0min face: Can't set bc_info->bc_loop_bounds[which_gz][face] = { i0min,i0max, ... } since it's not const :(\n",
    "    for(int i=0;i<6;i++) { bcstruct->bc_info.bc_loop_bounds[which_gz][face][i] = x0min_face_range[i]; }\n",
    "    face++;\n",
    "    // x0max face: Set loop bounds & allocate memory for outer_bc_array:\n",
    "    for(int i=0;i<6;i++) { bcstruct->bc_info.bc_loop_bounds[which_gz][face][i] = x0max_face_range[i]; }\n",
    "    face++;\n",
    "    ////////////////////////\n",
    "\n",
    "    ////////////////////////\n",
    "    // x1min and x1max faces: Allocate memory for outer_bc_array and set bc_loop_bounds:\n",
    "    //                        Note that x1min and x1max faces have exactly the same size.\n",
    "    //                   Also, note that face/2 --v   offsets this factor of 2 ------------------------------------------v\n",
    "    bcstruct->pure_outer_bc_array[3*which_gz + face/2] = (outerpt_bc_struct *restrict)malloc(sizeof(outerpt_bc_struct) * 2 *\n",
    "                                                                                             ((x1min_face_range[1]-x1min_face_range[0]) *\n",
    "                                                                                              (x1min_face_range[3]-x1min_face_range[2]) *\n",
    "                                                                                              (x1min_face_range[5]-x1min_face_range[4])));\n",
    "    // x1min face: Can't set bc_info->bc_loop_bounds[which_gz][face] = { i0min,i0max, ... } since it's not const :(\n",
    "    for(int i=0;i<6;i++) { bcstruct->bc_info.bc_loop_bounds[which_gz][face][i] = x1min_face_range[i]; }\n",
    "    face++;\n",
    "    // x1max face: Set loop bounds & allocate memory for outer_bc_array:\n",
    "    for(int i=0;i<6;i++) { bcstruct->bc_info.bc_loop_bounds[which_gz][face][i] = x1max_face_range[i]; }\n",
    "    face++;\n",
    "    ////////////////////////\n",
    "\n",
    "\n",
    "    ////////////////////////\n",
    "    // x2min and x2max faces: Allocate memory for outer_bc_array and set bc_loop_bounds:\n",
    "    //                        Note that x2min and x2max faces have exactly the same size.\n",
    "    //                   Also, note that face/2 --v   offsets this factor of 2 ------------------------------------------v\n",
    "    bcstruct->pure_outer_bc_array[3*which_gz + face/2] = (outerpt_bc_struct *restrict)malloc(sizeof(outerpt_bc_struct) * 2 *\n",
    "                                                                                             ((x2min_face_range[1]-x2min_face_range[0]) *\n",
    "                                                                                              (x2min_face_range[3]-x2min_face_range[2]) *\n",
    "                                                                                              (x2min_face_range[5]-x2min_face_range[4])));\n",
    "    // x2min face: Can't set bc_info->bc_loop_bounds[which_gz][face] = { i0min,i0max, ... } since it's not const :(\n",
    "    for(int i=0;i<6;i++) { bcstruct->bc_info.bc_loop_bounds[which_gz][face][i] = x2min_face_range[i]; }\n",
    "    face++;\n",
    "    // x2max face: Set loop bounds & allocate memory for outer_bc_array:\n",
    "    for(int i=0;i<6;i++) { bcstruct->bc_info.bc_loop_bounds[which_gz][face][i] = x2max_face_range[i]; }\n",
    "    face++;\n",
    "    ////////////////////////\n",
    "  }\n",
    "\n",
    "  for(int which_gz=0;which_gz<NGHOSTS;which_gz++) for(int dirn=0;dirn<3;dirn++) {\n",
    "      int idx2d = 0;\n",
    "      // LOWER FACE: dirn=0 -> x0min; dirn=1 -> x1min; dirn=2 -> x2min\n",
    "      {\n",
    "        const int face = dirn*2;\n",
    "#define IDX2D_BCS(i0,i0min,i0max, i1,i1min,i1max ,i2,i2min,i2max)       \\\n",
    "        ( ((i0)-(i0min)) + ((i0max)-(i0min)) * ( ((i1)-(i1min)) + ((i1max)-(i1min)) * ((i2)-(i2min)) ) )\n",
    "        const int FACEX0=(face==0) - (face==1); // +1 if face==0 (x0min) ; -1 if face==1 (x0max). Otherwise 0.\n",
    "        const int FACEX1=(face==2) - (face==3); // +1 if face==2 (x1min) ; -1 if face==3 (x1max). Otherwise 0.\n",
    "        const int FACEX2=(face==4) - (face==5); // +1 if face==4 (x2min) ; -1 if face==5 (x2max). Otherwise 0.\n",
    "        LOOP_NOOMP(i0,bcstruct->bc_info.bc_loop_bounds[which_gz][face][0],bcstruct->bc_info.bc_loop_bounds[which_gz][face][1],\n",
    "                   i1,bcstruct->bc_info.bc_loop_bounds[which_gz][face][2],bcstruct->bc_info.bc_loop_bounds[which_gz][face][3],\n",
    "                   i2,bcstruct->bc_info.bc_loop_bounds[which_gz][face][4],bcstruct->bc_info.bc_loop_bounds[which_gz][face][5]) {\n",
    "          REAL x0x1x2_inbounds[3];\n",
    "          int i0i1i2_inbounds[3];\n",
    "          EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt(params, xx, i0,i1,i2, x0x1x2_inbounds,i0i1i2_inbounds);\n",
    "          if(i0 == i0i1i2_inbounds[0] && i1==i0i1i2_inbounds[1] && i2==i0i1i2_inbounds[2]) {\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i0 = i0;\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i1 = i1;\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i2 = i2;\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX0 = FACEX0;\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX1 = FACEX1;\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX2 = FACEX2;\n",
    "            idx2d++;\n",
    "          }\n",
    "        }\n",
    "      }\n",
    "      // UPPER FACE: dirn=0 -> x0max; dirn=1 -> x1max; dirn=2 -> x2max\n",
    "      {\n",
    "        const int face = dirn*2+1;\n",
    "        const int FACEX0=(face==0) - (face==1); // +1 if face==0 ; -1 if face==1. Otherwise 0.\n",
    "        const int FACEX1=(face==2) - (face==3); // +1 if face==2 ; -1 if face==3. Otherwise 0.\n",
    "        const int FACEX2=(face==4) - (face==5); // +1 if face==4 ; -1 if face==5. Otherwise 0.\n",
    "        LOOP_NOOMP(i0,bcstruct->bc_info.bc_loop_bounds[which_gz][face][0],bcstruct->bc_info.bc_loop_bounds[which_gz][face][1],\n",
    "                   i1,bcstruct->bc_info.bc_loop_bounds[which_gz][face][2],bcstruct->bc_info.bc_loop_bounds[which_gz][face][3],\n",
    "                   i2,bcstruct->bc_info.bc_loop_bounds[which_gz][face][4],bcstruct->bc_info.bc_loop_bounds[which_gz][face][5]) {\n",
    "          REAL x0x1x2_inbounds[3];\n",
    "          int i0i1i2_inbounds[3];\n",
    "          EigenCoord_set_x0x1x2_inbounds__i0i1i2_inbounds_single_pt(params, xx, i0,i1,i2, x0x1x2_inbounds,i0i1i2_inbounds);\n",
    "          if(i0 == i0i1i2_inbounds[0] && i1==i0i1i2_inbounds[1] && i2==i0i1i2_inbounds[2]) {\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i0 = i0;\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i1 = i1;\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i2 = i2;\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX0 = FACEX0;\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX1 = FACEX1;\n",
    "            bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX2 = FACEX2;\n",
    "            idx2d++;\n",
    "          }\n",
    "        }\n",
    "      }\n",
    "      bcstruct->bc_info.num_pure_outer_boundary_points[which_gz][dirn] = idx2d;\n",
    "    }\n",
    "\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        prefunc=prefunc,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='apply_bcs_inner_only'></a>\n",
    "\n",
    "# Step 5: `apply_bcs_inner_only()`: Inner boundary condition C function \\[Back to [top](#toc)\\]\n",
    "$$\\label{apply_bcs_inner_only}$$\n",
    "\n",
    "After `bcstruct_set_up()`, `bcstruct` contains all information needed to fill inner boundary points. This includes\n",
    "\n",
    "1. The interior or outer boundary point to which the inner point maps.\n",
    "1. The tensor parity multiplier (+1 or -1) that the value of the gridfunction will need to be multiplied.\n",
    "\n",
    "This function will generally need to be called after outer boundary conditions have been applied, as inner boundary points can depend on outer boundary points having been filled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_apply_bcs_inner_only(rel_path_to_Cparams=os.path.join(\".\")):\n",
    "    includes = [os.path.join(rel_path_to_Cparams, \"NRPy_basic_defines.h\")]\n",
    "    desc = r\"\"\"\n",
    "Apply BCs to inner boundary points only,\n",
    "using data stored in bcstruct->inner_bc_array.\n",
    "These structs are set in bcstruct_set_up().\n",
    "Inner boundary points map to either the grid\n",
    "interior (\"pure inner\") or to pure outer\n",
    "boundary points (\"inner maps to outer\").\n",
    "\"\"\"\n",
    "    c_type = \"void\"\n",
    "    name = \"apply_bcs_inner_only\"\n",
    "    params = \"const paramstruct *restrict params, const bc_struct *restrict bcstruct, REAL *restrict gfs\"\n",
    "    body = r\"\"\"\n",
    "  // Unpack bc_info from bcstruct\n",
    "  const bc_info_struct *bc_info = &bcstruct->bc_info;\n",
    "\n",
    "  // collapse(2) results in a nice speedup here, esp in 2D. Two_BHs_collide goes from\n",
    "  //    5550 M/hr to 7264 M/hr on a Ryzen 9 5950X running on all 16 cores with core affinity.\n",
    "#pragma omp parallel for collapse(2)  // spawn threads and distribute across them\n",
    "  for(int which_gf=0;which_gf<NUM_EVOL_GFS;which_gf++) {\n",
    "    for(int pt=0;pt<bc_info->num_inner_boundary_points;pt++) {\n",
    "      const int dstpt = bcstruct->inner_bc_array[pt].dstpt;\n",
    "      const int srcpt = bcstruct->inner_bc_array[pt].srcpt;\n",
    "      gfs[IDX4ptS(which_gf, dstpt)] = bcstruct->inner_bc_array[pt].parity[evol_gf_parity[which_gf]] * gfs[IDX4ptS(which_gf, srcpt)];\n",
    "    } // END for(int pt=0;pt<num_inner_pts;pt++)\n",
    "  } // END for(int which_gf=0;which_gf<NUM_EVOL_GFS;which_gf++)\n",
    "\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='outer_bcs'></a>\n",
    "\n",
    "# Step 6: Outer boundary condition C functions \\[Back to [top](#toc)\\]\n",
    "$$\\label{outer_bcs}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='extrapolation_bcs'></a>\n",
    "\n",
    "## Step 6.a: `\"extrapolation\"` outer boundary conditions: apply quadratic polynomial extrapolation \\[Back to [top](#toc)\\]\n",
    "$$\\label{extrapolation_bcs}$$\n",
    "\n",
    "As an option, quadratic extrapolation may be applied to each outer boundary point `(i0,i1,i2)`, as follows.\n",
    "\n",
    "Suppose the outer boundary point is at the `i0=max(i0)` face. Then we fit *known* data at `i0-3`, `i0-2`, and `i0-1` [i.e., $f_{-3}=f(x_{i0-3}=x_{-3})$, $f_{-2}=f(x_{i0-2}=x_{-2})$, and $f_{-1}=f(x_{i0-1}=x_{-1})$] to the unique quadratic polynomial:\n",
    "\n",
    "\\begin{align}\n",
    "f_{-3} &= c_2 x_{-3}^2 + c_1 x_{-3} + c_0 \\\\\n",
    "f_{-2} &= c_2 x_{-2}^2 + c_1 x_{-2} + c_0 \\\\\n",
    "f_{-1} &= c_2 x_{-1}^2 + c_1 x_{-1} + c_0 \\\\\n",
    "\\end{align}\n",
    "\n",
    "We wish to extrapolate to $f_0=f(x_0)$. Since our grid has uniform spacing, \n",
    "\n",
    "* $x_{-3}=x_0-3\\Delta x$, \n",
    "* $x_{-2}=x_0-2\\Delta x$, and\n",
    "* $x_{-1}=x_0-\\Delta x$.\n",
    "\n",
    "The extrapolated value $f_0$ cannot depend on the choice of the fiducial $x_0$ (i.e., it will hold for *any* choice of $x_0$), so without loss of generality we will set $x_0=0$:\n",
    "\n",
    "$$\n",
    "\\mathbf{A c} =\n",
    "\\left[\n",
    "\\begin{array}{ccc}\n",
    " 1 & x_{-3}  & x_{-3}^2 \\\\\n",
    " 1 & x_{-2}  & x_{-2}^2 \\\\\n",
    " 1 & x_{-1}  & x_{-1}^2 \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    " c_0 \\\\\n",
    " c_1 \\\\\n",
    " c_2 \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{array}{ccc}\n",
    " 1 & -3 \\Delta x  & 9 \\Delta x^2 \\\\\n",
    " 1 & -2 \\Delta x  & 4 \\Delta x^2 \\\\\n",
    " 1 & -1 \\Delta x  &   \\Delta x^2 \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    " c_0 \\\\\n",
    " c_1 \\\\\n",
    " c_2 \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{array}{ccc}\n",
    " 1 & -3  & 9 \\\\\n",
    " 1 & -2  & 4 \\\\\n",
    " 1 & -1  & 1  \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    " c_0 \\\\\n",
    " c_1 \\Delta x \\\\\n",
    " c_2 \\Delta x^2 \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    " f_{-3} \\\\\n",
    " f_{-2} \\\\\n",
    " f_{-1} \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "= \\mathbf{f}.\n",
    "$$\n",
    "\n",
    "This is known as the [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix) for the quadratic polynomial, and its solution for $c_0$, $c_1$, and $c_2$ will be unique. But before we invert the matrix, we note that as we wish to solve for $f(x_0)=f(0)$ (again, the extrapolated value $f(x_0)$ cannot depend on the choice of the fiducial $x_0$; i.e., it will hold for *any* choice of $x_0$ so without loss of generality we will set $x_0=0$), our quadratic polynomial simplifies to:\n",
    "\n",
    "$$\n",
    "f_{0} = c_2 x_{0}^2 + c_1 x_{0} + c_0 = c_0.\n",
    "$$\n",
    "\n",
    "Thus we need only extract the value of $c_0$, which is done in the next cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "⎡ 1   -3   3 ⎤\n",
      "⎢            ⎥\n",
      "⎢3/2  -4  5/2⎥\n",
      "⎢            ⎥\n",
      "⎣1/2  -1  1/2⎦\n"
     ]
    }
   ],
   "source": [
    "from sympy import symbols,Matrix,factor,pretty_print,simplify,latex\n",
    "MM = Matrix([[1,-3,9],\n",
    "             [1,-2,4],\n",
    "             [1,-1,1]])\n",
    "# print(latex(factor(MM.inv())))\n",
    "pretty_print(factor(MM.inv()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus we get \n",
    "\n",
    "$$\n",
    "f_0 = c_0 = f_{-3} - 3 f_{-2} + 3 f_{-1}.\n",
    "$$\n",
    "To determine the coefficient at the `i0=min(i0)` face, the above analysis can be repeated:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "⎡ 1    -3   3  ⎤\n",
      "⎢              ⎥\n",
      "⎢-3/2  4   -5/2⎥\n",
      "⎢              ⎥\n",
      "⎣1/2   -1  1/2 ⎦\n"
     ]
    }
   ],
   "source": [
    "from sympy import symbols,Matrix,factor,pretty_print,simplify,latex\n",
    "MM = Matrix([[1,+3,9],\n",
    "             [1,+2,4],\n",
    "             [1,+1,1]])\n",
    "# print(latex(factor(MM.inv())))\n",
    "pretty_print(factor(MM.inv()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and we find the $c_0$ coefficient is basically the same as on the `i0=max(i0)` face; just replace $f_{-3}\\to f_{+3}$, etc:\n",
    "\n",
    "$$\n",
    "f_0 = c_0 = f_{+3} - 3 f_{+2} + 3 f_{+1}.\n",
    "$$\n",
    "\n",
    "The resulting extrapolation algorithm appears in the `apply_bcs_outerextrap_and_inner()` function below. I.e.,\n",
    "\n",
    "```c\n",
    "// *** Apply 2nd-order polynomial extrapolation BCs to all outer boundary points. ***\n",
    "gfs[IDX4ptS(which_gf, idx_offset0)] =\n",
    "  +3.0*gfs[IDX4ptS(which_gf, idx_offset1)]\n",
    "  -3.0*gfs[IDX4ptS(which_gf, idx_offset2)]\n",
    "  +1.0*gfs[IDX4ptS(which_gf, idx_offset3)];\n",
    "```\n",
    "\n",
    "This function is meant to be called within a [Method of Lines timestepping algorithm](Tutorial-Method_of_Lines-C_Code_Generation.ipynb), at the very end of each MoL substep. As its name implies, it also updates the inner boundary points, which must be updated after the outer boundary points as inner points can map to outer boundary points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_apply_bcs_outerextrap_and_inner(rel_path_to_Cparams=os.path.join(\".\")):\n",
    "    includes = [os.path.join(rel_path_to_Cparams, \"NRPy_basic_defines.h\"),\n",
    "                os.path.join(rel_path_to_Cparams, \"NRPy_function_prototypes.h\")]\n",
    "    desc = r\"\"\"\n",
    "Suppose the outer boundary point is at the i0=max(i0) face. Then we fit known data at i0-3, i0-2, and i0-1\n",
    "  to the unique quadratic polynomial that passes through those points, and fill the data at\n",
    "  i0 with the value implied from the polynomial.\n",
    "As derived in nrpytutorial's Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb,\n",
    "  the coefficients must be f_{i0} = f_{i0-3} - 3 f_{i0-2} + 3 f_{i0-1}.\n",
    "  To check these coefficients are correct, consider\n",
    "  * f(x0 = constant. Then f_{i0} = f_{i0-3} <- CHECK!\n",
    "  * f(x) = x. WOLOG suppose x0=0. Then f_{i0} = (-3dx) - 3(-2dx) + 3(-dx) = + dx(-3+6-3) = 0 <- CHECK!\n",
    "  * f(x) = x^2. WOLOG suppose x0=0. Then f_{i0} = (-3dx)^2 - 3(-2dx)^2 + 3(-dx)^2 = + dx^2(9-12+3) = 0 <- CHECK!\n",
    "\"\"\"\n",
    "    c_type = \"void\"\n",
    "    name = \"apply_bcs_outerextrap_and_inner\"\n",
    "    params = \"const paramstruct *restrict params, const bc_struct *restrict bcstruct, REAL *restrict gfs\"\n",
    "    body = r\"\"\"\n",
    "  // Unpack bc_info from bcstruct\n",
    "  const bc_info_struct *bc_info = &bcstruct->bc_info;\n",
    "\n",
    "  ////////////////////////////////////////////////////////\n",
    "  // STEP 1 of 2: Apply BCs to pure outer boundary points.\n",
    "  //              By \"pure\" we mean that these points are\n",
    "  //              on the outer boundary and not also on\n",
    "  //              an inner boundary.\n",
    "  //              Here we fill in the innermost ghost zone\n",
    "  //              layer first and move outward. At each\n",
    "  //              layer, we fill in +/- x0 faces first,\n",
    "  //              then +/- x1 faces, finally +/- x2 faces,\n",
    "  //              filling in the edges as we go.\n",
    "  // Spawn N OpenMP threads, either across all cores, or according to e.g., taskset.\n",
    "#pragma omp parallel\n",
    "  {\n",
    "    for(int which_gz=0;which_gz<NGHOSTS;which_gz++) for(int dirn=0;dirn<3;dirn++) {\n",
    "        // This option results in about 1.6% slower runtime for SW curvilinear at 64x24x24 on 8-core Ryzen 9 4900HS\n",
    "        //#pragma omp for collapse(2)\n",
    "        //for(int which_gf=0;which_gf<NUM_EVOL_GFS;which_gf++) for(int idx2d=0;idx2d<bc_info->num_pure_outer_boundary_points[which_gz][dirn];idx2d++) {\n",
    "        //  {\n",
    "        // Don't spawn a thread if there are no boundary points to fill; results in a nice little speedup.\n",
    "        if(bc_info->num_pure_outer_boundary_points[which_gz][dirn] > 0) {\n",
    "#pragma omp for  // threads have been spawned; here we distribute across them\n",
    "          for(int idx2d=0;idx2d<bc_info->num_pure_outer_boundary_points[which_gz][dirn];idx2d++) {\n",
    "            const short i0 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i0;\n",
    "            const short i1 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i1;\n",
    "            const short i2 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i2;\n",
    "            const short FACEX0 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX0;\n",
    "            const short FACEX1 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX1;\n",
    "            const short FACEX2 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX2;\n",
    "            const int idx_offset0 = IDX3S(i0,i1,i2);\n",
    "            const int idx_offset1 = IDX3S(i0+1*FACEX0,i1+1*FACEX1,i2+1*FACEX2);\n",
    "            const int idx_offset2 = IDX3S(i0+2*FACEX0,i1+2*FACEX1,i2+2*FACEX2);\n",
    "            const int idx_offset3 = IDX3S(i0+3*FACEX0,i1+3*FACEX1,i2+3*FACEX2);\n",
    "            for(int which_gf=0;which_gf<NUM_EVOL_GFS;which_gf++) {\n",
    "              // *** Apply 2nd-order polynomial extrapolation BCs to all outer boundary points. ***\n",
    "              gfs[IDX4ptS(which_gf, idx_offset0)] =\n",
    "                +3.0*gfs[IDX4ptS(which_gf, idx_offset1)]\n",
    "                -3.0*gfs[IDX4ptS(which_gf, idx_offset2)]\n",
    "                +1.0*gfs[IDX4ptS(which_gf, idx_offset3)];\n",
    "            }\n",
    "          }\n",
    "        }\n",
    "      }\n",
    "  }\n",
    "\n",
    "  ///////////////////////////////////////////////////////\n",
    "  // STEP 2 of 2: Apply BCs to inner boundary points.\n",
    "  //              These map to either the grid interior\n",
    "  //              (\"pure inner\") or to pure outer boundary\n",
    "  //              points (\"inner maps to outer\"). Those\n",
    "  //              that map to outer require that outer be\n",
    "  //              populated first; hence this being\n",
    "  //              STEP 2 OF 2.\n",
    "  apply_bcs_inner_only(params, bcstruct, gfs);\n",
    "\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='radiation_bcs'></a>\n",
    "\n",
    "## Step 6.b: `\"radiation\"` outer boundary conditions \\[Back to [top](#toc)\\]\n",
    "$$\\label{radiation_bcs}$$\n",
    "\n",
    "<a id='radiation_bcs_theory'></a>\n",
    "\n",
    "### Step 6.b.i: Core *ansatz*, and implications \\[Back to [top](#toc)\\]\n",
    "$$\\label{radiation_bcs_theory}$$\n",
    "\n",
    "Here we describe outgoing radiation (a.k.a., Sommerfeld) boundary conditions.\n",
    "\n",
    "Our implementation follows that of the [Einstein Toolkit](https://einsteintoolkit.org/)'s [`NewRad` thorn](https://einsteintoolkit.org/thornguide/EinsteinEvolve/NewRad/documentation.html), but extends to arbitrary finite-difference order.\n",
    "\n",
    "The basic *ansatz* is that at **each outer boundary point**, an arbitrary field $f$ behaves locally as an outgoing spherical wave:\n",
    "\n",
    "$$\n",
    "f(r,t) = f_{r\\to\\infty} + \\frac{w(r - ct)}{r} + \\frac{K}{r^2},\n",
    "$$\n",
    "where\n",
    "\n",
    "* $f_{r\\to\\infty}$ is the (constant) value the field $f$ tends to as $r\\to\\infty$.\n",
    "* $w(r - ct)/r$ is a solution to the spherically symmetric wave equation ($\\partial_t^2 (r w) - c^2 \\partial_r^2 (r w) = 0$) for an *outgoing* spherical wave.\n",
    "* $K/r^2$ corrects for the next-leading-order radial falloff, where $K$ is a constant determined at each boundary point $(r,\\theta,\\phi)$ at each time $t$ by analyzing the behavior of $f(r)$ just inside the outer boundary of the numerical domain. The full numerical implementation of this term is described below.\n",
    "\n",
    "This boundary condition approach thus contains two free parameters: $c$ (the wavespeed for field $f$) and $f_{r\\to\\infty}$. Both of these should be known for any given field prior to attempting a solution.\n",
    "\n",
    "For convenience, we apply this boundary condition not to $f$ itself but to $\\partial_t f$, which will always be computed in a [Method of Lines](Tutorial-Method_of_Lines-C_Code_Generation.ipynb) approach within `NRPy+`.\n",
    "\n",
    "Taking the first time derivative we get:\n",
    "\n",
    "$$\n",
    "\\partial_t f = -c \\frac{w'(r - ct)}{r}.\n",
    "$$\n",
    "\n",
    "As $w$ represents an outgoing wave, in which temporal and spatial derivatives are directly related, we will find it convenient to also compute the radial derivative $\\partial_r f$ as well:\n",
    "\n",
    "$$\n",
    "\\partial_r f = \\frac{w'(r - ct)}{r} - \\frac{w(r - ct)}{r^2} - 2\\frac{K}{r^3}.\n",
    "$$\n",
    "\n",
    "Combining these two equations we get:\n",
    "\n",
    "$$\n",
    "\\partial_t f = -c \\left(\\partial_r f + \\frac{w(r - ct)}{r^2} + 2\\frac{K}{r^3}\\right).\n",
    "$$\n",
    "\n",
    "Next we use the *ansatz* to compute $w(r - ct)/r^2$:\n",
    "\n",
    "\\begin{align}\n",
    "f(r,t) &= f_{r\\to\\infty} + \\frac{w(r - ct)}{r} + \\frac{K}{r^2} \\\\\n",
    "\\implies \\frac{w(r - ct)}{r^2} &= \\frac{f - f_{r\\to\\infty}}{r} - \\frac{K}{r^3}.\n",
    "\\end{align}\n",
    "\n",
    "This enables us to rewrite $\\partial_t f$ as\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_t f &= -c \\left(\\partial_r f + \\frac{f - f_{r\\to\\infty}}{r} - \\frac{K}{r^3} + 2\\frac{K}{r^3}\\right) \\\\\n",
    "&= -c \\left(\\partial_r f + \\frac{f - f_{r\\to\\infty}}{r} + \\frac{K}{r^3}\\right) \\\\\n",
    "&= -\\frac{c}{r} \\left[r \\partial_r f + \\left(f - f_{r\\to\\infty}\\right)\\right] + \\frac{k}{r^3},\n",
    "\\end{align}\n",
    "where $k=-Kc$ just re-expresses the unknown function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='radiation_bcs_numerical'></a>\n",
    "\n",
    "### Step 6.b.ii: Numerical implementation \\[Back to [top](#toc)\\]\n",
    "$$\\label{radiation_bcs_numerical}$$\n",
    "\n",
    "\n",
    "We start with the equation derived in the previous section:\n",
    "$$\n",
    "\\partial_t f = -\\frac{c}{r} \\left[r \\partial_r f + \\left(f - f_{r\\to\\infty}\\right)\\right] + \\frac{k}{r^3}.\n",
    "$$\n",
    "\n",
    "First note that the right-hand side of this equation is applied *at the current time* (i.e., the same time at which $\\partial_t f$ is evaluated in the Method of Lines timestepping). I.e., it is applied prior to the MoL update. Further, at the current time, $f$ is known at *all* gridpoints (interior and boundary).\n",
    "\n",
    "As this boundary condition must be applicable to *any* curvilinear coordinate system $x^i$ (and not just spherical coordinates), $\\partial_r f$ must be computed via\n",
    "\n",
    "$$\n",
    "\\partial_r f = \\frac{\\partial x^i}{\\partial r} \\partial_{i} f = \\frac{\\partial x^i}{\\partial r} \\frac{\\partial f}{\\partial x^{i}},\n",
    "$$\n",
    "\n",
    "where $\\partial_i f$ is evaluated using finite-difference derivatives, suitably upwinded to avoid extending beyond the domain of the grid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='radiation_bcs_numerical_inv_jacobian'></a>\n",
    "\n",
    "### Step 6.b.iii: The $\\partial_r f$ term: Computing $\\frac{\\partial x^i}{\\partial r}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{radiation_bcs_numerical_inv_jacobian}$$\n",
    "\n",
    "The term $\\partial x^i/\\partial r$ cannot be immediately computed, as $x^i$ is never written as a function of $r$. However, in all curvilinear coordinates supported by NRPy+, $r$, $\\theta$, and $\\phi$ must be specified as functions of $x^i$. Thus we have exact expressions for the Jacobian:\n",
    "\n",
    "$$\n",
    "J_i^j = \\frac{\\partial x^j_{\\rm Sph}}{\\partial x^i_{\\rm Curv}},\n",
    "$$\n",
    "\n",
    "and $\\partial x^i/\\partial r$ can be computed from the inverse of this matrix (computed within `NRPy+`):\n",
    "\n",
    "$$\n",
    "(J^{-1})_j^i = \\frac{\\partial x^j_{\\rm Curv}}{\\partial x^i_{\\rm Sph}},\n",
    "$$\n",
    "and only looking at the $i=0$ component (as $x^i_{\\rm Sph}=r$).\n",
    "\n",
    "For example let's set the coordinate system to Cylindrical, as exact expressions do exist for cylindrical coordinates in terms of spherical coordinates for this case:\n",
    "\n",
    "* $\\rho = r \\sin\\theta$\n",
    "* $\\phi = \\phi$\n",
    "* $z = r \\cos\\theta$\n",
    "\n",
    "Thus\n",
    "\n",
    "$$\n",
    "\\partial_r \\rho = \\sin\\theta = \\frac{\\rho}{\\sqrt{\\rho^2 + z^2}},\n",
    "$$\n",
    "\n",
    "where the second equality can be computed directly from the inverse Jacobian. Note also that the negative root is not considered, as $\\theta \\in [0,\\pi]$.\n",
    "\n",
    "Again expressions for cylindrical coordinates in terms of spherical coordinates *do not exist* within `NRPy+` (as they are not generally easy to describe).\n",
    "\n",
    "Next, let's check that NRPy+ yields the correct result for this case, choosing Cylindrical coordinates, where $x^i=$ `(xx0,xx1,xx2)` =$(\\rho,\\phi,z)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "partial rho / partial r  =  xx0/sqrt(xx0**2 + xx2**2)\n",
      "partial phi / partial r  =  0\n",
      "partial z   / partial r  =  xx2/sqrt(xx0**2 + xx2**2)\n"
     ]
    }
   ],
   "source": [
    "par.set_parval_from_str(\"reference_metric::CoordSystem\", \"Cylindrical\")\n",
    "rfm.reference_metric()\n",
    "\n",
    "# Jac_dUrfm_dDSphUD is the inverse Jacobian\n",
    "Jac_dUSph_dDrfmUD, Jac_dUrfm_dDSphUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Spherical()\n",
    "\n",
    "print(\"partial rho / partial r  =  \" + str(sp.simplify(Jac_dUrfm_dDSphUD[0][0])))\n",
    "print(\"partial phi / partial r  =  \" + str(sp.simplify(Jac_dUrfm_dDSphUD[1][0])))\n",
    "print(\"partial z   / partial r  =  \" + str(sp.simplify(Jac_dUrfm_dDSphUD[2][0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define a NRPy+-generated inlined C function, `r_and_partial_xi_partial_r_derivs()` for computing $\\frac{\\partial x^i}{\\partial r}$. It turns out to be convenient for this function to also compute $r$, since both are needed:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_Cfunction_r_and_partial_xi_partial_r_derivs():\n",
    "    desc = \"Compute r(xx0,xx1,xx2) and partial_r x^i.\"\n",
    "    c_type = \"static inline void\"\n",
    "    name = \"r_and_partial_xi_partial_r_derivs\"\n",
    "    params = \"\"\"const paramstruct *restrict params,const REAL xx0,const REAL xx1,const REAL xx2,\n",
    "                                                     REAL *r,\n",
    "                                                     REAL *partial_x0_partial_r,REAL *partial_x1_partial_r,REAL *partial_x2_partial_r\"\"\"\n",
    "    Jac_dUSph_dDrfmUD, Jac_dUrfm_dDSphUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Spherical()\n",
    "    body = outputC([rfm.xxSph[0],\n",
    "                    Jac_dUrfm_dDSphUD[0][0], # sp.simplify(expr) is too slow here for SinhCylindrical\n",
    "                    Jac_dUrfm_dDSphUD[1][0], # sp.simplify(expr) is too slow here for SinhCylindrical\n",
    "                    Jac_dUrfm_dDSphUD[2][0]],# sp.simplify(expr) is too slow here for SinhCylindrical\n",
    "                   [\"*r\", \"*partial_x0_partial_r\", \"*partial_x1_partial_r\", \"*partial_x2_partial_r\"], filename=\"returnstring\",\n",
    "                   params=\"preindent=1,outCverbose=False,includebraces=False\")\n",
    "    rel_path_to_Cparams = os.path.join(\".\")\n",
    "\n",
    "    _prototype, func = Cfunction(\n",
    "        includes=[],   desc=desc, c_type=c_type, name=name, params=params,    body=body,\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)\n",
    "    return func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='radiation_bcs_numerical_partial_i_f'></a>\n",
    "\n",
    "### Step 6.b.iv: The $\\partial_r f$ term: Computing $\\partial_i f$ with arbitrary-offset finite-difference derivatives \\[Back to [top](#toc)\\]\n",
    "$$\\label{radiation_bcs_numerical_partial_i_f}$$\n",
    "\n",
    "Coefficients computed here can be and have been validated against the tables in this [Wikipedia article](https://en.wikipedia.org/wiki/Finite_difference_coefficient)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "At FD order=4 & stencil [-4, -3, -2, -1, 0], we get coeffs = [1/4, -4/3, 3, -4, 25/12]\n",
      "At FD order=4 & stencil [-3, -2, -1, 0, 1], we get coeffs = [-1/12, 1/2, -3/2, 5/6, 1/4]\n",
      "At FD order=4 & stencil [-2, -1, 0, 1, 2], we get coeffs = [1/12, -2/3, 0, 2/3, -1/12]\n",
      "At FD order=4 & stencil [-1, 0, 1, 2, 3], we get coeffs = [-1/4, -5/6, 3/2, -1/2, 1/12]\n",
      "At FD order=4 & stencil [0, 1, 2, 3, 4], we get coeffs = [-25/12, 4, -3, 4/3, -1/4]\n"
     ]
    }
   ],
   "source": [
    "import finite_difference as fin\n",
    "# import sympy as sp\n",
    "\n",
    "def get_arb_offset_FD_coeffs_indices(FDORDER, offset, deriv):\n",
    "    # deriv = 1 <-- 1st derivative\n",
    "    Minv = fin.setup_FD_matrix__return_inverse(FDORDER+1, offset)\n",
    "    indices = []\n",
    "    coeffs = []\n",
    "    for i in range(FDORDER+1):\n",
    "        indices.append(i-int(FDORDER/2) + offset)\n",
    "        coeffs.append(Minv[i, deriv])\n",
    "    return coeffs, indices\n",
    "\n",
    "FDORDER=4\n",
    "for offset in range(-int(FDORDER/2), int(FDORDER/2)+1):\n",
    "    coeffs, indices = get_arb_offset_FD_coeffs_indices(FDORDER, offset, 1)\n",
    "    print(\"At FD order=\" + str(FDORDER) + \" & stencil \" + str(indices) + \", we get coeffs = \" + str(coeffs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define a NRPy+-generated inlined C function, `FD1_arbitrary_upwind_xN_dirn()` for computing $\\partial_i f$ with an arbitrary upwind/downwind:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_Cfunction_FD1_arbitrary_upwind(dirn, radiation_BC_FD_order=-1):\n",
    "    default_FDORDER = par.parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\")\n",
    "    if radiation_BC_FD_order == -1:\n",
    "        radiation_BC_FD_order = default_FDORDER\n",
    "\n",
    "    par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", radiation_BC_FD_order)\n",
    "\n",
    "    includes = []\n",
    "    desc = \"Compute 1st derivative finite-difference derivative with arbitrary upwind\"\n",
    "    c_type = \"static inline REAL\"\n",
    "    name = \"FD1_arbitrary_upwind_x\"+str(dirn)+\"_dirn\"\n",
    "    params = \"\"\"const paramstruct *restrict params, const REAL *restrict gf,\n",
    "                                                const int i0,const int i1,const int i2, const int offset\"\"\"\n",
    "    body = r\"\"\"switch(offset) {\n",
    "\"\"\"\n",
    "    tmp_list = []\n",
    "    for offset in range(0, int(radiation_BC_FD_order / 2) + 1):\n",
    "        tmp_list.append(offset)\n",
    "        if offset > 0:\n",
    "            tmp_list.append(-offset)\n",
    "\n",
    "    for offset in tmp_list:\n",
    "        body += \"case \" + str(offset) + \":\\n\"\n",
    "        body += \"  return (\"\n",
    "        coeffs, indices = get_arb_offset_FD_coeffs_indices(radiation_BC_FD_order, offset, 1)\n",
    "        for i, coeff in enumerate(coeffs):\n",
    "            if coeff == 0:\n",
    "                continue  # skip this iteration if coeff=0\n",
    "            offset = str(indices[i])\n",
    "            if i > 0:\n",
    "                body += \"          \"\n",
    "            if offset == \"0\":\n",
    "                body += \"+\"+str(sp.ccode(coeff))+\"*gf[IDX3S(i0,i1,i2)]\\n\"\n",
    "            else:\n",
    "                if dirn == 0:\n",
    "                    body += \"+\"+str(sp.ccode(coeff))+\"*gf[IDX3S(i0+\"+offset+\",i1,i2)]\\n\"\n",
    "                elif dirn == 1:\n",
    "                    body += \"+\"+str(sp.ccode(coeff))+\"*gf[IDX3S(i0,i1+\"+offset+\",i2)]\\n\"\n",
    "                elif dirn == 2:\n",
    "                    body += \"+\"+str(sp.ccode(coeff))+\"*gf[IDX3S(i0,i1,i2+\"+offset+\")]\\n\"\n",
    "        body = body[:-1].replace(\"+-\", \"-\") + \") * invdx\"+str(dirn)+\";\\n\"\n",
    "    body += \"\"\"}\n",
    "return 0.0 / 0.0;  // poison output if offset computed incorrectly\n",
    "\"\"\"\n",
    "    rel_path_to_Cparams = os.path.join(\".\")\n",
    "\n",
    "    _prototype, func = Cfunction(includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=indent_Ccode(body, \"  \"),\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)\n",
    "\n",
    "    par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", default_FDORDER)\n",
    "\n",
    "    return func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='radiation_bcs_numerical_compute_partial_r_f'></a>\n",
    "\n",
    "### Step 6.b.v: The $\\partial_r f$ term: Putting it all together in `compute_partial_r_f()` \\[Back to [top](#toc)\\]\n",
    "$$\\label{radiation_bcs_numerical_compute_partial_r_f}$$\n",
    "\n",
    "The NRPy+ generated inlined C function `compute_partial_r_f()` evaluates $\\partial_r f$ via\n",
    "\n",
    "$$\n",
    "\\partial_r f = \\frac{\\partial x^i}{\\partial r}\\partial_i f,\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "* $\\frac{\\partial x^i}{\\partial r}$ is computed with C function `r_and_partial_xi_partial_r_derivs()`, and\n",
    "* $\\partial_i f$ is computed with C function `FD1_arbitrary_upwind_xN_dirn()`;\n",
    "\n",
    "both of these C functions are constructed above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_Cfunction_compute_partial_r_f(radiation_BC_FD_order=-1):\n",
    "    desc = \"Compute \\partial_r f\"\n",
    "    c_type = \"static inline REAL\"\n",
    "    name = \"compute_partial_r_f\"\n",
    "    params = \"\"\"const paramstruct *restrict params, REAL *restrict xx[3], const REAL *restrict gfs,\n",
    "                                       const int which_gf, const int dest_i0,const int dest_i1,const int dest_i2,\n",
    "                                       const int FACEi0,const int FACEi1,const int FACEi2,\n",
    "                                       const REAL partial_x0_partial_r, const REAL partial_x1_partial_r, const REAL partial_x2_partial_r\"\"\"\n",
    "    Jac_dUSph_dDrfmUD, Jac_dUrfm_dDSphUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Spherical()\n",
    "\n",
    "    default_FDORDER = par.parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\")\n",
    "    if radiation_BC_FD_order == -1:\n",
    "        radiation_BC_FD_order = default_FDORDER\n",
    "\n",
    "    body = r\"\"\"  ///////////////////////////////////////////////////////////\n",
    "\n",
    "  // FD1_stencil_radius = radiation_BC_FD_order/2 = \"\"\" + str(int(radiation_BC_FD_order/2)) + r\"\"\"\n",
    "  const int FD1_stencil_radius = \"\"\" + str(int(radiation_BC_FD_order/2)) + r\"\"\";\n",
    "\n",
    "  const int ntot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\n",
    "\n",
    "  ///////////////////////////////////////////////////////////\n",
    "  // Next we'll compute partial_xi f, using a maximally-centered stencil.\n",
    "  //   The {i0,i1,i2}_offset parameters set the offset of the maximally-centered\n",
    "  //   stencil, such that an offset=0 implies a centered stencil.\n",
    "\n",
    "  // CHECK: Nxx_plus_2NGHOSTS0=10; FD1_stencil_radius=2. Then Nxx_plus_2NGHOSTS0-FD1_stencil_radius-1 = 7\n",
    "  //  if dest_i0 = 9, we get i0_offset=7-9=-2, so the (4th order) deriv\n",
    "  //  stencil is: -4,-3,-2,-1,0\n",
    "\n",
    "  // CHECK: if FD1_stencil_radius=2 and dest_i0 = 1, we get i0_offset = FD1_stencil_radius-dest_i0 = 1,\n",
    "  //  so the (4th order) deriv stencil is: -1,0,1,2,3\n",
    "\n",
    "  // CHECK: if FD1_stencil_radius=2 and dest_i0 = 0, we get i0_offset = FD1_stencil_radius-1 = 2,\n",
    "  //  so the (4th order) deriv stencil is: 0,1,2,3,4\n",
    "\"\"\"\n",
    "    for i in range(3):\n",
    "        si = str(i)\n",
    "        if check_zero(Jac_dUrfm_dDSphUD[i][0]):\n",
    "            body += \"  const REAL partial_x\"+si+\"_f=0.0;\\n\"\n",
    "        else:\n",
    "            body += \"  int i\"+si+\"_offset = FACEi\"+si+\";  // Shift stencil away from the face we're updating.\\n\"\n",
    "            body += \"  // Next adjust i\"+si+\"_offset so that FD stencil never goes out of bounds.\\n\"\n",
    "            body += \"  if(dest_i\"+si+\" < FD1_stencil_radius) i\"+si+\"_offset = FD1_stencil_radius-dest_i\"+si+\";\\n\"\n",
    "            body += \"  else if(dest_i\"+si+\" > (Nxx_plus_2NGHOSTS\"+si+\"-FD1_stencil_radius-1)) i\"+si+\"_offset = (Nxx_plus_2NGHOSTS\"+si+\"-FD1_stencil_radius-1) - dest_i\"+si+\";\\n\"\n",
    "            body += \"  const REAL partial_x\"+si+\"_f=FD1_arbitrary_upwind_x\"+si+\"_dirn(params,&gfs[which_gf*ntot],dest_i0,dest_i1,dest_i2,i\"+si+\"_offset);\\n\"\n",
    "    body += \"  return partial_x0_partial_r*partial_x0_f + partial_x1_partial_r*partial_x1_f + partial_x2_partial_r*partial_x2_f;\\n\"\n",
    "    rel_path_to_Cparams = os.path.join(\".\")\n",
    "\n",
    "    _prototype, func = Cfunction(\n",
    "        includes=[],   desc=desc, c_type=c_type, name=name, params=params,    body=body,\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)\n",
    "    return func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='radiation_bcs_evaluating_k'></a>\n",
    "\n",
    "### Step 6.b.vi: Evaluating the $k$ in the $k/r^3$ term \\[Back to [top](#toc)\\]\n",
    "$$\\label{radiation_bcs_evaluating_k}$$\n",
    "\n",
    "First we note that in the *ansatz*, if $f(r,t)$ perfectly satisfied the outgoing wave equation, then $k=0$ and\n",
    "\n",
    "$$\n",
    "\\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave} = -\\frac{c}{r} \\left[r \\partial_r f + \\left(f - f_{r\\to\\infty}\\right)\\right].\n",
    "$$\n",
    "\n",
    "However, numerical error, errors associated with outer boundaries being placed at finite (as opposed to infinite) radii, and higher-order radial falloffs will invariably lead to $k$ being *nonzero*. \n",
    "\n",
    "Given that $\\partial_t f$ is evaluated *at all points in the interior* of the grid, the difference between perfect satisfaction of the radially outgoing wave equation and the actual solution is *known* at these points. Let's call this difference $\\xi$:\n",
    "\n",
    "$$\n",
    "\\xi = \\partial_t f - \\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave} \\equiv \\frac{k}{r^3}.\n",
    "$$\n",
    "\n",
    "We compute $\\xi_{\\rm int} = \\xi(r_{\\rm int})$ at a neighboring interior point at radius $r=r_{\\rm int}$, so that\n",
    "\n",
    "$$\n",
    "\\xi_{\\rm int} = [\\partial_t f]_{\\rm int} - \\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave,\\ int} \\equiv \\frac{k}{r_{\\rm int}^3}.\n",
    "$$\n",
    "\n",
    "In this way, we obtain $k$:\n",
    "\n",
    "$$\n",
    "k = r_{\\rm int}^3 \\left([\\partial_t f]_{\\rm int} - \\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave,\\ int}\\right)\n",
    "$$\n",
    "\n",
    "To determine the appropriate interior point, we simply keep track of the current boundary face being updated and choose the nearest neighbor in the direction of the grid center (assumed to be at $r=0$).\n",
    "\n",
    "As computing $k$ depends on all other functions, we incorporate its computation into the core `apply_bcs_outerradiation_and_inner()` routine - the subject of the next subsection."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='radiation_bcs_apply_bcs_radiation'></a>\n",
    "\n",
    "### Step 6.b.vii: Putting it all together: `radiation_bcs_single_pt()` \\[Back to [top](#toc)\\]\n",
    "$$\\label{radiation_bcs_apply_bcs_radiation}$$\n",
    "\n",
    "`radiation_bcs_single_pt()` combines all the above algorithms to evaluate radiation boundary conditions at a single point.\n",
    "\n",
    "Recall that at each point on the boundary, we have\n",
    "$$\n",
    "\\partial_t f = \\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave} + \\frac{k}{r^3},\n",
    "$$\n",
    "where\n",
    "$$\n",
    "\\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave} = -\\frac{c}{r} \\left[r \\partial_r f + \\left(f - f_{r\\to\\infty}\\right)\\right],\n",
    "$$\n",
    "and\n",
    "$$\n",
    "k = r_{\\rm int}^3 \\left([\\partial_t f]_{\\rm int} - \\left[\\partial_t f\\right]_{\\rm Outgoing\\ Wave,\\ int}\\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_Cfunction_radiation_bcs(radiation_BC_FD_order=-1):\n",
    "    includes = []\n",
    "    prefunc = \"\"\n",
    "    Jac_dUSph_dDrfmUD, Jac_dUrfm_dDSphUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Spherical()\n",
    "    for i in range(3):\n",
    "        # Do not generate FD1_arbitrary_upwind_xj_dirn() if the symbolic expression for dxj/dr == 0!\n",
    "        if not check_zero(Jac_dUrfm_dDSphUD[i][0]):\n",
    "            prefunc += setup_Cfunction_FD1_arbitrary_upwind(dirn=i, radiation_BC_FD_order=radiation_BC_FD_order)\n",
    "    prefunc += setup_Cfunction_r_and_partial_xi_partial_r_derivs()\n",
    "    prefunc += setup_Cfunction_compute_partial_r_f(radiation_BC_FD_order=radiation_BC_FD_order)\n",
    "    desc = r\"\"\"*** Apply radiation BCs to all outer boundaries. ***\n",
    "\"\"\"\n",
    "    c_type = \"static inline REAL\"\n",
    "    name = \"radiation_bcs\"\n",
    "    params = \"\"\"const paramstruct *restrict params, const bc_struct *restrict bcstruct,REAL *restrict xx[3],\n",
    "                                 const REAL *restrict gfs, REAL *restrict gfs_rhss,\n",
    "                                 const int which_gf, const REAL gf_wavespeed, const REAL gf_f_infinity,\n",
    "                                 const int dest_i0,const int dest_i1,const int dest_i2,\n",
    "                                 const short FACEi0,const short FACEi1,const short FACEi2\"\"\"\n",
    "    body = r\"\"\"// Nearest \"interior\" neighbor of this gridpoint, based on current face\n",
    "const int dest_i0_int=dest_i0+1*FACEi0, dest_i1_int=dest_i1+1*FACEi1, dest_i2_int=dest_i2+1*FACEi2;\n",
    "REAL r, partial_x0_partial_r,partial_x1_partial_r,partial_x2_partial_r;\n",
    "REAL r_int, partial_x0_partial_r_int,partial_x1_partial_r_int,partial_x2_partial_r_int;\n",
    "r_and_partial_xi_partial_r_derivs(params,xx[0][dest_i0],xx[1][dest_i1],xx[2][dest_i2],\n",
    "                                  &r, &partial_x0_partial_r, &partial_x1_partial_r,  &partial_x2_partial_r);\n",
    "r_and_partial_xi_partial_r_derivs(params, xx[0][dest_i0_int], xx[1][dest_i1_int], xx[2][dest_i2_int],\n",
    "                                  &r_int, &partial_x0_partial_r_int, &partial_x1_partial_r_int, &partial_x2_partial_r_int);\n",
    "const REAL partial_r_f     = compute_partial_r_f(params,xx,gfs, which_gf,dest_i0,    dest_i1,    dest_i2,\n",
    "                                                 FACEi0,FACEi1,FACEi2,\n",
    "                                                 partial_x0_partial_r    ,partial_x1_partial_r    ,partial_x2_partial_r);\n",
    "const REAL partial_r_f_int = compute_partial_r_f(params,xx,gfs, which_gf,dest_i0_int,dest_i1_int,dest_i2_int,\n",
    "                                                 FACEi0,FACEi1,FACEi2,\n",
    "                                                 partial_x0_partial_r_int,partial_x1_partial_r_int,partial_x2_partial_r_int);\n",
    "\n",
    "const int idx3 = IDX3S(dest_i0,dest_i1,dest_i2);\n",
    "const int idx3_int = IDX3S(dest_i0_int,dest_i1_int,dest_i2_int);\n",
    "\n",
    "const REAL partial_t_f_int = gfs_rhss[IDX4ptS(which_gf, idx3_int)];\n",
    "\n",
    "const REAL c = gf_wavespeed;\n",
    "const REAL f_infinity = gf_f_infinity;\n",
    "const REAL f     = gfs[IDX4ptS(which_gf, idx3)];\n",
    "const REAL f_int = gfs[IDX4ptS(which_gf, idx3_int)];\n",
    "const REAL partial_t_f_int_outgoing_wave = -c * (partial_r_f_int + (f_int - f_infinity) / r_int);\n",
    "\n",
    "const REAL k = r_int*r_int*r_int * (partial_t_f_int - partial_t_f_int_outgoing_wave);\n",
    "\n",
    "const REAL rinv = 1.0 / r;\n",
    "const REAL partial_t_f_outgoing_wave = -c * (partial_r_f + (f - f_infinity) * rinv);\n",
    "\n",
    "return partial_t_f_outgoing_wave + k * rinv*rinv*rinv;\n",
    "\"\"\"\n",
    "    rel_path_to_Cparams = os.path.join(\".\")\n",
    "\n",
    "    _prototype, func = Cfunction(\n",
    "        includes=includes, prefunc=prefunc,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=indent_Ccode(body, \"  \"),\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)\n",
    "    return func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='apply_bcs_outerradiation_and_inner'></a>\n",
    "\n",
    "### Step 6.b.viii: `apply_bcs_outerradiation_and_inner()`: Apply radiation BCs at all outer boundary points and inner BCs at all inner boundary points \\[Back to [top](#toc)\\]\n",
    "$$\\label{apply_bcs_outerradiation_and_inner}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_apply_bcs_outerradiation_and_inner(rel_path_to_Cparams=os.path.join(\".\"), radiation_BC_FD_order=2):\n",
    "    includes = [os.path.join(rel_path_to_Cparams, \"NRPy_basic_defines.h\"),\n",
    "                os.path.join(rel_path_to_Cparams, \"NRPy_function_prototypes.h\")]\n",
    "    prefunc = setup_Cfunction_radiation_bcs(radiation_BC_FD_order=radiation_BC_FD_order)\n",
    "    desc = \"\"\n",
    "    c_type = \"void\"\n",
    "    name = \"apply_bcs_outerradiation_and_inner\"\n",
    "    params = \"\"\"const paramstruct *restrict params, const bc_struct *restrict bcstruct, REAL *restrict xx[3],\n",
    "                                        const REAL custom_wavespeed[NUM_EVOL_GFS],\n",
    "                                        const REAL custom_f_infinity[NUM_EVOL_GFS],\n",
    "                                        REAL *restrict gfs, REAL *restrict rhs_gfs\"\"\"\n",
    "    body = r\"\"\"\n",
    "  // Unpack bc_info from bcstruct\n",
    "  const bc_info_struct *bc_info = &bcstruct->bc_info;\n",
    "\n",
    "  ////////////////////////////////////////////////////////\n",
    "  // STEP 1 of 2: Apply BCs to pure outer boundary points.\n",
    "  //              By \"pure\" we mean that these points are\n",
    "  //              on the outer boundary and not also on\n",
    "  //              an inner boundary.\n",
    "  //              Here we fill in the innermost ghost zone\n",
    "  //              layer first and move outward. At each\n",
    "  //              layer, we fill in +/- x0 faces first,\n",
    "  //              then +/- x1 faces, finally +/- x2 faces,\n",
    "  //              filling in the edges as we go.\n",
    "  // Spawn N OpenMP threads, either across all cores, or according to e.g., taskset.\n",
    "#pragma omp parallel\n",
    "  {\n",
    "    for(int which_gz=0;which_gz<NGHOSTS;which_gz++) for(int dirn=0;dirn<3;dirn++) {\n",
    "        // This option results in about 1.6% slower runtime for SW curvilinear at 64x24x24 on 8-core Ryzen 9 4900HS\n",
    "        //#pragma omp for collapse(2)\n",
    "        //for(int which_gf=0;which_gf<NUM_EVOL_GFS;which_gf++) for(int idx2d=0;idx2d<bc_info->num_pure_outer_boundary_points[which_gz][dirn];idx2d++) {\n",
    "        //  {\n",
    "        // Don't spawn a thread if there are no boundary points to fill; results in a nice little speedup.\n",
    "        if(bc_info->num_pure_outer_boundary_points[which_gz][dirn] > 0) {\n",
    "#pragma omp for  // threads have been spawned; here we distribute across them\n",
    "          for(int idx2d=0;idx2d<bc_info->num_pure_outer_boundary_points[which_gz][dirn];idx2d++) {\n",
    "            const short i0 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i0;\n",
    "            const short i1 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i1;\n",
    "            const short i2 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].i2;\n",
    "            const short FACEX0 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX0;\n",
    "            const short FACEX1 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX1;\n",
    "            const short FACEX2 = bcstruct->pure_outer_bc_array[dirn + (3*which_gz)][idx2d].FACEX2;\n",
    "            const int idx3 = IDX3S(i0,i1,i2);\n",
    "            for(int which_gf=0;which_gf<NUM_EVOL_GFS;which_gf++) {\n",
    "              // *** Apply radiation BCs to all outer boundary points. ***\n",
    "              rhs_gfs[IDX4ptS(which_gf, idx3)] = radiation_bcs(params, bcstruct, xx, gfs, rhs_gfs, which_gf,\n",
    "                                                               custom_wavespeed[which_gf], custom_f_infinity[which_gf],\n",
    "                                                               i0,i1,i2, FACEX0,FACEX1,FACEX2);\n",
    "            }\n",
    "          }\n",
    "        }\n",
    "      }\n",
    "  }\n",
    "\n",
    "  ///////////////////////////////////////////////////////\n",
    "  // STEP 2 of 2: Apply BCs to inner boundary points.\n",
    "  //              These map to either the grid interior\n",
    "  //              (\"pure inner\") or to pure outer boundary\n",
    "  //              points (\"inner maps to outer\"). Those\n",
    "  //              that map to outer require that outer be\n",
    "  //              populated first; hence this being\n",
    "  //              STEP 2 OF 2.\n",
    "  apply_bcs_inner_only(params, bcstruct, rhs_gfs); // <- apply inner BCs to RHS gfs only\n",
    "\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        prefunc=prefunc,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='start2finish'></a>\n",
    "\n",
    "# Step 7: `CurviBC_Playground.c`: Start-to-Finish C code module for testing & validating curvilinear boundary conditions \\[Back to [top](#toc)\\]\n",
    "$$\\label{start2finish}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='register_gfs'></a>\n",
    "\n",
    "## Step 7.a: Register gridfunctions of all 10 parity types \\[Back to [top](#toc)\\]\n",
    "$$\\label{register_gfs}$$ \n",
    "\n",
    "Here we register within NRPy+ one gridfunction per each of the 10 parity conditions. These will be output automatically to `NRPy_basic_defines.h` below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:15.269753Z",
     "iopub.status.busy": "2021-09-08T15:47:15.269360Z",
     "iopub.status.idle": "2021-09-08T15:47:15.271299Z",
     "shell.execute_reply": "2021-09-08T15:47:15.270905Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 7.a: Register gridfunctions of all 10 parity types\n",
    "\n",
    "# 6 gridfunctions, corresponding to all unique rank-2 tensor components:\n",
    "ranktwosymmDD = ixp.register_gridfunctions_for_single_rank2(\"EVOL\",\"ranktwosymmDD\", \"sym01\", f_infinity=0.0, wavespeed=1.0)\n",
    "# 3 gridfunctions, corresponding to all unique rank-1 tensor components:\n",
    "rankoneU = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"rankoneU\", f_infinity=0.0, wavespeed=1.0)\n",
    "# 1 rank-0 (scalar) gridfunction\n",
    "# WAS: rankzero = ixp.gri.register_gridfunctions(\"EVOL\",\"rankzero\", f_infinity=1.0, wavespeed=sp.sqrt(2.0))\n",
    "# NOW: consistent with ScalarWaveCurvilinear:\n",
    "rankzero = ixp.gri.register_gridfunctions(\"EVOL\",\"rankzero\", f_infinity=1.0, wavespeed=1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='validate'></a>\n",
    "\n",
    "## Step 7.b: Set up test data for Curvilinear Boundary Conditions code validation \\[Back to [top](#toc)\\]\n",
    "$$\\label{validate}$$\n",
    "\n",
    "We will validate this curvilinear boundary condition module by comparing its results with the original (trusted) SENR code, as follows:\n",
    "\n",
    "* **Discrete data test**:\n",
    "    1. Fill all 10 gridfunctions at each gridpoint with the unique gridpoint integer index `IDX3S(i0,i1,i2)`\n",
    "    1. Apply curvilinear boundary conditions\n",
    "    1. Compare output data at all gridpoints with those from the original SENR code. Agreement should be perfect.\n",
    "\n",
    "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:\n",
    "\n",
    "* **Smooth data test** (TODO):\n",
    "    1. Fill all 10 gridfunctions with data that are smooth in the Cartesian basis.\n",
    "    1. Apply Jacobian transformation to all data points, to convert to curvilinear basis\n",
    "    1. Apply curvilinear boundary conditions\n",
    "    1. Apply Jacobian transformation to all data points, to convert back to Cartesian basis\n",
    "    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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:15.274546Z",
     "iopub.status.busy": "2021-09-08T15:47:15.274158Z",
     "iopub.status.idle": "2021-09-08T15:47:15.276042Z",
     "shell.execute_reply": "2021-09-08T15:47:15.275654Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_CurviBC_Discrete_initial_data():\n",
    "    add_to_Cfunction_dict(\n",
    "        includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"],\n",
    "        desc     =\"Populate in_gfs[] with discrete, regular data\",\n",
    "        c_type     =\"void\",\n",
    "        name     =\"CurviBC_Discrete_initial_data\",\n",
    "        params   =\"const paramstruct *restrict params,REAL *restrict in_gfs\",\n",
    "        body     =\"\"\"\n",
    "#pragma omp parallel for\n",
    "    for(int i=0;i<Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2*NUM_EVOL_GFS;i++) {\n",
    "        in_gfs[i] = (REAL)i;\n",
    "    }\n",
    "\"\"\", rel_path_to_Cparams = os.path.join(\".\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "import loop as lp\n",
    "import ScalarWave.InitialData as swid\n",
    "def add_to_Cfunction_dict_CurviBC_Smooth_ScalarWave_initial_data():\n",
    "    swid.SphericalGaussian(CoordSystem=CoordSystem)\n",
    "    add_to_Cfunction_dict(\n",
    "        includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"],\n",
    "        desc     =\"Populate in_gfs[] and in_gfs_rhss[] with smooth, regular data\",\n",
    "        c_type     =\"void\",\n",
    "        name     =\"CurviBC_Smooth_ScalarWave_initial_data\",\n",
    "        params   =\"\"\"const paramstruct *restrict params,REAL *restrict xx[3],\n",
    "                                            REAL *restrict in_gfs,REAL *restrict in_gfs_rhss\"\"\",\n",
    "        body    = \"\"\"\n",
    "  const REAL sigma = 3.0;\n",
    "  const REAL time = 0.0;\n",
    "\"\"\" + lp.simple_loop('AllPoints, Read_xxs',\n",
    "                                  outputC([swid.uu_ID,swid.vv_ID],\n",
    "              [\"in_gfs[IDX4S(RANKZEROGF,i0,i1,i2)]\",\"in_gfs_rhss[IDX4S(RANKZEROGF,i0,i1,i2)]\"],\n",
    "               filename=\"returnstring\", params=\"outCverbose=False,includebraces=False\")),\n",
    "        rel_path_to_Cparams = os.path.join(\".\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mainc'></a>\n",
    "\n",
    "## Step 7.c: `CurviBC_Playground`'s `main.c` Code \\[Back to [top](#toc)\\]\n",
    "$$\\label{mainc}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:15.281826Z",
     "iopub.status.busy": "2021-09-08T15:47:15.281400Z",
     "iopub.status.idle": "2021-09-08T15:47:15.283297Z",
     "shell.execute_reply": "2021-09-08T15:47:15.282912Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_main__CurviBC_Playground():\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    desc = \"\"\"// main() function:\n",
    "// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates\n",
    "// Step 1: Write test data to gridfunctions\n",
    "// Step 2: Overwrite all data in ghost zones with NaNs\n",
    "// Step 3: Apply curvilinear boundary conditions\n",
    "// Step 4: Print gridfunction data after curvilinear boundary conditions have been applied\n",
    "// Step 5: Free all allocated memory\n",
    "\"\"\"\n",
    "    c_type = \"int\"\n",
    "    name = \"main\"\n",
    "    params = \"int argc, const char *argv[]\"\n",
    "    body  = r\"\"\"  // Step 0a: Read command-line input, error out if nonconformant\n",
    "  if(argc != 2) {\n",
    "    fprintf(stderr,\"Error: Expected one command-line argument: ./CurviBC_Playground [test type: Smooth or Discrete],\\n\");\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  griddata_struct griddata;\n",
    "  set_Cparameters_to_default(&griddata.params);\n",
    "\n",
    "  char CoordSystem_name[50];\n",
    "  snprintf(CoordSystem_name, 50, \"\"\"\n",
    "    body += \"\\\"\"+CoordSystem+\"\\\"\"\n",
    "    body += r\"\"\");\n",
    "\n",
    "  // Step 0b: Set number of gridpoints...\n",
    "  const int Nxx[3] = { 4, 4, 4 };\n",
    "\n",
    "  // Step 0c: Set test type to Smooth or Discrete\n",
    "  char test_type[100];\n",
    "  snprintf(test_type, 100, \"%s\", argv[1]);\n",
    "  if(strncmp(\"Smooth\",  test_type, 100) != 0 &&\n",
    "     strncmp(\"Discrete\",test_type, 100) != 0) {\n",
    "    fprintf(stderr,\"Error: test type = %s not supported. Choose Smooth or Discrete (CASE SENSITIVE).\\n\",test_type);\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  // Step 0d: Uniform coordinate grids are stored to *xx[3]\n",
    "  // Step 0d.i: Set bcstruct\n",
    "  {\n",
    "    int EigenCoord;\n",
    "    EigenCoord = 1;\n",
    "    // Step 0d.ii: Call set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
    "    //             params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
    "    //             chosen Eigen-CoordSystem.\n",
    "    set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &griddata.params, griddata.xx);\n",
    "    // Step 0e: Find ghostzone mappings; set up bcstruct\n",
    "    bcstruct_set_up(&griddata.params, griddata.xx, &griddata.bcstruct);\n",
    "    // Step 0e.i: Free allocated space for xx[][] array\n",
    "    for(int i=0;i<3;i++) free(griddata.xx[i]);\n",
    "\n",
    "    // Step 0f: Call set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
    "    //          params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
    "    //          chosen (non-Eigen) CoordSystem.\n",
    "    EigenCoord = 0;\n",
    "    set_Nxx_dxx_invdx_params__and__xx(EigenCoord, Nxx, &griddata.params, griddata.xx);\n",
    "  }\n",
    "\n",
    "  // Step 0g: Set all C parameters \"blah\" for params.blah, including\n",
    "  //          Nxx_plus_2NGHOSTS0 = params.Nxx_plus_2NGHOSTS0, etc.\n",
    "  const int Nxx_plus_2NGHOSTS0 = griddata.params.Nxx_plus_2NGHOSTS0;\n",
    "  const int Nxx_plus_2NGHOSTS1 = griddata.params.Nxx_plus_2NGHOSTS1;\n",
    "  const int Nxx_plus_2NGHOSTS2 = griddata.params.Nxx_plus_2NGHOSTS2;\n",
    "\n",
    "  // Step 0h: Allocate memory for gridfunctions and store to local pointers test_gfs & test_gfs_rhss\n",
    "  griddata.test_gfs      = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2);\n",
    "  griddata.test_gfs_rhss = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2);\n",
    "  REAL *restrict test_gfs      = griddata.test_gfs;\n",
    "  REAL *restrict test_gfs_rhss = griddata.test_gfs_rhss;\n",
    "\n",
    "  // Step 1: Write test data to gridfunctions\n",
    "  if(strncmp(\"Discrete\", test_type, 100)==0) {\n",
    "    CurviBC_Discrete_initial_data(&griddata.params, test_gfs);\n",
    "    CurviBC_Discrete_initial_data(&griddata.params, test_gfs_rhss);\n",
    "  } else {\n",
    "    fprintf(stderr, \"Sorry, curvilinear boundary conditions test = %s not yet supported. Feel free to contribute!\\n\",test_type);\n",
    "  }\n",
    "\n",
    "  // Step 2: Overwrite all data in ghost zones with NaNs\n",
    "  LOOP_REGION(0,Nxx_plus_2NGHOSTS0, 0,Nxx_plus_2NGHOSTS1, 0,Nxx_plus_2NGHOSTS2) {\n",
    "    for(int gf=0;gf<NUM_EVOL_GFS;gf++) {\n",
    "      const int idx4 = IDX4S(gf,i0,i1,i2);\n",
    "      if(i0 < NGHOSTS || i0 >= Nxx_plus_2NGHOSTS0-NGHOSTS) test_gfs_rhss[idx4] = +(0.0 / 0.0);\n",
    "      if(i1 < NGHOSTS || i1 >= Nxx_plus_2NGHOSTS1-NGHOSTS) test_gfs_rhss[idx4] = +(0.0 / 0.0);\n",
    "      if(i2 < NGHOSTS || i2 >= Nxx_plus_2NGHOSTS2-NGHOSTS) test_gfs_rhss[idx4] = +(0.0 / 0.0);\n",
    "    }\n",
    "  }\n",
    "\n",
    "  // Step 3: Apply outer and inner boundary conditions\n",
    "  if(griddata.params.outer_bc_type == EXTRAPOLATION_OUTER_BCS) {\n",
    "    // Warning: Extrapolation BCs should generally be applied to test_gfs, not test_gfs_rhss.\n",
    "    //          We apply them to rhss here just to make the coding simpler (gf_rhss have the\n",
    "    //          same data in the Discrete test as gfs).\n",
    "    apply_bcs_outerextrap_and_inner(&griddata.params, &griddata.bcstruct, test_gfs_rhss);\n",
    "  } else if(griddata.params.outer_bc_type == RADIATION_OUTER_BCS) {\n",
    "    apply_bcs_outerradiation_and_inner(&griddata.params, &griddata.bcstruct, griddata.xx,\n",
    "                                       gridfunctions_wavespeed,gridfunctions_f_infinity,\n",
    "                                       test_gfs,test_gfs_rhss);\n",
    "  }\n",
    "\n",
    "  // Step 4: Print gridfunction data after curvilinear boundary conditions have been applied:\n",
    "  char filename[120];  sprintf(filename,\"out4x4x4-%s-NGHOSTS4oFD.txt\", CoordSystem_name);\n",
    "  FILE *outfile = fopen(filename, \"w\");\n",
    "\n",
    "  LOOP_REGION(0,Nxx_plus_2NGHOSTS0, 0,Nxx_plus_2NGHOSTS1, 0,Nxx_plus_2NGHOSTS2) {\n",
    "    fprintf(outfile, \"%d %d %d | \", i0,i1,i2);\n",
    "    for(int gf=0;gf<NUM_EVOL_GFS;gf++) {\n",
    "      const int idx4 = IDX4S(gf,i0,i1,i2);\n",
    "      if(!isnan(test_gfs_rhss[idx4])) {\n",
    "        fprintf(outfile, \"%d \", (int)test_gfs_rhss[idx4]);\n",
    "      } else {\n",
    "        fprintf(stderr, \"ERROR: found NaN %d %d %d %d %d\\n\", gf, i0,i1,i2, NUM_EVOL_GFS);\n",
    "        //exit(1);\n",
    "      }\n",
    "    }\n",
    "    fprintf(outfile, \"\\n\");\n",
    "  }\n",
    "  fclose(outfile);\n",
    "\n",
    "  // Step 5: Free all allocated memory\n",
    "  free(griddata.bcstruct.inner_bc_array);\n",
    "  for(int ng=0;ng<NGHOSTS*3;ng++) free(griddata.bcstruct.pure_outer_bc_array[ng]);\n",
    "  free(griddata.test_gfs);  free(griddata.test_gfs_rhss);\n",
    "  for(int i=0;i<3;i++) free(griddata.xx[i]);\n",
    "  return 0;\n",
    "\"\"\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body, enableCparameters=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='curvibc_setupall'></a>\n",
    "\n",
    "## Step 7.d: `CurviBoundaryConditions_register_C_functions_and_NRPy_basic_defines()`: Add all CurviBC C codes to C function dictionary, and add CurviBC definitions to `NRPy_basic_defines.h` \\[Back to [top](#toc)\\]\n",
    "$$\\label{curvibc_setupall}$$\n",
    "\n",
    "Populate `NRPy_basic_defines.h` with `bcstruct` and friends, as well as parity types for registered gridfunctions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:15.286703Z",
     "iopub.status.busy": "2021-09-08T15:47:15.286318Z",
     "iopub.status.idle": "2021-09-08T15:47:15.288180Z",
     "shell.execute_reply": "2021-09-08T15:47:15.287794Z"
    }
   },
   "outputs": [],
   "source": [
    "# Only call this after ALL gridfunctions have been registered!\n",
    "def CurviBoundaryConditions_register_NRPy_basic_defines(verbose=True):\n",
    "    # Then set up the dictionary entry for CurviBC in NRPy_basic_defines\n",
    "    Nbd_str  = NRPy_basic_defines_CurviBC_data_structures()\n",
    "    Nbd_str += NRPy_basic_defines_set_gridfunction_defines_with_parity_types(verbose=verbose)\n",
    "    outC_NRPy_basic_defines_h_dict[\"CurviBoundaryConditions\"] = Nbd_str\n",
    "\n",
    "    # Register griddata_struct variables for this module,\n",
    "    #   where griddata_struct is declared in NRPy_basic_defines.h\n",
    "    gri.glb_griddata_struct_list += [gri.glb_griddata(__name__, \"bc_struct bcstruct;\")]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:15.286703Z",
     "iopub.status.busy": "2021-09-08T15:47:15.286318Z",
     "iopub.status.idle": "2021-09-08T15:47:15.288180Z",
     "shell.execute_reply": "2021-09-08T15:47:15.287794Z"
    }
   },
   "outputs": [],
   "source": [
    "def CurviBoundaryConditions_register_C_functions(rel_path_to_Cparams=os.path.join(\"./\"),\n",
    "                                                 radiation_BC_FD_order=4):\n",
    "    add_to_Cfunction_dict_bcstruct_set_up(rel_path_to_Cparams=rel_path_to_Cparams)\n",
    "    add_to_Cfunction_dict_apply_bcs_outerradiation_and_inner(rel_path_to_Cparams=rel_path_to_Cparams,\n",
    "                                                             radiation_BC_FD_order=radiation_BC_FD_order)\n",
    "    add_to_Cfunction_dict_apply_bcs_inner_only(rel_path_to_Cparams=rel_path_to_Cparams)\n",
    "    add_to_Cfunction_dict_apply_bcs_outerextrap_and_inner(rel_path_to_Cparams=rel_path_to_Cparams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='add_cfunction_dicts'></a>\n",
    "\n",
    "# Step 8: Add all C functions to `Cfunction_dict`, also output `NRPy_basic_defines.h` and `NRPy_function_prototypes.h` \\[Back to [top](#toc)\\]\n",
    "$$\\label{add_cfunction_dicts}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:15.295762Z",
     "iopub.status.busy": "2021-09-08T15:47:15.294840Z",
     "iopub.status.idle": "2021-09-08T15:47:18.611539Z",
     "shell.execute_reply": "2021-09-08T15:47:18.611295Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evolved gridfunction \"rankoneU0\" has parity type 1.\n",
      "Evolved gridfunction \"rankoneU1\" has parity type 2.\n",
      "Evolved gridfunction \"rankoneU2\" has parity type 3.\n",
      "Evolved gridfunction \"ranktwosymmDD00\" has parity type 4.\n",
      "Evolved gridfunction \"ranktwosymmDD01\" has parity type 5.\n",
      "Evolved gridfunction \"ranktwosymmDD02\" has parity type 6.\n",
      "Evolved gridfunction \"ranktwosymmDD11\" has parity type 7.\n",
      "Evolved gridfunction \"ranktwosymmDD12\" has parity type 8.\n",
      "Evolved gridfunction \"ranktwosymmDD22\" has parity type 9.\n",
      "Evolved gridfunction \"rankzero\" has parity type 0.\n"
     ]
    }
   ],
   "source": [
    "import outputC as outC\n",
    "# Then we set the coordinate system for the numerical grid\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\", CoordSystem)\n",
    "rfm.reference_metric()\n",
    "# Step 5: Generate & register set_Nxx_dxx_invdx_params__and__xx(), which sets\n",
    "#         params Nxx,Nxx_plus_2NGHOSTS,dxx,invdx, and xx[] for the\n",
    "#         chosen (not necessarily Eigen-) CoordSystem.\n",
    "rfm.register_NRPy_basic_defines(enable_rfm_precompute=False)\n",
    "rfm.add_to_Cfunc_dict_set_Nxx_dxx_invdx_params__and__xx()\n",
    "rfm.register_NRPy_basic_defines()\n",
    "\n",
    "CurviBoundaryConditions_register_NRPy_basic_defines()\n",
    "CurviBoundaryConditions_register_C_functions()\n",
    "\n",
    "add_to_Cfunction_dict_apply_bcs_inner_only()\n",
    "add_to_Cfunction_dict_apply_bcs_outerextrap_and_inner()\n",
    "add_to_Cfunction_dict_apply_bcs_outerradiation_and_inner(radiation_BC_FD_order=RADIATION_BC_FD_ORDER)\n",
    "\n",
    "\n",
    "outC.outputC_register_C_functions_and_NRPy_basic_defines()  # #define M_PI,  etc.\n",
    "# Declare paramstruct, register set_Cparameters_to_default(),\n",
    "#   and output declare_Cparameters_struct.h and set_Cparameters[].h:\n",
    "outC.NRPy_param_funcs_register_C_functions_and_NRPy_basic_defines(os.path.join(Ccodesrootdir))\n",
    "par.register_NRPy_basic_defines()  # add `paramstruct params` to griddata struct.\n",
    "\n",
    "gri.register_C_functions_and_NRPy_basic_defines(list_of_extras_in_griddata_struct=\n",
    "                                                [\"REAL *restrict test_gfs\",\"REAL *restrict test_gfs_rhss\"])  # #define IDX3S(),  etc.\n",
    "fin.register_C_functions_and_NRPy_basic_defines(NGHOSTS_account_for_onezone_upwind=True)  # #define NGHOSTS, etc.\n",
    "# Comment out above line and uncomment below line to confirm independent Python module agrees.\n",
    "# import CurviBoundaryConditions.CurviBoundaryConditions as CBC\n",
    "# CBC.CurviBoundaryConditions_register_C_functions_and_NRPy_basic_defines()\n",
    "\n",
    "# initial data function:\n",
    "add_to_Cfunction_dict_CurviBC_Discrete_initial_data()\n",
    "add_to_Cfunction_dict_CurviBC_Smooth_ScalarWave_initial_data()\n",
    "\n",
    "# main function:\n",
    "add_to_Cfunction_dict_main__CurviBC_Playground()\n",
    "\n",
    "outC.construct_NRPy_basic_defines_h(Ccodesrootdir)\n",
    "outC.construct_NRPy_function_prototypes_h(Ccodesrootdir)\n",
    "# The following is run from inside cmd.new_C_compile() in the next code cell.\n",
    "# outC.construct_Makefile_from_outC_function_dict(Ccodesrootdir, \"CurviBC_Playground\")\n",
    "# print(outC.outC_function_dict[\"set_Nxx_dxx_invdx_params__and__xx\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='senr_compare'></a>\n",
    "\n",
    "# Step 9: Validation: Compile & compare with original (trusted) SENR results \\[Back to [top](#toc)\\]\n",
    "$$\\label{senr_compare}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:18.614645Z",
     "iopub.status.busy": "2021-09-08T15:47:18.614349Z",
     "iopub.status.idle": "2021-09-08T15:47:19.240208Z",
     "shell.execute_reply": "2021-09-08T15:47:19.239773Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(EXEC): Executing `make -j18`...\n",
      "(BENCH): Finished executing in 0.41 seconds.\n",
      "Finished compilation.\n",
      "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./CurviBC_Playground Discrete`...\n",
      "(BENCH): Finished executing in 0.20 seconds.\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd\n",
    "# from outputC import construct_Makefile_from_outC_function_dict\n",
    "# construct_Makefile_from_outC_function_dict(Ccodesrootdir, \"CurviBC_Playground\", compiler_opt_option=\"fast\")\n",
    "cmd.new_C_compile(Ccodesrootdir, \"CurviBC_Playground\", compiler_opt_option=\"fast\") # fastdebug or debug also supported\n",
    "os.chdir(Ccodesrootdir)\n",
    "cmd.Execute(\"CurviBC_Playground\", \"Discrete\")\n",
    "#cmd.Execute(\"CurviBC_Playground\", \"4 4 4 Discrete\", \"out4x4x4-Spherical-NGHOSTS4oFD.txt\")\n",
    "os.chdir(\"..\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:19.244765Z",
     "iopub.status.busy": "2021-09-08T15:47:19.244251Z",
     "iopub.status.idle": "2021-09-08T15:47:19.246131Z",
     "shell.execute_reply": "2021-09-08T15:47:19.246437Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SinhSpherical boundary condition comparison test between this tutorial notebook & trusted original SENR code: PASSED\n"
     ]
    }
   ],
   "source": [
    "import filecmp\n",
    "if outer_bcs_type == \"extrapolation\":\n",
    "    if \"Cylindrical\" in CoordSystem:\n",
    "        if filecmp.cmp(os.path.join(Ccodesrootdir, \"out4x4x4-\"+CoordSystem+\"-NGHOSTS4oFD.txt\"),\n",
    "                       os.path.join(\"CurviBoundaryConditions\", \"SENRout4x4x4-Cylindrical_NGHOSTS4oFD.txt\")) == False:\n",
    "            print(\"ERROR: \"+CoordSystem+\" boundary conditions malfunction!\")\n",
    "            sys.exit(1)\n",
    "    elif \"Spherical\" in CoordSystem:\n",
    "        if filecmp.cmp(os.path.join(Ccodesrootdir, \"out4x4x4-\"+CoordSystem+\"-NGHOSTS4oFD.txt\"),\n",
    "                       os.path.join(\"CurviBoundaryConditions\", \"SENRout4x4x4-Spherical_NGHOSTS4oFD.txt\")) == False:\n",
    "            print(\"ERROR: \"+CoordSystem+\" boundary conditions malfunction!\")\n",
    "            sys.exit(1)\n",
    "    else:\n",
    "        print(\"ERROR: \"+CoordSystem+\" coordinate system comparison unavailable!\")\n",
    "        sys.exit(1)\n",
    "    print(CoordSystem + \" boundary condition comparison test between this tutorial notebook & trusted original SENR code: PASSED\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 10: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "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\n",
    "[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.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-09-08T15:47:19.249219Z",
     "iopub.status.busy": "2021-09-08T15:47:19.248835Z",
     "iopub.status.idle": "2021-09-08T15:47:22.888138Z",
     "shell.execute_reply": "2021-09-08T15:47:22.888544Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Start_to_Finish-Curvilinear_BCs.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-Start_to_Finish-Curvilinear_BCs.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Start_to_Finish-Curvilinear_BCs\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}