Tutorial-ETK_thorn-Interpolation_to_Arbitrary_Grids_multi_order.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",
    "# `interp_arbgrid_MO_ETK`: An Einstein Toolkit module for interpolation to arbitrary grids, at multiple interpolation orders, in Cartesian basis. \n",
    "\n",
    "## (Includes notes on transformations to other coordinate bases.)\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module is designed to interpolate arbitrary quantities on [Einstein Toolkit](https://einsteintoolkit.org/) Adaptive-Mesh Refinement (AMR) grids (using the [Carpet](https://carpetcode.org/) AMR infrastructure) to numerical grids with arbitrary sampling.\n",
    "\n",
    "**Validation Status:** <font color='red'><b> In progress</b></font>\n",
    "\n",
    "**Validation Notes:** This module is currently undergoing validation testing.\n",
    "\n",
    "\n",
    "## Introduction: \n",
    "\n",
    "Given some set of $N$ quantities $\\mathbf{Q}=\\{Q_0,Q_1,Q_2,...,Q_{N-2},Q_{N-1}\\}$, this module performs the following for each $Q_i$:\n",
    "\n",
    "1. Evaluate $Q_i$ at all gridpoints that are not ghost zones. Sometimes $Q_i$ is computed using finite difference derivatives, so this is necessary.\n",
    "1. Call upon Carpet's interpolation and interprocessor synchronization functions to fill in $Q_i$ at all ghost zones, *except* at the outer boundary. We do not generally trust $Q_i$ at the outer boundary due to errors associated with the approximate outer boundary conditions. \n",
    "1. At this point, $Q_i$ is set at all gridpoints except ghost zones at the outer boundary. Interpolate $Q_i$ to the desired output grids, **maintaining the Cartesian basis for all vectors and tensors**, and append the result to a file.\n",
    "\n",
    "This tutorial notebook takes a three-part structure. First, all the needed core Einstein Toolkit (ETK) C routines for interpolation are presented. Second, NRPy+ is used to output gridfunctions needed on the output grids. Third, the needed files for interfacing this module with the rest of the Einstein Toolkit (ccl files) are specified."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents: \n",
    "$$\\label{toc}$$\n",
    "\n",
    "1. [Step 1](#etkmodule): Setting up the Core C Code for the Einstein Toolkit Module\n",
    "    1. [Step 1.a](#etk_interp): Low-Level Einstein Toolkit Interpolation Function\n",
    "    1. [Step 1.b](#fileformat): Outputting to File\n",
    "    1. [Step 1.c](#maininterpolator): The Main Interpolator Driver Function\n",
    "    1. [Step 1.d](#standalonerandompoints): Standalone C code to output random points data \n",
    "1. [Step 2](#nrpy): Using NRPy+ to Generate C Code for Needed Gridfunctions\n",
    "    1. [Step 2.a](#nrpy_list_of_funcs_interp): Set up NRPy-based `list_of_functions_to_interpolate.h`\n",
    "        1. [Step 2.a.i](#nrpygrmhd): All GRMHD quantities, except vector potential $A_i$\n",
    "        1. [Step 2.a.ii](#unstaggera): Unstagger $A_i$ and add to \"list of functions to interpolate\"\n",
    "        1. [Step 2.a.iii](#nrpy4metric): Compute all 10 components of the 4-metric $g_{\\mu\\nu}$\n",
    "        1. [Step 2.a.iv](#nrpy4christoffels_cartesian):Compute all 40 4-Christoffels $\\Gamma^{\\mu}_{\\nu\\delta}$\n",
    "        1. [Step 2.a.v](#nrpy4christoffels_spherical):  Notes on computing all 40 4-Christoffels $\\Gamma^{\\mu}_{\\nu\\delta}$ in the Spherical basis\n",
    "        1. [Step 2.a.vi](#nrpybasisxform): Notes on basis transforming all Cartesian basis quantities to spherical\n",
    "        1. [Step 2.a.vii](#psi4andfriends): Output Weyl scalars $\\psi_0$ through $\\psi_4$, as well as Weyl invariants $J$ and $I$, from the `WeylScal4` ETK thorn\n",
    "        1. [Step 2.a.viii](#constraints_gfs): Hamiltonian and momentum constraints\n",
    "        1. [Step 2.a.ix](#nuc_eos_gfs): Nuclear (tabulated) equation of state quantities\n",
    "    1. [Step 2.b](#nrpy_c_calling_function): C code calling function for the NRPy+ C output\n",
    "    1. [Step 2.c](#nrpygetgfname): The `get_gf_name()` function\n",
    "    1. [Step 2.d](#nrpy_interp_counter): C Code for Initializing and incrementing `InterpCounter`\n",
    "    1. [Step 2.e](#validationagainstfm): Validation of interpolated data against exact Fishbone-Moncrief data\n",
    "1. [Step 3](#cclfiles): Interfacing with the rest of the Einstein Toolkit; Setting up CCL files\n",
    "    1. [Step 3.a](#makecodedefn): `make.code.defn`\n",
    "    1. [Step 3.b](#interfaceccl): `interface.ccl`\n",
    "    1. [Step 3.c](#paramccl): `param.ccl`\n",
    "    1. [Step 3.d](#scheduleccl): `schedule.ccl`\n",
    "1. [Step 4](#readingoutputfile): Python Script for Reading the Output File\n",
    "1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='etkmodule'></a>\n",
    "\n",
    "# Step 1: Setting up the Core C Code for the Einstein Toolkit Module \\[Back to [top](#toc)\\]\n",
    "$$\\label{etkmodule}$$\n",
    "\n",
    "First we set up the output directories for the ETK module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:10.855736Z",
     "iopub.status.busy": "2021-03-07T17:07:10.854802Z",
     "iopub.status.idle": "2021-03-07T17:07:10.862080Z",
     "shell.execute_reply": "2021-03-07T17:07:10.862737Z"
    }
   },
   "outputs": [],
   "source": [
    "import cmdline_helper as cmd  # NRPy+: Multi-platform Python command-line interface\n",
    "import shutil, os             # Standard Python modules for multiplatform OS-level functions, benchmarking\n",
    "\n",
    "# Create C code output directory:\n",
    "Ccodesdir = \"interp_arbgrid_MO_ETK\"\n",
    "# First remove C code output directory and all subdirectories if they exist\n",
    "# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
    "shutil.rmtree(Ccodesdir, ignore_errors=True)\n",
    "# Then create a fresh directory\n",
    "cmd.mkdir(Ccodesdir)\n",
    "cmd.mkdir(os.path.join(Ccodesdir,\"src/\"))\n",
    "cmd.mkdir(os.path.join(Ccodesdir,\"src\",\"standalone/\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='etk_interp'></a>\n",
    "\n",
    "## Step 1.a: Low-Level ETK Interpolation Function \\[Back to [top](#toc)\\]\n",
    "$$\\label{etk_interp}$$\n",
    "\n",
    "We start by writing the low-level interpolation function **`Interpolate_to_dest_grid()`**, which  to file. \n",
    "\n",
    "**`Interpolate_to_dest_grid()`** takes as input\n",
    "\n",
    "* **cctkGH**: Information about the underlying Cactus/Carpet grid hierarchy.\n",
    "* **interp_num_points**: Number of destination interpolation points\n",
    "* **point_x_temp, point_y_temp, point_z_temp**: Cartesian $(x,y,z)$ location for each of the **interp_num_points** interpolation points.\n",
    "* **input_array_names[1]**: List of input gridfunction names to interpolate. We will do this only one gridfunction at a time, for gridfunction $Q_i$, as described above.\n",
    "\n",
    "**`Interpolate_to_dest_grid()`** outputs:\n",
    "\n",
    "* **output_f[1]**: The gridfunction **input_array_names[1]** interpolated to the set of **interp_num_points** specified in the input."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:10.870699Z",
     "iopub.status.busy": "2021-03-07T17:07:10.869830Z",
     "iopub.status.idle": "2021-03-07T17:07:10.873357Z",
     "shell.execute_reply": "2021-03-07T17:07:10.873905Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/src/Interpolate_to_dest_grid.h\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/Interpolate_to_dest_grid.h\n",
    "\n",
    "void Interpolate_to_dest_grid(const cGH *restrict cctkGH,\n",
    "                              const CCTK_INT interp_num_points,\n",
    "                              const CCTK_INT interp_order,\n",
    "                              char interpolator_name[100],\n",
    "                              const CCTK_REAL *restrict point_x_temp,\n",
    "                              const CCTK_REAL *restrict point_y_temp,\n",
    "                              const CCTK_REAL *restrict point_z_temp,\n",
    "                              const CCTK_STRING input_array_names[1],\n",
    "                              CCTK_REAL *output_f[1]) {\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  CCTK_INT ierr;\n",
    "  const CCTK_INT NUM_INPUT_ARRAYS=1;\n",
    "  const CCTK_INT NUM_OUTPUT_ARRAYS=1;\n",
    "\n",
    "  CCTK_STRING coord_system = \"cart3d\";\n",
    "\n",
    "  // Set up handles\n",
    "  const CCTK_INT coord_system_handle = CCTK_CoordSystemHandle(coord_system);\n",
    "  if (coord_system_handle < 0) {\n",
    "    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,\n",
    "               \"can't get coordinate system handle for coordinate system \\\"%s\\\"!\",\n",
    "               coord_system);\n",
    "  }\n",
    "\n",
    "  const CCTK_INT operator_handle = CCTK_InterpHandle(interpolator_name);\n",
    "  if (operator_handle < 0)\n",
    "    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,\n",
    "               \"couldn't find interpolator \\\"%s\\\"!\",\n",
    "               interpolator_name);\n",
    "\n",
    "  char interp_order_string[10];\n",
    "  snprintf(interp_order_string, 10, \"order=%d\", interp_order);\n",
    "  CCTK_STRING interpolator_pars = interp_order_string;\n",
    "  CCTK_INT param_table_handle = Util_TableCreateFromString(interpolator_pars);\n",
    "  if (param_table_handle < 0) {\n",
    "    CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,\n",
    "               \"bad interpolator parameter(s) \\\"%s\\\"!\",\n",
    "               interpolator_pars);\n",
    "  }\n",
    "\n",
    "  CCTK_INT operand_indices[NUM_INPUT_ARRAYS]; //NUM_OUTPUT_ARRAYS + MAX_NUMBER_EXTRAS];\n",
    "  for(int i = 0 ; i < NUM_INPUT_ARRAYS  ; i++) {\n",
    "    operand_indices[i] = i;\n",
    "  }\n",
    "  Util_TableSetIntArray(param_table_handle, NUM_OUTPUT_ARRAYS,\n",
    "                        operand_indices, \"operand_indices\");\n",
    "\n",
    "  CCTK_INT operation_codes[NUM_INPUT_ARRAYS];\n",
    "  for(int i = 0 ; i < NUM_INPUT_ARRAYS  ; i++) {\n",
    "    operation_codes[i] = 0;\n",
    "  }\n",
    "  Util_TableSetIntArray(param_table_handle, NUM_OUTPUT_ARRAYS,\n",
    "                        operation_codes, \"operation_codes\");\n",
    "\n",
    "  const void* interp_coords[3] = { (const void *) point_x_temp,\n",
    "                                   (const void *) point_y_temp,\n",
    "                                   (const void *) point_z_temp };\n",
    "\n",
    "  CCTK_INT input_array_indices[NUM_INPUT_ARRAYS];\n",
    "  for(int i = 0 ; i < NUM_INPUT_ARRAYS ; i++) {\n",
    "    input_array_indices[i] = CCTK_VarIndex(input_array_names[i]);\n",
    "    if(input_array_indices[i] < 0) {\n",
    "      CCTK_VWarn(0, __LINE__, __FILE__, CCTK_THORNSTRING,\n",
    "                 \"COULD NOT FIND VARIABLE '%s'.\",\n",
    "                 input_array_names[i]);\n",
    "      exit(1);\n",
    "    }\n",
    "  }\n",
    "\n",
    "  CCTK_INT output_array_types[NUM_OUTPUT_ARRAYS];\n",
    "  for(int i = 0 ; i < NUM_OUTPUT_ARRAYS ; i++) {\n",
    "    output_array_types[i] = CCTK_VARIABLE_REAL;\n",
    "  }\n",
    "\n",
    "  void* output_arrays[NUM_OUTPUT_ARRAYS] = { (void *) output_f[0] };\n",
    "\n",
    "  // actual interpolation call\n",
    "  ierr = CCTK_InterpGridArrays(cctkGH,\n",
    "                               3, // number of dimensions\n",
    "                               operator_handle,\n",
    "                               param_table_handle,\n",
    "                               coord_system_handle,\n",
    "                               interp_num_points,\n",
    "                               CCTK_VARIABLE_REAL,\n",
    "                               interp_coords,\n",
    "                               NUM_INPUT_ARRAYS, // Number of input arrays\n",
    "                               input_array_indices,\n",
    "                               NUM_OUTPUT_ARRAYS, // Number of output arrays\n",
    "                               output_array_types,\n",
    "                               output_arrays);\n",
    "\n",
    "  if (ierr<0) {\n",
    "    CCTK_WARN(1,\"interpolation screwed up\");\n",
    "    Util_TableDestroy(param_table_handle);\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  ierr = Util_TableDestroy(param_table_handle);\n",
    "  if (ierr != 0) {\n",
    "    CCTK_WARN(1,\"Could not destroy table\");\n",
    "    exit(1);\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:10.880829Z",
     "iopub.status.busy": "2021-03-07T17:07:10.879924Z",
     "iopub.status.idle": "2021-03-07T17:07:10.883201Z",
     "shell.execute_reply": "2021-03-07T17:07:10.883666Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/src/interpolate_set_of_points_in_file.h\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/interpolate_set_of_points_in_file.h\n",
    "\n",
    "#define ALLOCATE_2D_GENERIC(type,array,ni,nj) type **array=(type **)malloc(ni * sizeof(type *)); \\\n",
    "  for(int cc = 0; cc < ni; cc++)                 array[cc]=(type * )malloc(nj * sizeof(type));\n",
    "#define FREE_2D_GENERIC(type,array,ni,nj) for(int cc = 0; cc < ni;cc++) free((void *)array[cc]); \\\n",
    "  /**/                                                                  free((void *)array);\n",
    "#include \"output_to_file.h\"\n",
    "\n",
    "// Calls the above function and output_to_file().\n",
    "void interpolate_set_of_points_in_file(CCTK_ARGUMENTS,char filename_basename[100],char gf_name[100],char interpolator_name[100],int num_interp_orders,int *interp_orders_list) {\n",
    "\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS; // Needed for x_center,y_center,z_center\n",
    "  // Set up output array:\n",
    "  // The name of the input gridfunction is always \"interp_arbgrid_MO_ETK::interped_gf\":\n",
    "  const CCTK_STRING input_array_names[1] = { \"interp_arbgrid_MO_ETK::interped_gf\" };\n",
    "  CCTK_REAL *points_dest_grid_x,*points_dest_grid_y,*points_dest_grid_z; // Coordinates of points of destination grid\n",
    "  CCTK_REAL **output_f; // Output to be written to dataset, will be filled with NaNs at points out of bounds\n",
    "  // For benchmarking purposes:\n",
    "  time_t start_timer,end_timer;\n",
    "  time(&start_timer); // Resolution of one second...\n",
    "  CCTK_REAL time_in_seconds;\n",
    "\n",
    "  int num_dest_grid_points;\n",
    "  if(CCTK_MyProc(cctkGH)==0) {\n",
    "    // Step 1: Read list of desired interpolation destination points from file:\n",
    "\n",
    "    // Step 1.a: Read integer at top of file indicating number of points.\n",
    "    int num_dest_grid_pointsx,num_dest_grid_pointsy,num_dest_grid_pointsz;\n",
    "    char pointsx_filename[100]; snprintf(pointsx_filename,100,\"%s-x.dat\",filename_basename);\n",
    "    printf(\"Reading list of x data points from file %s...\\n\",pointsx_filename);\n",
    "    FILE *pointsx_file = fopen(pointsx_filename, \"rb\");\n",
    "    if(!pointsx_file) { printf(\"Error: Unable to open %s\\n\",pointsx_filename); exit(1); }\n",
    "    fread(&num_dest_grid_pointsx, sizeof(int), 1, pointsx_file);\n",
    "\n",
    "    char pointsy_filename[100]; snprintf(pointsy_filename,100,\"%s-y.dat\",filename_basename);\n",
    "    printf(\"Reading list of y data points from file %s...\\n\",pointsy_filename);\n",
    "    FILE *pointsy_file = fopen(pointsy_filename, \"rb\");\n",
    "    if(!pointsy_file) { printf(\"Error: Unable to open %s\\n\",pointsy_filename); exit(1); }\n",
    "    fread(&num_dest_grid_pointsy, sizeof(int), 1, pointsy_file);\n",
    "\n",
    "    char pointsz_filename[100]; snprintf(pointsz_filename,100,\"%s-z.dat\",filename_basename);\n",
    "    printf(\"Reading list of z data points from file %s...\\n\",pointsz_filename);\n",
    "    FILE *pointsz_file = fopen(pointsz_filename, \"rb\");\n",
    "    if(!pointsz_file) { printf(\"Error: Unable to open %s\\n\",pointsz_filename); exit(1); }\n",
    "    fread(&num_dest_grid_pointsz, sizeof(int), 1, pointsz_file);\n",
    "\n",
    "    // Step 1.a.i: Sanity check: make sure that num_dest_grid_pointsx == num_dest_grid_pointsy == num_dest_grid_pointsz\n",
    "    if(num_dest_grid_pointsx != num_dest_grid_pointsy || num_dest_grid_pointsy != num_dest_grid_pointsz) {\n",
    "      printf(\"Error: Failed sanity check. Number of interpolation points different in %s-{x,y,z}.dat data files!\\n\",\n",
    "             filename_basename);\n",
    "      exit(1);\n",
    "    } else {\n",
    "      // If sanity check passes:\n",
    "      num_dest_grid_points = num_dest_grid_pointsx;\n",
    "    } // END sanity check\n",
    "\n",
    "      // Step 1.b: Allocate memory for destination grids and interpolation output\n",
    "    if(num_dest_grid_points <= 0 || num_dest_grid_points > 2000000000) {\n",
    "      printf(\"Error: Failed sanity check. Number of interpolation points was found to be: %d\",num_dest_grid_points);\n",
    "      exit(1);\n",
    "    } // END sanity check\n",
    "    points_dest_grid_x = (CCTK_REAL  *)malloc(sizeof(CCTK_REAL)*num_dest_grid_points);\n",
    "    points_dest_grid_y = (CCTK_REAL  *)malloc(sizeof(CCTK_REAL)*num_dest_grid_points);\n",
    "    points_dest_grid_z = (CCTK_REAL  *)malloc(sizeof(CCTK_REAL)*num_dest_grid_points);\n",
    "    output_f = (CCTK_REAL **)malloc(1 * sizeof(CCTK_REAL *));\n",
    "    for(int cc = 0; cc < 1; cc++) output_f[cc]=(CCTK_REAL *)malloc(num_dest_grid_points * sizeof(CCTK_REAL));\n",
    "\n",
    "    // Step 1.c: Store cell-centered points to allocated memory.\n",
    "    fread(points_dest_grid_x, sizeof(CCTK_REAL), num_dest_grid_points, pointsx_file);\n",
    "    fread(points_dest_grid_y, sizeof(CCTK_REAL), num_dest_grid_points, pointsy_file);\n",
    "    fread(points_dest_grid_z, sizeof(CCTK_REAL), num_dest_grid_points, pointsz_file);\n",
    "    int magic_numberx; fread(&magic_numberx, sizeof(int), 1, pointsx_file);\n",
    "    int magic_numbery; fread(&magic_numbery, sizeof(int), 1, pointsy_file);\n",
    "    int magic_numberz; fread(&magic_numberz, sizeof(int), 1, pointsz_file);\n",
    "    int correct_magicnum = -349289480;\n",
    "    if(magic_numberx != correct_magicnum || magic_numbery != correct_magicnum || magic_numberz != correct_magicnum) {\n",
    "      printf(\"Error: Failed sanity check. Magic numbers in x,y,z data files were: %d %d %d, respectively, but should have been: %d\",\n",
    "             magic_numberx,magic_numbery,magic_numberz,correct_magicnum);\n",
    "      exit(1);\n",
    "    }\n",
    "    fclose(pointsx_file);\n",
    "    fclose(pointsy_file);\n",
    "    fclose(pointsz_file);\n",
    "    time(&end_timer); time_in_seconds = difftime(end_timer,start_timer); time(&start_timer);\n",
    "    printf(\"Finished in %e seconds.\\n\",time_in_seconds);\n",
    "\n",
    "    // Step 1.d: Apply offset to x,y,z coordinates to ensure they are centered on (x_center,y_center,z_center)\n",
    "    //           For example, if a black hole is situated at (x,y,z) = (1,2,3), then we set\n",
    "    //           (x_center,y_center,z_center) = (1,2,3) in our ETK parameter file (i.e., with extension .par)\n",
    "    //           and if we desire a point at (x_dest,y_dest,z_dest) = (0,0,0) on the *destination* grid,\n",
    "    //           this will correspond to point (x_src,y_src,z_src) = (1,2,3) = (x_center,y_center,z_center)\n",
    "    //           on the source grid. Thus the translation between source and destination grids is given by\n",
    "    //           (x_src,y_src,z_src) = (x_dest+x_center, y_dest+y_center, z_dest+z_center),\n",
    "    //           where (x_src,y_src,z_src) = (points_dest_grid_x[i],points_dest_grid_y[i],points_dest_grid_z[i]) for point i.\n",
    "    for(int point=0;point<num_dest_grid_points;point++) {\n",
    "      points_dest_grid_x[point] += x_center;\n",
    "      points_dest_grid_y[point] += y_center;\n",
    "      points_dest_grid_z[point] += z_center;\n",
    "    }\n",
    "\n",
    "    // Step 1.e: Look for the points that are out of bounds and set their coordinates to out_of_bounds_interp_xyz.\n",
    "    //           At the end, we will replace the interpolation output at these points with NaNs.\n",
    "    printf(\"Looking for points that are out of bounds.\\n\");\n",
    "    int num_out_of_bounds_points = 0;\n",
    "    for(int i=0;i<num_dest_grid_points;i++){\n",
    "      if(fabs(points_dest_grid_x[i])>out_of_bounds_interp_xyz ||\n",
    "         fabs(points_dest_grid_y[i])>out_of_bounds_interp_xyz ||\n",
    "         fabs(points_dest_grid_z[i])>out_of_bounds_interp_xyz) {\n",
    "        points_dest_grid_x[i] = out_of_bounds_interp_xyz;\n",
    "        points_dest_grid_y[i] = out_of_bounds_interp_xyz;\n",
    "        points_dest_grid_z[i] = out_of_bounds_interp_xyz;\n",
    "        num_out_of_bounds_points++;\n",
    "      }\n",
    "    }\n",
    "    time(&end_timer); time_in_seconds = difftime(end_timer,start_timer); time(&start_timer);\n",
    "    printf(\"Finished in %e seconds, found %i points out of bounds.\\n\", time_in_seconds, num_out_of_bounds_points);\n",
    "\n",
    "  } // END if(CCTK_MyProc(cctkGH)==0)\n",
    "\n",
    "\n",
    "  // Step 1.f: Looping over interp order as desired, interpolate to destination points & output to file\n",
    "  for(int order_i=0; order_i<num_interp_orders; order_i++) {\n",
    "    int order = interp_orders_list[order_i];\n",
    "    printf(\"Interpolating\\033[1m %s \\033[0m... using interpolation order = %d\\n\",gf_name,order);\n",
    "    printf(\"and %s \\n\",interpolator_name);\n",
    "    if(CCTK_MyProc(cctkGH)==0) {\n",
    "      Interpolate_to_dest_grid(cctkGH, num_dest_grid_points, order, interpolator_name,\n",
    "                               points_dest_grid_x,points_dest_grid_y,points_dest_grid_z, input_array_names, output_f);\n",
    "\n",
    "      // Step 1.f.i: Sanity check -- check for bad point:\n",
    "#pragma omp parallel for\n",
    "      for(int i=0;i<num_dest_grid_points;i++) {\n",
    "        if(output_f[0][i] > 1e20) {\n",
    "          printf(\"BAD POINT: %s %d %e %e %e %e\\n\",gf_name,i,points_dest_grid_x[i],points_dest_grid_y[i],points_dest_grid_z[i], output_f[0][i]);\n",
    "          exit(1);\n",
    "        }\n",
    "      }\n",
    "      time(&end_timer); time_in_seconds = difftime(end_timer,start_timer); time(&start_timer);\n",
    "      printf(\"Finished in %d seconds. Next: Filling cells out of bounds with NaNs\\n\",time_in_seconds);\n",
    "      // Step 1.f.ii: Filling cells out of bounds with NaNs:\n",
    "      for(int i=0;i<num_dest_grid_points;i++){\n",
    "        if(fabs(points_dest_grid_x[i])==out_of_bounds_interp_xyz ||\n",
    "           fabs(points_dest_grid_y[i])==out_of_bounds_interp_xyz ||\n",
    "           fabs(points_dest_grid_z[i])==out_of_bounds_interp_xyz) {\n",
    "          output_f[0][i] = 0.0 / 0.0; // ie NaN\n",
    "        }\n",
    "      }\n",
    "      time(&end_timer); time_in_seconds = difftime(end_timer,start_timer); time(&start_timer);\n",
    "      printf(\"Finished in %e seconds. Next: Interpolate_to_dest_grid_main_function(): Outputting to file at iteration %d\\n\",time_in_seconds,cctk_iteration);\n",
    "      output_to_file(CCTK_PASS_CTOC,gf_name,&order,&num_dest_grid_points,output_f);\n",
    "      time(&end_timer); time_in_seconds = difftime(end_timer,start_timer); time(&start_timer);\n",
    "      printf(\"Finished in %e seconds. Interpolate_to_dest_grid_main_function(): Finished output to file at iteration %d\\n\",time_in_seconds,cctk_iteration);\n",
    "    } else {\n",
    "      // On all MPI processes that are nonzero, only call the interpolation function\n",
    "      //    to ensure the MPI calls from the actual interpolation (driven by proc==0) are seen.\n",
    "      // Setting num_dest_grid_points to zero results in a segfault on certain (ahem, Frontera)\n",
    "      //    systems. So we set num_dest_grid_points = 1 and interpolate to the origin\n",
    "      //    only for MPI processes that are nonzero, leaving the heavy lifting to MPI process 0.\n",
    "      num_dest_grid_points = 1;\n",
    "      points_dest_grid_x = (CCTK_REAL  *)malloc(sizeof(CCTK_REAL)*num_dest_grid_points);\n",
    "      points_dest_grid_y = (CCTK_REAL  *)malloc(sizeof(CCTK_REAL)*num_dest_grid_points);\n",
    "      points_dest_grid_z = (CCTK_REAL  *)malloc(sizeof(CCTK_REAL)*num_dest_grid_points);\n",
    "      output_f = (CCTK_REAL **)malloc(1 * sizeof(CCTK_REAL *));\n",
    "      for(int cc = 0; cc < 1; cc++) output_f[cc]=(CCTK_REAL *)malloc(num_dest_grid_points * sizeof(CCTK_REAL));\n",
    "\n",
    "      Interpolate_to_dest_grid(cctkGH, num_dest_grid_points, order, interpolator_name,\n",
    "                               points_dest_grid_x,points_dest_grid_y,points_dest_grid_z, input_array_names, output_f);\n",
    "      output_f[0][0] = 0.0;\n",
    "    } // END if(CCTK_MyProc(cctkGH)==0)\n",
    "  } // END for(int order_i=0; order_i<num_interp_orders; order_i++)\n",
    "\n",
    "  // Step 1.g: Free memory for destination grids and interpolation output\n",
    "  free(points_dest_grid_x);\n",
    "  free(points_dest_grid_y);\n",
    "  free(points_dest_grid_z);\n",
    "  FREE_2D_GENERIC(CCTK_REAL,output_f,1,num_dest_grid_points);\n",
    "} // END function\n",
    "#undef ALLOCATE_2D_GENERIC\n",
    "#undef FREE_2D_GENERIC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='fileformat'></a>\n",
    "\n",
    "## Step 1.b: Outputting to File (File format notes) \\[Back to [top](#toc)\\]\n",
    "$$\\label{fileformat}$$\n",
    "\n",
    "Since they take almost no space relative to the data chunks, we attach the entire metadata to each interpolated function that is output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:10.889506Z",
     "iopub.status.busy": "2021-03-07T17:07:10.888493Z",
     "iopub.status.idle": "2021-03-07T17:07:10.892672Z",
     "shell.execute_reply": "2021-03-07T17:07:10.891873Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/src/output_to_file.h\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/output_to_file.h\n",
    "\n",
    "#include \"define_NumInterpFunctions.h\"\n",
    "\n",
    "// output_to_file() starts order and InterpCounter both with the value 1\n",
    "void output_to_file(CCTK_ARGUMENTS,\n",
    "                    char gf_name[100],\n",
    "                    CCTK_INT *restrict order,\n",
    "                    CCTK_INT *restrict num_interp_points,\n",
    "                    CCTK_REAL *restrict output_f[1]) {\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  char filename[100];\n",
    "  sprintf (filename, \"%s/interp_dest_grids_MO.dat\", out_dir);\n",
    "  FILE *file;\n",
    "  if(*InterpCounter == 1 && *order==1) {\n",
    "    file = fopen (filename,\"w\");\n",
    "    printf(\"WRITING to file %s\\n\",filename);\n",
    "    // Write EOS information to the beginning of the file\n",
    "    char eos_info[256];\n",
    "    if( enable_nuc_eos ) {\n",
    "      // If using nuclear EOS, give the table name\n",
    "      sprintf(eos_info,\"Nuclear (tabulated) EOS from file %s\",nuceos_table_name);\n",
    "    }\n",
    "    else {\n",
    "      // If using hybrid EOS, then give Gamma_th and the number of polytropic pieces\n",
    "      sprintf(eos_info,\"Hybrid EOS with Gamma_th = %.15e and %d polytropic pieces\",Gamma_th,neos);\n",
    "    }\n",
    "    fwrite(eos_info, 256*sizeof(char), 1, file);\n",
    "  }\n",
    "  else {\n",
    "    file = fopen (filename,\"a+\");\n",
    "    printf(\"Appending to file %s\\n\",filename);\n",
    "  }\n",
    "  if (! file) {\n",
    "    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,\n",
    "                \"interp_dest_grid__ET_thorn: Cannot open output file '%s'\", filename);\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  fwrite(gf_name, 100*sizeof(char), 1, file);\n",
    "  fwrite(order, sizeof(CCTK_INT), 1, file);\n",
    "  fwrite(num_interp_points, sizeof(CCTK_INT),1,file);\n",
    "\n",
    "  CCTK_REAL magic_number = 1.130814081305130e-21;\n",
    "  fwrite(&magic_number, sizeof(CCTK_REAL), 1, file);\n",
    "  fwrite(&cctk_iteration, sizeof(CCTK_INT), 1, file);\n",
    "  fwrite(&cctk_time, sizeof(CCTK_REAL), 1, file);\n",
    "  for(CCTK_INT i=0;i<1;i++) {\n",
    "    fwrite(output_f[i], sizeof(CCTK_REAL)*(*num_interp_points), 1, file);\n",
    "  }\n",
    "\n",
    "  fclose(file);\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='maininterpolator'></a>\n",
    "\n",
    "## Step 1.c: The Main Interpolation Driver Function \\[Back to [top](#toc)\\]\n",
    "$$\\label{maininterpolator}$$\n",
    "\n",
    "The **`Interpolate_to_dest_grid_main_function()`** function calls the above functions as follows:\n",
    "\n",
    "1. **`Interpolate_to_dest_grid()`** ([Above section](#etk_interp)): Interpolates to destination grid and calls\n",
    "    1. **`output_to_file()`** ([Above section](#fileformat)): Outputs information about interpolation, as well as interpolation result, to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:10.900193Z",
     "iopub.status.busy": "2021-03-07T17:07:10.899051Z",
     "iopub.status.idle": "2021-03-07T17:07:10.903650Z",
     "shell.execute_reply": "2021-03-07T17:07:10.902975Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/src/main_function.cc\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/main_function.cc\n",
    "\n",
    "\n",
    "// Include needed ETK & C library header files:\n",
    "#include <stdio.h>\n",
    "#include <stdlib.h>\n",
    "#include <string.h>\n",
    "#include <math.h>\n",
    "#include <time.h> // for benchmarking\n",
    "// Needed for dealing with Cactus/ETK infrastructure\n",
    "#include \"cctk.h\"\n",
    "#include \"cctk_Arguments.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "// Needed for low-level interpolation functions\n",
    "#include \"util_Table.h\"\n",
    "#include \"util_String.h\"\n",
    "\n",
    "// Include locally-defined C++ functions:\n",
    "#include \"Interpolate_to_dest_grid.h\"\n",
    "#include \"get_gf_name.h\"\n",
    "#include \"interpolate_set_of_points_in_file.h\"\n",
    "\n",
    "void Interpolate_to_dest_grid_main_function(CCTK_ARGUMENTS) {\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  // Perform interpolation only at iteration == interp_out_iteration:\n",
    "  if(cctk_iteration != interp_out_iteration) return;\n",
    "  // Perform interpolation!\n",
    "  // Process zero (CCTK_MyProc(cctkGH)==0) is responsible for directing the interpolation.\n",
    "  //    All other processes must see the cctk_InterpGridArrays() within Interpolate_to_dest_grid(),\n",
    "  //    so that the MPI calls work properly, but these nonzero processes can call\n",
    "  //    Interpolate_to_dest_grid() with number of interpolated points set to zero, and\n",
    "  //    without needing a malloc().\n",
    "  char gf_name[100]; get_gf_name(*InterpCounter,gf_name);\n",
    "  char filename_basename[100];\n",
    "  char interpolator_name[100];\n",
    "\n",
    "  // The following \"if\" statement is destination-code dependent.\n",
    "  // In our case, since we are interpolating variables for harm3d,\n",
    "  // we interpolate the vector potential to the corners of each cell.\n",
    "  // Every other quantity is interpolated to the center of each cell.\n",
    " if(strncmp(gf_name,\"Unstaggered\",11) == 0){\n",
    "    sprintf(filename_basename,\"corner_points\");\n",
    "  }\n",
    "  else{\n",
    "    sprintf(filename_basename,\"cell_centered_points\");\n",
    "  }\n",
    "\n",
    "  int num_interp_orders,*interp_orders_list;\n",
    "  // 4-metric, 4-Christoffels and A_mu only output Hermite interpolation order==3\n",
    "  // so we secure continuity of first derivative.\n",
    "  if( (strncmp(gf_name,\"4-\",2) == 0) ||\n",
    "      (strncmp(gf_name,\"Unstaggered\",11) == 0) ||\n",
    "     (strncmp(gf_name,\"Constraints\",11) == 0)) {\n",
    "    num_interp_orders = 1;\n",
    "    sprintf(interpolator_name, \"Hermite polynomial interpolation\");\n",
    "    interp_orders_list = (int *)malloc(sizeof(int)*num_interp_orders);\n",
    "    interp_orders_list[0] = 3;\n",
    "  } else {\n",
    "    num_interp_orders = 3;\n",
    "    sprintf(interpolator_name, \"Lagrange polynomial interpolation\");\n",
    "    interp_orders_list = (int *)malloc(sizeof(int)*num_interp_orders);\n",
    "    int count = 0; for(int order=1;order<=4;order*=2) { interp_orders_list[count] = order; count++; }\n",
    "  }\n",
    "\n",
    "  interpolate_set_of_points_in_file(CCTK_PASS_CTOC,filename_basename,gf_name,interpolator_name,num_interp_orders,interp_orders_list);\n",
    "  free(interp_orders_list);\n",
    "\n",
    "  // Now perform interpolation of 4-metric on\n",
    "  //   faces (i-1/2,j,k), (i,j-1/2,k), (i,j,k-1/2) and corners (i-1/2,j-1/2,k-1/2)\n",
    "  if(strncmp(gf_name,\"4-metric\",8) == 0) {\n",
    "    num_interp_orders = 1;\n",
    "    sprintf(interpolator_name, \"Hermite polynomial interpolation\");\n",
    "    interp_orders_list = (int *)malloc(sizeof(int)*num_interp_orders);\n",
    "    interp_orders_list[0] = 3;\n",
    "\n",
    "    char gf_name_new[100];\n",
    "\n",
    "    sprintf(filename_basename,\"faceim_points\");\n",
    "    snprintf(gf_name_new,100,\"faceim (i-1/2,j,k): %s\",gf_name);\n",
    "    interpolate_set_of_points_in_file(CCTK_PASS_CTOC,filename_basename,gf_name_new,interpolator_name,num_interp_orders,interp_orders_list);\n",
    "\n",
    "    sprintf(filename_basename,\"facejm_points\");\n",
    "    snprintf(gf_name_new,100,\"facejm (i,j-1/2,k): %s\",gf_name);\n",
    "    interpolate_set_of_points_in_file(CCTK_PASS_CTOC,filename_basename,gf_name_new,interpolator_name,num_interp_orders,interp_orders_list);\n",
    "\n",
    "    sprintf(filename_basename,\"facekm_points\");\n",
    "    snprintf(gf_name_new,100,\"facekm (i,j,k-1/2): %s\",gf_name);\n",
    "    interpolate_set_of_points_in_file(CCTK_PASS_CTOC,filename_basename,gf_name_new,interpolator_name,num_interp_orders,interp_orders_list);\n",
    "\n",
    "    sprintf(filename_basename,\"corner_points\");\n",
    "    snprintf(gf_name_new,100,\"cornr (i-1/2,j-1/2,k-1/2): %s\",gf_name);\n",
    "    interpolate_set_of_points_in_file(CCTK_PASS_CTOC,filename_basename,gf_name_new,interpolator_name,num_interp_orders,interp_orders_list);\n",
    "  } // END if(strncmp(gf_name,\"4-metric\",8) == 0)\n",
    "} // END function\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='standalonerandompoints'></a>\n",
    "\n",
    "## Step 1.d: Standalone C code to output random points data \\[Back to [top](#toc)\\]\n",
    "$$\\label{standalonerandompoints}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:10.910393Z",
     "iopub.status.busy": "2021-03-07T17:07:10.909429Z",
     "iopub.status.idle": "2021-03-07T17:07:10.913131Z",
     "shell.execute_reply": "2021-03-07T17:07:10.913594Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/src/standalone/standalone_C_code_genpoints.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/standalone/standalone_C_code_genpoints.c\n",
    "\n",
    "// Part P1: Import needed header files\n",
    "#include \"stdio.h\"\n",
    "#include \"stdlib.h\"\n",
    "#include \"math.h\"\n",
    "\n",
    "const double xyzmin = -1000.0;\n",
    "const double xyzmax =  1000.0;\n",
    "\n",
    "void write_to_xyz_files(int num_interp_points, char filename_basename[100]) {\n",
    "  char filenamex[100],filenamey[100],filenamez[100];\n",
    "  snprintf(filenamex,100,\"%s-x.dat\",filename_basename);\n",
    "  snprintf(filenamey,100,\"%s-y.dat\",filename_basename);\n",
    "  snprintf(filenamez,100,\"%s-z.dat\",filename_basename);\n",
    "  FILE *filex = fopen(filenamex,\"wb\");\n",
    "  FILE *filey = fopen(filenamey,\"wb\");\n",
    "  FILE *filez = fopen(filenamez,\"wb\");\n",
    "\n",
    "  // Write file headers:\n",
    "  fwrite(&num_interp_points, sizeof(int), 1, filex);\n",
    "  fwrite(&num_interp_points, sizeof(int), 1, filey);\n",
    "  fwrite(&num_interp_points, sizeof(int), 1, filez);\n",
    "\n",
    "  // Write guts of file:\n",
    "  for(int ii=0;ii<num_interp_points;ii++) {\n",
    "    double rngx = xyzmin + (xyzmax - xyzmin)*drand48(); // drand48() returns between 0.0 & 1.0\n",
    "    double rngy = xyzmin + (xyzmax - xyzmin)*drand48();\n",
    "    double rngz = xyzmin + (xyzmax - xyzmin)*drand48();\n",
    "    fwrite(&rngx, sizeof(double), 1, filex);\n",
    "    fwrite(&rngy, sizeof(double), 1, filey);\n",
    "    fwrite(&rngz, sizeof(double), 1, filez);\n",
    "  }\n",
    "\n",
    "  // Write magic number as file footers:\n",
    "  int magic_number = -349289480;\n",
    "  fwrite(&magic_number, sizeof(int), 1, filex);\n",
    "  fwrite(&magic_number, sizeof(int), 1, filey);\n",
    "  fwrite(&magic_number, sizeof(int), 1, filez);\n",
    "\n",
    "  // Close files.\n",
    "  fclose(filex);\n",
    "  fclose(filey);\n",
    "  fclose(filez);\n",
    "}\n",
    "\n",
    "int main(int argc, const char *argv[]) {\n",
    "\n",
    "  // Step 0a: Read command-line input, error out if nonconformant\n",
    "  if(argc != 2 || atoi(argv[1]) < 1) {\n",
    "    printf(\"Error: Expected one command-line argument: ./standalone_C_code_genpoints [num_interp_points],\\n\");\n",
    "    exit(1);\n",
    "  }\n",
    "\n",
    "  const int num_interp_points = atoi(argv[1]);\n",
    "\n",
    "  char filename_basename[100];\n",
    "  sprintf(filename_basename,\"cell_centered_points\");\n",
    "  write_to_xyz_files(num_interp_points, filename_basename);\n",
    "  sprintf(filename_basename,\"faceim_points\");\n",
    "  write_to_xyz_files(num_interp_points, filename_basename);\n",
    "  sprintf(filename_basename,\"facejm_points\");\n",
    "  write_to_xyz_files(num_interp_points, filename_basename);\n",
    "  sprintf(filename_basename,\"facekm_points\");\n",
    "  write_to_xyz_files(num_interp_points, filename_basename);\n",
    "  sprintf(filename_basename,\"corner_points\");\n",
    "  write_to_xyz_files(num_interp_points, filename_basename);\n",
    "\n",
    "  return 0;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy'></a>\n",
    "\n",
    "# Step 2: Use NRPy+ C Output to Set All Output Gridfunctions \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:10.925126Z",
     "iopub.status.busy": "2021-03-07T17:07:10.924018Z",
     "iopub.status.idle": "2021-03-07T17:07:11.247771Z",
     "shell.execute_reply": "2021-03-07T17:07:11.248295Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Import needed NRPy+ parameters\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import grid as gri               # NRPy+: Functions having to do with numerical grids\n",
    "import finite_difference as fin  # NRPy+: Finite difference C code generation module\n",
    "from outputC import lhrh         # NRPy+: Core C code output module\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import loop as lp                # NRPy+: Generate C code loops\n",
    "\n",
    "par.set_parval_from_str(\"grid::GridFuncMemAccess\",\"ETK\")\n",
    "\n",
    "from collections import namedtuple\n",
    "gf_interp = namedtuple('gf_interp', 'gf_description')\n",
    "gf_interp_list = []\n",
    "gf_interp_list.append(gf_interp(\"dummy -- used because this is a 1-offset array\"))\n",
    "\n",
    "interped_gf = gri.register_gridfunctions(\"AUX\",\"interped_gf\")\n",
    "\n",
    "def interp_fileout(which_InterpCounter, expression, filename):\n",
    "    kernel = fin.FD_outputC(\"returnstring\",lhrh(lhs=gri.gfaccess(\"out_gfs\",\"interped_gf\"),rhs=expression),\"outCverbose=False\")\n",
    "    output_type=\"a\"\n",
    "    if which_InterpCounter == 1:\n",
    "        # Write the file header, which includes #define's for GAMMA_SPEED_LIMIT and TINYDOUBLE:\n",
    "        with open(filename, \"w\") as file:\n",
    "            file.write(\"// Parameters needed to ensure velocity computations are robust:\\n\")\n",
    "            file.write(\"#define GAMMA_SPEED_LIMIT 20\\n\")\n",
    "            file.write(\"#define TINYDOUBLE        1e-100\\n\\n\")\n",
    "    compute_xx0xx1xx2 = \"\"\n",
    "    if \"SPHERICAL\" in gf_interp_list[which_InterpCounter].gf_description:\n",
    "        compute_xx0xx1xx2 = \"\"\"\n",
    "// ONLY NEEDED/USED IF CONVERTING TO SPHERICAL BASIS:\n",
    "const double Cartx = x[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)] - x_center;\n",
    "const double Carty = y[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)] - y_center;\n",
    "const double Cartz = z[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)] - z_center;\n",
    "\n",
    "const double xx0 = sqrt(Cartx*Cartx + Carty*Carty + Cartz*Cartz);\n",
    "const double xx1 = acos(Cartz/xx0);\n",
    "const double xx2 = atan2(Carty,Cartx);\\n\n",
    "\"\"\"\n",
    "    with open(filename, output_type) as file:\n",
    "        if( \"temperature\" in str(interp_expr) or\n",
    "            \"Y_e\" in str(interp_expr) or\n",
    "            \"eps\" in str(interp_expr) or\n",
    "            \"entropy\"in str(interp_expr) ):\n",
    "            file.write(\"if(*InterpCounter == \"+str(which_InterpCounter)+\" && enable_nuc_eos) {\\n\")\n",
    "        else:\n",
    "            file.write(\"if(*InterpCounter == \"+str(which_InterpCounter)+\") {\\n\")\n",
    "        file.write(\"    // Interpolating: \"+gf_interp_list[which_InterpCounter].gf_description+\"\\n\")\n",
    "        file.write(lp.loop([\"i2\",\"i1\",\"i0\"],\n",
    "                           [\"cctk_nghostzones[2]\",\"cctk_nghostzones[1]\",\"cctk_nghostzones[0]\"],\\\n",
    "                           [\"cctk_lsh[2]-cctk_nghostzones[2]\",\n",
    "                            \"cctk_lsh[1]-cctk_nghostzones[1]\",\n",
    "                            \"cctk_lsh[0]-cctk_nghostzones[0]\"],\\\n",
    "                           [\"1\",\"1\",\"1\"],\\\n",
    "                           [\"#pragma omp parallel for\",\"\",\"\"],\"   \",\n",
    "                           compute_xx0xx1xx2+kernel))\n",
    "        file.write(\"}\\n\")\n",
    "    # If successful, return incremented which_InterpCounter:\n",
    "    return which_InterpCounter+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy_list_of_funcs_interp'></a>\n",
    "\n",
    "## Step 2.a: Set up NRPy-based `list_of_functions_to_interpolate.h` \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy_list_of_funcs_interp}$$\n",
    "\n",
    "First specify NRPy+ output file and initialize `which_InterpCounter`, which keeps track of the number of interpolated functions on the grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:11.252796Z",
     "iopub.status.busy": "2021-03-07T17:07:11.252142Z",
     "iopub.status.idle": "2021-03-07T17:07:11.254402Z",
     "shell.execute_reply": "2021-03-07T17:07:11.254935Z"
    }
   },
   "outputs": [],
   "source": [
    "NRPyoutfilename = os.path.join(Ccodesdir,\"src\",\"list_of_functions_to_interpolate.h\")\n",
    "\n",
    "which_InterpCounter = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpygrmhd'></a>\n",
    "\n",
    "### Step 2.a.i: GRMHD quantities \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpygrmhd}$$\n",
    "\n",
    "These include\n",
    "* $\\rho_b$, the baryonic density (e.g. $\\text{HydroBase::rho}$ or $\\text{IllinoisGRMHD::rho_b}$)\n",
    "* $P$, the total gas pressure (e.g. $\\text{HydroBase::press}$ or $\\text{IllinoisGRMHD::P}$)\n",
    "* $\\Gamma v_{(n)}^i$, the Valencia 3-velocity times the Lorentz factor (i.e., the HydroBase 3-gridfuntion $\\verb|vel|$, multiplied by the Lorentz factor). This definition of velocity has the advantage that after interpolation, it will not violate $u^\\mu u_\\mu = -1$. In terms of the IllinoisGRMHD 3-velocity $v^i = u^i / u^0$, the Valencia 3-velocity is given by (Eq. 11 of [Etienne *et al*](https://arxiv.org/pdf/1501.07276.pdf)):\n",
    "$$\n",
    "v_{(n)}^i = \\frac{1}{\\alpha} \\left(v^i + \\beta^i\\right).\n",
    "$$\n",
    "Further, $\\Gamma = \\alpha u^0$ is given by (as shown [here](Tutorial-u0_smallb_Poynting-Cartesian.ipynb)):\n",
    "$$\n",
    "\\Gamma = \\alpha u^0 = \\sqrt{\\frac{1}{1 - \\gamma_{ij}v^i_{(n)}v^j_{(n)}}}.\n",
    "$$\n",
    "Therefore, $\\Gamma v_{(n)}^i$ is given by\n",
    "$$\n",
    "\\Gamma v_{(n)}^i = \\frac{1}{\\alpha} \\left(v^i + \\beta^i\\right) \\sqrt{\\frac{1}{1 - \\gamma_{ij}v^i_{(n)}v^j_{(n)}}}.\n",
    "$$\n",
    "* $A_i$, the *unstaggered* magnetic vector potential.\n",
    "* $B^i$, the *unstaggered* magnetic field vector (output only for validation purposes)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:11.265836Z",
     "iopub.status.busy": "2021-03-07T17:07:11.265122Z",
     "iopub.status.idle": "2021-03-07T17:07:11.356611Z",
     "shell.execute_reply": "2021-03-07T17:07:11.355986Z"
    }
   },
   "outputs": [],
   "source": [
    "# INPUT GRIDFUNCTIONS: The AUX or EVOL designation is *not* used in diagnostic modules.\n",
    "\n",
    "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUX\",\"gammaDD\", \"sym01\")\n",
    "alpha = gri.register_gridfunctions(\"AUX\",\"alpha\")\n",
    "betaU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"betaU\")\n",
    "\n",
    "# Add a constant beta offset, to account for linear\n",
    "#      (i.e., constant velocity) coordinate drift.\n",
    "# Note that beta_offsetU's are set in param.ccl.\n",
    "# As beta_offsetU is constant in space, it has no\n",
    "#      impact on betaU_dD's.\n",
    "beta_offsetU0,beta_offsetU1,beta_offsetU2 = par.Cparameters(\"REAL\",\"modulenamedoesntmatter\",\n",
    "                                                            [\"beta_offsetU0\",\"beta_offsetU1\",\"beta_offsetU2\"],\n",
    "                                                            [0.0,0.0,0.0])\n",
    "betaU[0] += beta_offsetU0\n",
    "betaU[1] += beta_offsetU1\n",
    "betaU[2] += beta_offsetU2\n",
    "\n",
    "# Tensors are given in Cartesian basis:\n",
    "# Derivatives of metric\n",
    "gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\")\n",
    "betaU_dD   = ixp.declarerank2(\"betaU_dD\",\"nosym\")\n",
    "alpha_dD   = ixp.declarerank1(\"alpha_dD\")\n",
    "\n",
    "DIM=3\n",
    "\n",
    "IGMvU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"IGMvU\")\n",
    "BU    = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"BU\")\n",
    "\n",
    "gf_interp_list.append(gf_interp(\"IGM density primitive\"))\n",
    "rho_b       = gri.register_gridfunctions(\"AUX\",\"rho_b\")\n",
    "interp_expr = rho_b\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "gf_interp_list.append(gf_interp(\"IGM pressure primitive\"))\n",
    "P = gri.register_gridfunctions(\"AUX\",\"P\")\n",
    "interp_expr = P\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we implement:\n",
    "$$\n",
    "v_{(n)}^i = \\frac{1}{\\alpha} \\left(v^i + \\beta^i\\right),\n",
    "$$\n",
    "and\n",
    "$$\n",
    "\\Gamma v_{(n)}^i = \\sqrt{\\frac{1}{1 - \\gamma_{ij}v^i_{(n)}v^j_{(n)}}} v_{(n)}^i.\n",
    "$$\n",
    "\n",
    "We use expressions for $v_{(n)}^i$, $\\Gamma v_{(n)}^i$, and $v^i$ as implemented [in the GRHD equations notebook](Tutorial-GRHD_Equations-Cartesian.ipynb#convertvtou) and corresponding [GRHD.equations Python module](../edit/GRHD/equations.py). These expressions enforce a speed limit on $\\Gamma$ to ensure e.g., the denominator within the radical in the above expression for $\\Gamma v_{(n)}^i$ is never negative, which would result in `NaN`s. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:11.364586Z",
     "iopub.status.busy": "2021-03-07T17:07:11.363578Z",
     "iopub.status.idle": "2021-03-07T17:07:11.454362Z",
     "shell.execute_reply": "2021-03-07T17:07:11.453784Z"
    }
   },
   "outputs": [],
   "source": [
    "import GRHD.equations as Ge\n",
    "Ge.u4U_in_terms_of_vU__rescale_vU_by_applying_speed_limit(alpha, betaU, gammaDD, IGMvU)\n",
    "\n",
    "# ValenciavU = ixp.zerorank1()\n",
    "# for i in range(DIM):\n",
    "#     ValenciavU[i] = 1/alpha * (IGMvU[i] + betaU[i])\n",
    "# v_dot_v = sp.sympify(0)\n",
    "# for i in range(DIM):\n",
    "#     for j in range(DIM):\n",
    "#         v_dot_v += gammaDD[i][j]*ValenciavU[i]*ValenciavU[j]\n",
    "# u4Uzero = sp.sqrt(1/(1 - v_dot_v))/alpha # u^0 = LorentzGamma/alpha\n",
    "\n",
    "u4Uzero = Ge.u4U_ito_vU[0]\n",
    "gf_interp_list.append(gf_interp(\"u^0: zero (time) component of 4-velocity\"))\n",
    "interp_expr = u4Uzero\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "Gamma_times_ValenciavU = ixp.zerorank1()\n",
    "Gamma      = alpha*Ge.u4U_ito_vU[0]\n",
    "ValenciavU = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    ValenciavU[i] = (1/alpha * (Ge.rescaledvU[i] + betaU[i])) # simplify?\n",
    "\n",
    "for i in range(DIM):\n",
    "    Gamma_times_ValenciavU[i] = Gamma*ValenciavU[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we'll basis transform $W^i = \\Gamma v_{\\rm n}^i$ from Cartesian to the spherical basis:\n",
    "\n",
    "Within [`reference_metric.py`](../edit/reference_metric.py), the `compute_Jacobian_and_inverseJacobian_tofrom_Cartesian()` function defines Jacobians relative to the center of the source (reference metric) grid, at a point $x^j_{\\rm src}=$(`xx0,xx1,xx2`)${}_{\\rm src}$ on the source grid:\n",
    "$$\n",
    "{\\rm Jac\\_dUCart\\_dDsrcUD[i][j]} = \\frac{\\partial x^i_{\\rm Cart}}{\\partial x^j_{\\rm src}},\n",
    "$$\n",
    "\n",
    "via exact differentiation (courtesy SymPy), and the inverse Jacobian\n",
    "$$\n",
    "{\\rm Jac\\_dUsrc\\_dDCartUD[i][j]} = \\frac{\\partial x^i_{\\rm src}}{\\partial x^j_{\\rm Cart}},\n",
    "$$\n",
    "\n",
    "using NRPy+'s `generic_matrix_inverter3x3()` function. \n",
    "\n",
    "In terms of these, the transformation of $W^i$ from Cartesian coordinates to `\"reference_metric::CoordSystem=Spherical\"`may be written:\n",
    "\n",
    "\\begin{align}\n",
    "W^i_{\\rm Sph} &= \\frac{\\partial x^i_{\\rm Sph}}{\\partial x^\\ell_{\\rm Cart}} W^\\ell_{\\rm Cart}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:11.461298Z",
     "iopub.status.busy": "2021-03-07T17:07:11.460552Z",
     "iopub.status.idle": "2021-03-07T17:07:12.240456Z",
     "shell.execute_reply": "2021-03-07T17:07:12.239736Z"
    }
   },
   "outputs": [],
   "source": [
    "import reference_metric as rfm    # NRPy+: Reference metric support\n",
    "\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Spherical\")\n",
    "rfm.reference_metric()\n",
    "\n",
    "# Step 2.a: Construct Jacobian & Inverse Jacobians:\n",
    "Jac_dUCart_dDrfmUD,Jac_dUrfm_dDCartUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Cartesian()\n",
    "\n",
    "Sph_basis_Gamma_times_ValenciavU = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    for l in range(DIM):\n",
    "        Sph_basis_Gamma_times_ValenciavU[i] += Jac_dUrfm_dDCartUD[i][l] * Gamma_times_ValenciavU[l]\n",
    "\n",
    "for i in range(DIM):\n",
    "    gf_interp_list.append(gf_interp(\"*SPHERICAL BASIS* Lorentz factor, times Valencia vU\"+str(i)))\n",
    "    interp_expr = Sph_basis_Gamma_times_ValenciavU[i]\n",
    "    which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:12.299807Z",
     "iopub.status.busy": "2021-03-07T17:07:12.263835Z",
     "iopub.status.idle": "2021-03-07T17:07:12.685637Z",
     "shell.execute_reply": "2021-03-07T17:07:12.685008Z"
    }
   },
   "outputs": [],
   "source": [
    "for i in range(DIM):\n",
    "    gf_interp_list.append(gf_interp(\"(speed-limited) Valencia 3-velocity vU\"+str(i)))\n",
    "    interp_expr = ValenciavU[i]\n",
    "    which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "for i in range(DIM):\n",
    "    if i==0:\n",
    "        gf_interp_list.append(gf_interp(\"Local grid resolution dx=dy=dz\"))\n",
    "        invdx0 = sp.symbols('invdx0', real=True)\n",
    "        interp_expr = 1/invdx0\n",
    "    else:\n",
    "        gf_interp_list.append(gf_interp(\"(speed-limited) IGM 3-velocity vU\"+str(i)+\" = u^i divided by u^0\"))\n",
    "        interp_expr = Ge.u4U_ito_vU[i]\n",
    "    which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "# For testing:\n",
    "# gf_interp_list.append(gf_interp(\"Lorentz factor\"))\n",
    "# interp_expr = v_dot_v\n",
    "# which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "# for i in range(DIM):\n",
    "#     gf_interp_list.append(gf_interp(\"Valencia vU\"+str(i)))\n",
    "#     interp_expr = Valenciav[i]\n",
    "#     which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "for i in range(DIM):\n",
    "    gf_interp_list.append(gf_interp(\"IGM magnetic field component B\"+str(i)))\n",
    "    interp_expr = BU[i]\n",
    "    which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='unstaggera'></a>\n",
    "\n",
    "### Step 2.a.ii: Unstagger $A_i$ and add to \"list of functions to interpolate\" \\[Back to [top](#toc)\\]\n",
    "$$\\label{unstaggera}$$\n",
    "\n",
    "First generate the C code needed to unstagger the A-fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:12.691469Z",
     "iopub.status.busy": "2021-03-07T17:07:12.690631Z",
     "iopub.status.idle": "2021-03-07T17:07:12.693950Z",
     "shell.execute_reply": "2021-03-07T17:07:12.694545Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/src/unstagger_A_fields.cc\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/unstagger_A_fields.cc\n",
    "\n",
    "// Include needed ETK & C library header files:\n",
    "#include <stdio.h>\n",
    "#include <stdlib.h>\n",
    "#include <string.h>\n",
    "#include <math.h>\n",
    "// Needed for dealing with Cactus/ETK infrastructure\n",
    "#include \"cctk.h\"\n",
    "#include \"cctk_Arguments.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "\n",
    "void unstagger_A_fields(CCTK_ARGUMENTS) {\n",
    "    DECLARE_CCTK_ARGUMENTS;\n",
    "    DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "    // Set Ai_unstaggered = Ai and exit the function if A fields are unstaggered already.\n",
    "    if(A_fields_are_staggered == 0) {\n",
    "#pragma omp parallel for\n",
    "        for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {\n",
    "            int index=CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "            Ax_unstaggered[index] = Ax[index];\n",
    "            Ay_unstaggered[index] = Ay[index];\n",
    "            Az_unstaggered[index] = Az[index];\n",
    "        }\n",
    "        return;\n",
    "    }\n",
    "    printf(\"Unstaggering A fields on grid with dx = %e!\\n\",CCTK_DELTA_SPACE(0));\n",
    "    // If A fields are staggered (IllinoisGRMHD-style), then unstagger them:\n",
    "    // First unstagger A_x, which is defined at (i, j+1/2, k+1/2). Unstaggering\n",
    "    //   is as simple as A_x(i,j,k) = 1/4 * (A_x(i,j-1/2,k-1/2)+A_x(i,j-1/2,k+1/2)+A_x(i,j+1/2,k-1/2)+A_x(i,j+1/2,k+1/2))\n",
    "#pragma omp parallel for\n",
    "    for(int k=1;k<cctk_lsh[2];k++) for(int j=1;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {\n",
    "        int index=CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "        Ax_unstaggered[index] = 0.25*(Ax[CCTK_GFINDEX3D(cctkGH,i,j,k)]     + Ax[CCTK_GFINDEX3D(cctkGH,i,j-1,k)] +\n",
    "                                      Ax[CCTK_GFINDEX3D(cctkGH,i,j-1,k-1)] + Ax[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]);\n",
    "    }\n",
    "#pragma omp parallel for\n",
    "    for(int k=1;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=1;i<cctk_lsh[0];i++) {\n",
    "        int index=CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "        Ay_unstaggered[index] = 0.25*(Ay[CCTK_GFINDEX3D(cctkGH,i,j,k)]     + Ay[CCTK_GFINDEX3D(cctkGH,i-1,j,k)] +\n",
    "                                      Ay[CCTK_GFINDEX3D(cctkGH,i-1,j,k-1)] + Ay[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]);\n",
    "    }\n",
    "#pragma omp parallel for\n",
    "    for(int k=0;k<cctk_lsh[2];k++) for(int j=1;j<cctk_lsh[1];j++) for(int i=1;i<cctk_lsh[0];i++) {\n",
    "        int index=CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "        Az_unstaggered[index] = 0.25*(Az[CCTK_GFINDEX3D(cctkGH,i,j,k)]     + Az[CCTK_GFINDEX3D(cctkGH,i-1,j,k)] +\n",
    "                                      Az[CCTK_GFINDEX3D(cctkGH,i-1,j-1,k)] + Az[CCTK_GFINDEX3D(cctkGH,i,j-1,k)]);\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we instruct NRPy+ to interpolate the unstaggered gridfunctions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:12.705658Z",
     "iopub.status.busy": "2021-03-07T17:07:12.704820Z",
     "iopub.status.idle": "2021-03-07T17:07:12.715058Z",
     "shell.execute_reply": "2021-03-07T17:07:12.715588Z"
    }
   },
   "outputs": [],
   "source": [
    "Ax_unstaggered = gri.register_gridfunctions(\"AUX\",\"Ax_unstaggered\")\n",
    "Ay_unstaggered = gri.register_gridfunctions(\"AUX\",\"Ay_unstaggered\")\n",
    "Az_unstaggered = gri.register_gridfunctions(\"AUX\",\"Az_unstaggered\")\n",
    "\n",
    "gf_interp_list.append(gf_interp(\"Unstaggered vector potential component Ax\"))\n",
    "interp_expr = Ax_unstaggered\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "gf_interp_list.append(gf_interp(\"Unstaggered vector potential component Ay\"))\n",
    "interp_expr = Ay_unstaggered\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "gf_interp_list.append(gf_interp(\"Unstaggered vector potential component Az\"))\n",
    "interp_expr = Az_unstaggered\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy4metric'></a>\n",
    "\n",
    "### Step 2.a.iii: Compute all 10 components of the 4-metric $g_{\\mu\\nu}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy4metric}$$\n",
    "\n",
    "We are given $\\gamma_{ij}$, $\\alpha$, and $\\beta^i$ from ADMBase, and the 4-metric is given in terms of these quantities as\n",
    "$$\n",
    "g_{\\mu\\nu} = \\begin{pmatrix} \n",
    "-\\alpha^2 + \\beta^k \\beta_k & \\beta_i \\\\\n",
    "\\beta_j & \\gamma_{ij}\n",
    "\\end{pmatrix}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:12.761512Z",
     "iopub.status.busy": "2021-03-07T17:07:12.759214Z",
     "iopub.status.idle": "2021-03-07T17:07:12.824738Z",
     "shell.execute_reply": "2021-03-07T17:07:12.825360Z"
    }
   },
   "outputs": [],
   "source": [
    "# Eq. 2.121 in B&S\n",
    "betaD = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        betaD[i] += gammaDD[i][j]*betaU[j]\n",
    "\n",
    "# Now compute the beta contraction.\n",
    "beta2 = sp.sympify(0)\n",
    "for i in range(DIM):\n",
    "    beta2 += betaU[i]*betaD[i]\n",
    "\n",
    "# Eq. 2.122 in B&S\n",
    "g4DD = ixp.zerorank2(DIM=4)\n",
    "g4DD[0][0] = -alpha**2 + beta2\n",
    "for i in range(DIM):\n",
    "    g4DD[i+1][0] = g4DD[0][i+1] = betaD[i]\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        g4DD[i+1][j+1] = gammaDD[i][j]\n",
    "\n",
    "for mu in range(4):\n",
    "    for nu in range(mu,4):\n",
    "        gf_interp_list.append(gf_interp(\"4-metric component g4DD\"+str(mu)+str(nu)))\n",
    "        interp_expr = g4DD[mu][nu]\n",
    "        which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy4christoffels_cartesian'></a>\n",
    "\n",
    "### Step 2.a.iv: Compute all 40 4-Christoffels $\\Gamma^{\\mu}_{\\nu\\delta}$ in Cartesian coordinates \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy4christoffels_cartesian}$$\n",
    "\n",
    "$$\n",
    "\\Gamma^{\\mu}_{\\nu\\delta} = \\frac{1}{2} g^{\\mu\\eta} \\left(g_{\\eta\\nu,\\delta} + g_{\\eta\\delta,\\nu} - g_{\\nu\\delta,\\eta}  \\right)\n",
    "$$\n",
    "\n",
    "Recall that $g_{\\mu\\nu}$ is given from $\\gamma_{ij}$, $\\alpha$, and $\\beta^i$ via\n",
    "$$\n",
    "g_{\\mu\\nu} = \\begin{pmatrix} \n",
    "-\\alpha^2 + \\beta^k \\beta_k & \\beta_i \\\\\n",
    "\\beta_j & \\gamma_{ij}\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "The derivatives $g_{\\mu\\nu,\\eta}$ are then computed in terms of finite-difference derivatives of the input ADM gridfunctions $\\gamma_{ij}$, $\\alpha$, and $\\beta^i$, **assuming that the 4-metric is static, so that $\\partial_t g_{\\mu\\nu}=0$ for all $\\mu$ and $\\nu$**.\n",
    "\n",
    "To compute $g^{\\mu\\nu}$, we use the standard formula (Eq. 4.49 in [Gourgoulhon](https://arxiv.org/pdf/gr-qc/0703035.pdf)):\n",
    "$$\n",
    "g^{\\mu\\nu} = \\begin{pmatrix} \n",
    "-\\frac{1}{\\alpha^2} & \\frac{\\beta^i}{\\alpha^2} \\\\\n",
    "\\frac{\\beta^i}{\\alpha^2} & \\gamma^{ij} - \\frac{\\beta^i\\beta^j}{\\alpha^2}\n",
    "\\end{pmatrix},\n",
    "$$\n",
    "where $\\gamma^{ij}$ is given by the inverse of $\\gamma_{ij}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:12.865137Z",
     "iopub.status.busy": "2021-03-07T17:07:12.839188Z",
     "iopub.status.idle": "2021-03-07T17:07:12.882955Z",
     "shell.execute_reply": "2021-03-07T17:07:12.882336Z"
    }
   },
   "outputs": [],
   "source": [
    "betaDdD = ixp.zerorank2()\n",
    "\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            # Recall that betaD[i] = gammaDD[i][j]*betaU[j] (Eq. 2.121 in B&S)\n",
    "            betaDdD[i][k] += gammaDD_dD[i][j][k]*betaU[j] + gammaDD[i][j]*betaU_dD[j][k]\n",
    "\n",
    "# Eq. 2.122 in B&S\n",
    "g4DDdD = ixp.zerorank3(DIM=4)\n",
    "for i in range(DIM):\n",
    "    # Recall that g4DD[0][0] = -alpha^2 + betaU[i]*betaD[i]\n",
    "    g4DDdD[0][0][i+1] += -2*alpha*alpha_dD[i]\n",
    "    for j in range(DIM):\n",
    "        g4DDdD[0][0][i+1] += betaU_dD[j][i]*betaD[j] + betaU[j]*betaDdD[j][i]\n",
    "\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        # Recall that g4DD[i][0] = g4DD[0][i] = betaD[i]\n",
    "        g4DDdD[i+1][0][j+1] = g4DDdD[0][i+1][j+1] = betaDdD[i][j]\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            # Recall that g4DD[i][j] = gammaDD[i][j]\n",
    "            g4DDdD[i+1][j+1][k+1] = gammaDD_dD[i][j][k]\n",
    "\n",
    "gammaUU, dummyDET = ixp.symm_matrix_inverter3x3(gammaDD)\n",
    "\n",
    "g4UU = ixp.zerorank2(DIM=4)\n",
    "g4UU[0][0] = -1 / alpha**2\n",
    "for i in range(DIM):\n",
    "    g4UU[0][i+1] = g4UU[i+1][0] = betaU[i]/alpha**2\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        g4UU[i+1][j+1] = gammaUU[i][j] - betaU[i]*betaU[j]/alpha**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we are to compute:\n",
    "$$\n",
    "\\Gamma^{\\mu}_{\\nu\\delta} = \\frac{1}{2} g^{\\mu\\eta} \\left(g_{\\eta\\nu,\\delta} + g_{\\eta\\delta,\\nu} - g_{\\nu\\delta,\\eta}  \\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:12.957431Z",
     "iopub.status.busy": "2021-03-07T17:07:12.921585Z",
     "iopub.status.idle": "2021-03-07T17:07:27.453676Z",
     "shell.execute_reply": "2021-03-07T17:07:27.453029Z"
    }
   },
   "outputs": [],
   "source": [
    "Gamma4UDD = ixp.zerorank3(DIM=4)\n",
    "for mu in range(4):\n",
    "    for nu in range(4):\n",
    "        for delta in range(4):\n",
    "            for eta in range(4):\n",
    "                Gamma4UDD[mu][nu][delta] += sp.Rational(1,2)*g4UU[mu][eta]*\\\n",
    "                (g4DDdD[eta][nu][delta] + g4DDdD[eta][delta][nu] - g4DDdD[nu][delta][eta])\n",
    "\n",
    "# Now output the 4-Christoffels to file:\n",
    "for mu in range(4):\n",
    "    for nu in range(4):\n",
    "        for delta in range(nu,4):\n",
    "            gf_interp_list.append(gf_interp(\"4-Christoffel GammaUDD\"+str(mu)+str(nu)+str(delta)))\n",
    "            interp_expr = Gamma4UDD[mu][nu][delta]\n",
    "            which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy4christoffels_spherical'></a>\n",
    "\n",
    "### Step 2.a.v: Notes on computing all 40 4-Christoffels $\\Gamma^{\\mu}_{\\nu\\delta}$ in the Spherical basis \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy4christoffels_spherical}$$\n",
    "\n",
    "As explained in [Eq. 3.15 of Carroll's lecture notes on GR](https://ned.ipac.caltech.edu/level5/March01/Carroll3/Carroll3.html), while connection coefficients (a.k.a. Christoffel symbols) are not tensors, differences in connection coefficients are tensors.\n",
    "\n",
    "Thus we may define\n",
    "\n",
    "$$\n",
    "\\Delta^\\mu_{\\nu\\delta} = \\Gamma^\\mu_{\\nu\\delta} - \\hat{\\Gamma}^\\mu_{\\nu\\delta},\n",
    "$$\n",
    "\n",
    "where for example $\\Gamma^\\mu_{\\nu\\delta}$ is the connection related to the curved spacetime 4-metric in some basis and $\\hat{\\Gamma}^\\mu_{\\nu\\delta}$ is the connection related to the flat spacetime 4-metric in the same basis.\n",
    "\n",
    "We are given the 4-metric data in Cartesian coordinates, for which $\\hat{\\Gamma}^\\mu_{\\nu\\delta}=0$. The basis transform to spherical coordinates is then straightforward:\n",
    "\n",
    "\\begin{align}\n",
    "\\Delta^\\mu_{\\text{Sph}\\ \\nu\\delta} &= \\Gamma^\\mu_{\\text{Sph}\\ \\nu\\delta} - \\hat{\\Gamma}^\\mu_{\\text{Sph}\\ \\nu\\delta} \\\\\n",
    "&= \\frac{\\partial x^\\mu_{\\rm Sph}}{\\partial x^\\alpha_{\\rm Cart}}\n",
    "\\frac{\\partial x^\\beta_{\\rm Cart}}{\\partial x^\\nu_{\\rm Sph}}\n",
    "\\frac{\\partial x^\\gamma_{\\rm Cart}}{\\partial x^\\delta_{\\rm Sph}} \\Delta^\\alpha_{\\text{Cart}\\ \\beta\\gamma} \\\\\n",
    "&= \\frac{\\partial x^\\mu_{\\rm Sph}}{\\partial x^\\alpha_{\\rm Cart}}\n",
    "\\frac{\\partial x^\\beta_{\\rm Cart}}{\\partial x^\\nu_{\\rm Sph}}\n",
    "\\frac{\\partial x^\\gamma_{\\rm Cart}}{\\partial x^\\delta_{\\rm Sph}} \\Gamma^\\alpha_{\\text{Cart}\\ \\beta\\gamma} \\\\\n",
    "\\implies \\Gamma^\\mu_{\\text{Sph}\\ \\nu\\delta} &= \\frac{\\partial x^\\mu_{\\rm Sph}}{\\partial x^\\alpha_{\\rm Cart}}\n",
    "\\frac{\\partial x^\\beta_{\\rm Cart}}{\\partial x^\\nu_{\\rm Sph}}\n",
    "\\frac{\\partial x^\\gamma_{\\rm Cart}}{\\partial x^\\delta_{\\rm Sph}} \\Gamma^\\alpha_{\\text{Cart}\\ \\beta\\gamma} +\n",
    "\\hat{\\Gamma}^\\mu_{\\text{Sph}\\ \\nu\\delta}\n",
    "\\end{align}\n",
    "\n",
    "**Define $\\hat{\\Gamma}^\\mu_{\\text{Sph}\\ \\nu\\delta}$.**\n",
    "\n",
    "By definition,\n",
    "$$\n",
    "\\hat{\\Gamma}^{\\mu}_{\\nu\\delta} = \\frac{1}{2} \\hat{g}^{\\mu\\eta} \\left(\\hat{g}_{\\eta\\nu,\\delta} + \\hat{g}_{\\eta\\delta,\\nu} - \\hat{g}_{\\nu\\delta,\\eta}  \\right).\n",
    "$$\n",
    "\n",
    "In static spherical coordinates, $\\hat{g}_{\\nu\\delta}$ is given by \n",
    "\n",
    "$$\n",
    "\\hat{g}_{\\mu\\nu} = \\begin{pmatrix} \n",
    "-1 & 0 \\\\\n",
    "0 & \\hat{\\gamma}_{ij}\n",
    "\\end{pmatrix},\n",
    "$$\n",
    "so the inverse is easy to compute:\n",
    "$$\n",
    "\\hat{g}^{\\mu\\nu} = \\begin{pmatrix} \n",
    "-1 & 0 \\\\\n",
    "0 & 1/\\hat{\\gamma}_{ij}\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "Here is the NRPy+ code implementation of $\\hat{g}_{\\mu\\nu}$, $\\hat{g}^{\\mu\\nu}$, and $\\hat{g}_{\\eta\\nu,\\delta}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.458692Z",
     "iopub.status.busy": "2021-03-07T17:07:27.457530Z",
     "iopub.status.idle": "2021-03-07T17:07:27.460414Z",
     "shell.execute_reply": "2021-03-07T17:07:27.461049Z"
    }
   },
   "outputs": [],
   "source": [
    "# import reference_metric as rfm\n",
    "# # Set the desired *output* coordinate system to Spherical:\n",
    "# #par.set_parval_from_str(\"reference_metric::CoordSystem\",\"NobleSphericalThetaOptionOne\")\n",
    "# par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Spherical\")\n",
    "# print(\"Calling reference_metric()...\")\n",
    "# rfm.reference_metric()\n",
    "# print(\"Just finished calling reference_metric()...\")\n",
    "\n",
    "\n",
    "# g4hatDD = ixp.zerorank2(DIM=4)\n",
    "# g4hatUU = ixp.zerorank2(DIM=4)\n",
    "# g4hatDD[0][0] = sp.sympify(-1)\n",
    "# g4hatUU[0][0] = sp.sympify(-1)\n",
    "# for j in range(3):\n",
    "#     g4hatDD[j+1][j+1] = rfm.ghatDD[j][j]\n",
    "#     g4hatUU[j+1][j+1] = 1/rfm.ghatDD[j][j]\n",
    "# g4hatDDdD = ixp.zerorank3(DIM=4)\n",
    "# for eta in range(4):\n",
    "#     for nu in range(4):\n",
    "#         for j in range(3): # Time derivatives are all zero, so g4hatDDdD[eta][nu][0] = 0 (as initialized).\n",
    "#             g4hatDDdD[eta][nu][j+1] = sp.diff(g4hatDD[eta][nu],rfm.xx[j])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we compute the 4-Christoffels $\\hat{\\Gamma}^\\mu_{\\text{Sph}\\ \\nu\\delta}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.466156Z",
     "iopub.status.busy": "2021-03-07T17:07:27.465072Z",
     "iopub.status.idle": "2021-03-07T17:07:27.467655Z",
     "shell.execute_reply": "2021-03-07T17:07:27.468242Z"
    }
   },
   "outputs": [],
   "source": [
    "# Gamma4hatSphUDD = ixp.zerorank3(DIM=4)\n",
    "# for mu in range(4):\n",
    "#     for nu in range(4):\n",
    "#         for delta in range(4):\n",
    "#             for eta in range(4):\n",
    "#                 Gamma4hatSphUDD[mu][nu][delta] += sp.Rational(1,2)*g4hatUU[mu][eta]* \\\n",
    "#                 ( g4hatDDdD[eta][nu][delta] + g4hatDDdD[eta][delta][nu] - g4hatDDdD[nu][delta][eta] )\n",
    "\n",
    "# # Here are the results, cf. Eq 18 of https://arxiv.org/pdf/1211.6632.pdf\n",
    "# sp.pretty_print(Gamma4hatSphUDD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, compute $\\Gamma^\\mu_{\\text{Sph}\\ \\nu\\delta}$. Recall from above that\n",
    "\\begin{align}\n",
    "\\Gamma^\\mu_{\\text{Sph}\\ \\nu\\delta} &= \\frac{\\partial x^\\mu_{\\rm Sph}}{\\partial x^\\alpha_{\\rm Cart}}\n",
    "\\frac{\\partial x^\\beta_{\\rm Cart}}{\\partial x^\\nu_{\\rm Sph}}\n",
    "\\frac{\\partial x^\\gamma_{\\rm Cart}}{\\partial x^\\delta_{\\rm Sph}} \\Gamma^\\alpha_{\\text{Cart}\\ \\beta\\gamma} +\n",
    "\\hat{\\Gamma}^\\mu_{\\text{Sph}\\ \\nu\\delta}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpybasisxform'></a>\n",
    "\n",
    "### Step 2.a.vi:  Notes on basis transforming all Cartesian basis quantities to spherical \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpybasisxform}$$\n",
    "\n",
    "All tensors and vectors are in the Cartesian coordinate basis $x^i_{\\rm Cart} = (x,y,z)$, but we need them in the curvilinear coordinate basis $x^i_{\\rm rfm}$=`(xx0,xx1,xx2)`=$(r,\\theta,\\phi)$ set by NRPy+'s `\"reference_metric::CoordSystem\"` variable (we'll set this parameter to `\"Spherical\"`). \n",
    "\n",
    "Empirically speaking, it is usually easier to write `(x(xx0,xx1,xx2),y(xx0,xx1,xx2),z(xx0,xx1,xx2))` than the inverse, so we will compute the Jacobian matrix\n",
    "\n",
    "$$\n",
    "{\\rm Jac\\_dUSph\\_dDrfmUD[i][j]} = \\frac{\\partial x^i_{\\rm Cart}}{\\partial x^j_{\\rm rfm}},\n",
    "$$\n",
    "\n",
    "via exact differentiation (courtesy SymPy), and the inverse Jacobian\n",
    "$$\n",
    "{\\rm Jac\\_dUrfm\\_dDSphUD[i][j]} = \\frac{\\partial x^i_{\\rm rfm}}{\\partial x^j_{\\rm Cart}},\n",
    "$$\n",
    "\n",
    "using NRPy+'s `generic\\_matrix\\_inverter3x3()` function. In terms of these, the transformation of vectors and rank-2 fully covariant tensors from Cartesian to `\"reference_metric::CoordSystem\"` (Spherical) coordinates may be written:\n",
    "\n",
    "\\begin{align}\n",
    "g^{\\rm rfm}_{\\mu\\nu} &= \n",
    "\\frac{\\partial x^{\\alpha}_{\\rm Cart}}{\\partial x^{\\mu}_{\\rm rfm}}\n",
    "\\frac{\\partial x^{\\beta}_{\\rm Cart}}{\\partial x^{\\nu}_{\\rm rfm}} g^{\\rm Cart}_{\\alpha \\beta} \\\\\n",
    "\\Gamma^\\mu_{\\text{Sph}\\ \\nu\\delta} &= \\frac{\\partial x^\\mu_{\\rm Sph}}{\\partial x^\\alpha_{\\rm Cart}}\n",
    "\\frac{\\partial x^\\beta_{\\rm Cart}}{\\partial x^\\nu_{\\rm Sph}}\n",
    "\\frac{\\partial x^\\gamma_{\\rm Cart}}{\\partial x^\\delta_{\\rm Sph}} \\Gamma^\\alpha_{\\text{Cart}\\ \\beta\\gamma} +\n",
    "\\hat{\\Gamma}^\\mu_{\\text{Sph}\\ \\nu\\delta}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.474628Z",
     "iopub.status.busy": "2021-03-07T17:07:27.473362Z",
     "iopub.status.idle": "2021-03-07T17:07:27.476413Z",
     "shell.execute_reply": "2021-03-07T17:07:27.476980Z"
    }
   },
   "outputs": [],
   "source": [
    "# Jac_dUCart_dDrfmUD = ixp.zerorank2()\n",
    "# for i in range(DIM):\n",
    "#     for j in range(DIM):\n",
    "#         Jac_dUCart_dDrfmUD[i][j] = sp.simplify(sp.diff(rfm.xx_to_Cart[i],rfm.xx[j]))\n",
    "\n",
    "# Jac_dUrfm_dDCartUD, dummyDET = ixp.generic_matrix_inverter3x3(Jac_dUCart_dDrfmUD)\n",
    "# Jac4_dUCart_dDrfmUD = ixp.zerorank2(DIM=4)\n",
    "# Jac4_dUrfm_dDCartUD = ixp.zerorank2(DIM=4)\n",
    "# for alp in range(4):\n",
    "#     for bet in range(4):\n",
    "#         if alp==0 or bet==0:\n",
    "#             Jac4_dUCart_dDrfmUD[alp][bet] = sp.sympify(1) # Time components unchanged\n",
    "#             Jac4_dUrfm_dDCartUD[alp][bet] = sp.sympify(1) # Time components unchanged\n",
    "#         else:\n",
    "#             Jac4_dUCart_dDrfmUD[alp][bet] = sp.simplify(Jac_dUCart_dDrfmUD[alp-1][bet-1])\n",
    "#             Jac4_dUrfm_dDCartUD[alp][bet] = sp.simplify(Jac_dUrfm_dDCartUD[alp-1][bet-1])\n",
    "\n",
    "# Gamma4SphUDD = ixp.zerorank3(DIM=4)\n",
    "# for mu in range(4):\n",
    "#     for nu in range(4):\n",
    "#         for delt in range(4):\n",
    "#             Gamma4SphUDD[mu][nu][delt] = Gamma4hatSphUDD[mu][nu][delt]\n",
    "#             for alp in range(4):\n",
    "#                 for bet in range(4):\n",
    "#                     for gam in range(4):\n",
    "#                         Gamma4SphUDD[mu][nu][delt] += \\\n",
    "#                 Jac4_dUrfm_dDCartUD[mu][alp]*Jac4_dUCart_dDrfmUD[bet][nu]*Jac4_dUCart_dDrfmUD[gam][delt] * \\\n",
    "#                           Gamma4UDD[alp][bet][gam]\n",
    "\n",
    "# # Now output the Spherical 4-Christoffels to file:\n",
    "# for mu in range(4):\n",
    "#     for nu in range(4):\n",
    "#         for delt in range(nu,4):\n",
    "#             gf_interp_list.append(gf_interp(\"4-Christoffel component in SPHERICAL BASIS: GammaSphUDD\"+str(mu)+str(nu)+str(delt)))\n",
    "#             interp_expr = Gamma4SphUDD[mu][nu][delt]\n",
    "#             which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Output the 4-metric in Spherical coordinates to file:\n",
    "\n",
    "\\begin{align}\n",
    "g^{\\rm rfm}_{\\mu\\nu} &= \n",
    "\\frac{\\partial x^{\\alpha}_{\\rm Cart}}{\\partial x^{\\mu}_{\\rm rfm}}\n",
    "\\frac{\\partial x^{\\beta}_{\\rm Cart}}{\\partial x^{\\nu}_{\\rm rfm}} g^{\\rm Cart}_{\\alpha \\beta}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.481376Z",
     "iopub.status.busy": "2021-03-07T17:07:27.480466Z",
     "iopub.status.idle": "2021-03-07T17:07:27.483441Z",
     "shell.execute_reply": "2021-03-07T17:07:27.482808Z"
    }
   },
   "outputs": [],
   "source": [
    "# g4SphDD = ixp.zerorank2(DIM=4)\n",
    "# for mu in range(4):\n",
    "#     for nu in range(4):\n",
    "#         for alp in range(4):\n",
    "#             for bet in range(4):\n",
    "#                 g4SphDD[mu][nu] += Jac4_dUCart_dDrfmUD[alp][mu]*Jac4_dUCart_dDrfmUD[bet][nu]*g4DD[alp][bet]\n",
    "\n",
    "# for mu in range(4):\n",
    "#     for nu in range(mu,4):\n",
    "#         gf_interp_list.append(gf_interp(\"4-metric component in SPHERICAL BASIS: g4SphDD\"+str(mu)+str(nu)))\n",
    "#         interp_expr = g4SphDD[mu][nu]\n",
    "#         which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next output the various GRMHD 3-vectors in the Spherical basis\n",
    "\n",
    "\\begin{align}\n",
    "v^i_{\\rm Sph} &= \n",
    "\\frac{\\partial x^{i}_{\\rm Sph}}{\\partial x^{j}_{\\rm Cart}} v^j_{\\rm Cart}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.488612Z",
     "iopub.status.busy": "2021-03-07T17:07:27.487505Z",
     "iopub.status.idle": "2021-03-07T17:07:27.490163Z",
     "shell.execute_reply": "2021-03-07T17:07:27.490875Z"
    }
   },
   "outputs": [],
   "source": [
    "# IGMvSphU                  = ixp.zerorank1()\n",
    "# ValenciavSphU             = ixp.zerorank1()\n",
    "# Gamma_times_ValenciavSphU = ixp.zerorank1()\n",
    "# BSphU                     = ixp.zerorank1()\n",
    "\n",
    "# for i in range(DIM):\n",
    "#     for j in range(DIM):\n",
    "#         IGMvSphU[i]                  += Jac_dUrfm_dDCartUD[i][j] * IGMvU[j]\n",
    "#         ValenciavSphU[i]             += Jac_dUrfm_dDCartUD[i][j] * ValenciavU[j]\n",
    "#         Gamma_times_ValenciavSphU[i] += Jac_dUrfm_dDCartUD[i][j] * Gamma_times_ValenciavU[j]\n",
    "#         BSphU[i]                     += Jac_dUrfm_dDCartUD[i][j] * BU[j]\n",
    "\n",
    "# for i in range(DIM):\n",
    "#     gf_interp_list.append(gf_interp(\"IGM 3-velocity vU\"+str(i)+\" = u^i/u^0 in SPHERICAL BASIS\"))\n",
    "#     interp_expr = IGMvSphU[i]\n",
    "#     which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "# for i in range(DIM):\n",
    "#     gf_interp_list.append(gf_interp(\"Valencia 3-velocity vU\"+str(i)+\" in SPHERICAL BASIS\"))\n",
    "#     interp_expr = ValenciavSphU[i]\n",
    "#     which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "# for i in range(DIM):\n",
    "#     gf_interp_list.append(gf_interp(\"Lorentz factor, times Valencia vU\"+str(i)+\" in SPHERICAL BASIS\"))\n",
    "#     interp_expr = Gamma_times_ValenciavSphU[i]\n",
    "#     which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "# for i in range(DIM):\n",
    "#     gf_interp_list.append(gf_interp(\"IGM magnetic field component B\"+str(i)+\" in SPHERICAL BASIS\"))\n",
    "#     interp_expr = BSphU[i]\n",
    "#     which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " <a id='psi4andfriends'></a>\n",
    "\n",
    "### Step 2.a.vii: ): Output Weyl scalars $\\psi_0$ through $\\psi_4$, as well as Weyl invariants $J$ and $I$, from the `WeylScal4` ETK thorn \\[Back to [top](#toc)\\]\n",
    "$$\\label{psi4andfriends}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.505221Z",
     "iopub.status.busy": "2021-03-07T17:07:27.504407Z",
     "iopub.status.idle": "2021-03-07T17:07:27.550079Z",
     "shell.execute_reply": "2021-03-07T17:07:27.549512Z"
    }
   },
   "outputs": [],
   "source": [
    "Weylgfs = [\"Psi0r\",\"Psi0i\",\"Psi1r\",\"Psi1i\",\"Psi2r\",\"Psi2i\",\"Psi3r\",\"Psi3i\",\"Psi4r\",\"Psi4i\",\n",
    "                                    \"curvIr\",\"curvIi\",\"curvJr\",\"curvJi\"]\n",
    "Psi0r,Psi0i,Psi1r,Psi1i,Psi2r,Psi2i,Psi3r,Psi3i,Psi4r,Psi4i,curvIr,curvIi,curvJr,curvJi = \\\n",
    "    gri.register_gridfunctions(\"AUX\",Weylgfs)\n",
    "count = 0\n",
    "for gf in [Psi0r,Psi0i,Psi1r,Psi1i,Psi2r,Psi2i,Psi3r,Psi3i,Psi4r,Psi4i,curvIr,curvIi,curvJr,curvJi]:\n",
    "    gf_interp_list.append(gf_interp(\"4-Weyl scalar or invariant \"+Weylgfs[count]))\n",
    "    interp_expr = gf\n",
    "    which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "    count = count + 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='constraints_gfs'></a>\n",
    "\n",
    "### Step 2.a.viii: Hamiltonian and momentum constraints \\[Back to [top](#toc)\\]\n",
    "$$\\label{constraints_gfs}$$\n",
    "\n",
    "We now setup the code to output the Hamiltonian and momentum constraint gridfunctions, that is:\n",
    "\n",
    "* $\\mathcal{H}$ (e.g., $\\text{ML_BSSN::H}$ or $\\text{Baikal::HGF}$)\n",
    "* $\\mathcal{M}^{0}$ (e.g., $\\text{ML_BSSN::M1}$ or $\\text{Baikal::MU0GF}$)\n",
    "* $\\mathcal{M}^{1}$ (e.g., $\\text{ML_BSSN::M2}$ or $\\text{Baikal::MU1GF}$)\n",
    "* $\\mathcal{M}^{2}$ (e.g., $\\text{ML_BSSN::M3}$ or $\\text{Baikal::MU2GF}$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "gf_interp_list.append(gf_interp(\"Constraints - Hamiltonian\"))\n",
    "H = gri.register_gridfunctions(\"AUX\",\"InterpH\")\n",
    "interp_expr = H\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "MU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"InterpMU\")\n",
    "\n",
    "# Sph_basis_MU = ixp.zerorank1()\n",
    "# for i in range(DIM):\n",
    "#     for l in range(DIM):\n",
    "#         Sph_basis_MU[i] += Jac_dUrfm_dDCartUD[i][l] * MU[l]\n",
    "\n",
    "for i in range(DIM):\n",
    "    gf_interp_list.append(gf_interp(\"Constraints - Momentum M^\"+chr(ord('x')+i)))\n",
    "#     interp_expr = Sph_basis_MU[i]\n",
    "    interp_expr = MU[i]\n",
    "    which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nuc_eos_gfs'></a>\n",
    "\n",
    "### Step 2.a.ix: Nuclear (tabulated) equation of state quantities \\[Back to [top](#toc)\\]\n",
    "$$\\label{nuc_eos_gfs}$$\n",
    "\n",
    "When using nuclear (tabulated) equations of state, additional gridfunctions are available for output. These are:\n",
    "\n",
    "* $T$, the temperature (e.g. $\\text{HydroBase::temperature}$ or $\\text{IllinoisGRMHD::igm_temperature}$)\n",
    "* $Y_{e}$, the electron fraction (e.g. $\\text{HydroBase::Y_e}$ or $\\text{IllinoisGRMHD::igm_Ye}$)\n",
    "* $\\epsilon$, the specific internal energy (e.g. $\\text{HydroBase::eps}$ or $\\text{IllinoisGRMHD::igm_eps}$)\n",
    "* $S$, the entropy (e.g. $\\text{HydroBase::entropy}$ or $\\text{IllinoisGRMHD::igm_entropy}$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "gf_interp_list.append(gf_interp(\"Temperature primitive\"))\n",
    "temperature = gri.register_gridfunctions(\"AUX\",\"temperature\")\n",
    "interp_expr = temperature\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "gf_interp_list.append(gf_interp(\"Electron fraction primitive\"))\n",
    "Ye = gri.register_gridfunctions(\"AUX\",\"Y_e\")\n",
    "interp_expr = Ye\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "gf_interp_list.append(gf_interp(\"Specific internal energy primitive\"))\n",
    "eps = gri.register_gridfunctions(\"AUX\",\"eps\")\n",
    "interp_expr = eps\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "gf_interp_list.append(gf_interp(\"Entropy primitive\"))\n",
    "entropy = gri.register_gridfunctions(\"AUX\",\"entropy\")\n",
    "interp_expr = entropy\n",
    "which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy_c_calling_function'></a>\n",
    "\n",
    "## Step 2.b: C code calling function for the NRPy+ C output \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy_c_calling_function}$$\n",
    "\n",
    "In the above blocks, we wrote and appended to a file `list_of_functions_to_interpolate.h`. Here we write the calling function for this C code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.556180Z",
     "iopub.status.busy": "2021-03-07T17:07:27.555295Z",
     "iopub.status.idle": "2021-03-07T17:07:27.559757Z",
     "shell.execute_reply": "2021-03-07T17:07:27.559097Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/src/construct_function_to_interpolate__store_to_interped_gf.cc\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/construct_function_to_interpolate__store_to_interped_gf.cc\n",
    "#include <stdio.h>\n",
    "#include <stdlib.h>\n",
    "#include \"cctk.h\"\n",
    "#include \"cctk_Arguments.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "\n",
    "// Set the gridfunction interped_gf, according to the interpolation counter variable interp_counter.\n",
    "//    For example, we might interpolate \"IllinoisGRMHD::rho_b\" if interp_counter==0. The following\n",
    "//    function takes care of these\n",
    "void list_of_functions_to_interpolate(cGH *restrict cctkGH,\n",
    "                                      const CCTK_INT *restrict cctk_lsh,\n",
    "                                      const CCTK_INT *restrict cctk_nghostzones,\n",
    "                                      const CCTK_REAL *restrict x,const CCTK_REAL *restrict y,const CCTK_REAL *restrict z,\n",
    "                                      const CCTK_REAL invdx0,const CCTK_REAL invdx1,const CCTK_REAL invdx2,\n",
    "                                      const CCTK_INT *restrict InterpCounter,\n",
    "                                      const CCTK_REAL *restrict rho_bGF,\n",
    "                                      const CCTK_REAL *restrict PGF,\n",
    "                                      const CCTK_REAL *restrict temperatureGF,\n",
    "                                      const CCTK_REAL *restrict Y_eGF,\n",
    "                                      const CCTK_REAL *restrict epsGF,\n",
    "                                      const CCTK_REAL *restrict entropyGF,\n",
    "                                      const CCTK_REAL *restrict IGMvU0GF,const CCTK_REAL *restrict IGMvU1GF,const CCTK_REAL *restrict IGMvU2GF,\n",
    "                                      const CCTK_REAL *restrict BU0GF,const CCTK_REAL *restrict BU1GF,const CCTK_REAL *restrict BU2GF,\n",
    "                                      const CCTK_REAL *restrict gammaDD00GF,const CCTK_REAL *restrict gammaDD01GF,const CCTK_REAL *restrict gammaDD02GF,\n",
    "                                      const CCTK_REAL *restrict gammaDD11GF,const CCTK_REAL *restrict gammaDD12GF,const CCTK_REAL *restrict gammaDD22GF,\n",
    "                                      const CCTK_REAL *restrict betaU0GF,const CCTK_REAL *restrict betaU1GF,const CCTK_REAL *restrict betaU2GF,\n",
    "                                      const CCTK_REAL *restrict alphaGF,\n",
    "                                      CCTK_REAL *restrict interped_gfGF,\n",
    "                                      CCTK_REAL *restrict Ax_unstaggeredGF,CCTK_REAL *restrict Ay_unstaggeredGF,CCTK_REAL *restrict Az_unstaggeredGF,\n",
    "                                      const CCTK_REAL *restrict Psi0rGF,const CCTK_REAL *restrict Psi0iGF,\n",
    "                                      const CCTK_REAL *restrict Psi1rGF,const CCTK_REAL *restrict Psi1iGF,\n",
    "                                      const CCTK_REAL *restrict Psi2rGF,const CCTK_REAL *restrict Psi2iGF,\n",
    "                                      const CCTK_REAL *restrict Psi3rGF,const CCTK_REAL *restrict Psi3iGF,\n",
    "                                      const CCTK_REAL *restrict Psi4rGF,const CCTK_REAL *restrict Psi4iGF,\n",
    "                                      const CCTK_REAL *restrict curvIrGF,const CCTK_REAL *restrict curvIiGF,\n",
    "                                      const CCTK_REAL *restrict curvJrGF,const CCTK_REAL *restrict curvJiGF,\n",
    "                                      const CCTK_REAL *restrict InterpHGF,\n",
    "                                      const CCTK_REAL *restrict InterpMU0GF,\n",
    "                                      const CCTK_REAL *restrict InterpMU1GF,\n",
    "                                      const CCTK_REAL *restrict InterpMU2GF) {\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "#include \"list_of_functions_to_interpolate.h\"\n",
    "}\n",
    "\n",
    "void construct_function_to_interpolate__store_to_interped_gf(CCTK_ARGUMENTS) {\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  printf(\"Called construct_function_to_interpolate__store_to_interped_gf() on grid with dx = %e!\\n\",CCTK_DELTA_SPACE(0));\n",
    "  const CCTK_REAL invdx0 = 1.0 / CCTK_DELTA_SPACE(0);\n",
    "  const CCTK_REAL invdx1 = 1.0 / CCTK_DELTA_SPACE(1);\n",
    "  const CCTK_REAL invdx2 = 1.0 / CCTK_DELTA_SPACE(2);\n",
    "\n",
    "  // .------------------------------------.\n",
    "  // | Constraint equations gridfunctions |\n",
    "  // .------------------------------------.\n",
    "  // We get the pointer using the CCTK_VarDataPtr() function because\n",
    "  // if the gridfunction does not exist in this version of IGM then\n",
    "  // we will get compilation errors.\n",
    "  int timelevel = 0;\n",
    "\n",
    "  // Hamiltonian constraint gridfunction\n",
    "  CCTK_REAL *InterpHGF = (CCTK_REAL*)(CCTK_VarDataPtr(cctkGH,timelevel,HVarString));\n",
    "  if( !InterpHGF ) CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, \"Hamiltonian constraint - Couldn't get data pointer of input array variable '%s'\", HVarString);\n",
    "\n",
    "  // Momentum constraint gridfunction (x-direction)\n",
    "  CCTK_REAL *InterpMU0GF = (CCTK_REAL*)(CCTK_VarDataPtr(cctkGH,timelevel,M0VarString));\n",
    "  if( !InterpMU0GF ) CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, \"Momentum constraint (x-direction) - Couldn't get data pointer of input array variable '%s'\", M0VarString);\n",
    "\n",
    "  // Momentum constraint gridfunction (y-direction)\n",
    "  CCTK_REAL *InterpMU1GF = (CCTK_REAL*)(CCTK_VarDataPtr(cctkGH,timelevel,M1VarString));\n",
    "  if( !InterpMU1GF ) CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, \"Momentum constraint (y-direction) - Couldn't get data pointer of input array variable '%s'\", M1VarString);\n",
    "\n",
    "  // Momentum constraint gridfunction (z-direction)\n",
    "  CCTK_REAL *InterpMU2GF = (CCTK_REAL*)(CCTK_VarDataPtr(cctkGH,timelevel,M2VarString));\n",
    "  if( !InterpMU2GF ) CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, \"Momentum constraint (z-direction) - Couldn't get data pointer of input array variable '%s'\", M2VarString);\n",
    "\n",
    "  // .-------------------------------------------------------.\n",
    "  // | Additional hydro quantities for nuclear/tabulated EOS |\n",
    "  // .-------------------------------------------------------.\n",
    "  CCTK_REAL *epsGF, *temperatureGF, *Y_eGF, *entropyGF;\n",
    "  if( enable_nuc_eos ) {\n",
    "    // Temperature gridfunction\n",
    "    temperatureGF = (CCTK_REAL*)(CCTK_VarDataPtr(cctkGH,timelevel,temperatureVarString));\n",
    "    if( !temperatureGF ) CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, \"temperature - Couldn't get data pointer of input array variable '%s'\", temperatureVarString);\n",
    "\n",
    "    // Electron fraction gridfunction\n",
    "    Y_eGF = (CCTK_REAL*)(CCTK_VarDataPtr(cctkGH,timelevel,Y_eVarString));\n",
    "    if( !Y_eGF ) CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, \"Y_e - Couldn't get data pointer of input array variable '%s'\", Y_eVarString);\n",
    "\n",
    "    // Specific internal energy gridfunction\n",
    "    epsGF = (CCTK_REAL*)(CCTK_VarDataPtr(cctkGH,timelevel,epsVarString));\n",
    "    if( !epsGF ) CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, \"eps - Couldn't get data pointer of input array variable '%s'\", epsVarString);\n",
    "\n",
    "    // Entropy gridfunction\n",
    "    entropyGF = (CCTK_REAL*)(CCTK_VarDataPtr(cctkGH,timelevel,entropyVarString));\n",
    "    if( !entropyGF ) CCTK_VError(__LINE__, __FILE__, CCTK_THORNSTRING, \"entropy - Couldn't get data pointer of input array variable '%s'\", entropyVarString);\n",
    "  }\n",
    "  else {\n",
    "    // If we are not using nuclear EOS, then set the gridfunction pointers to NULL.\n",
    "    temperatureGF = NULL;\n",
    "    Y_eGF         = NULL;\n",
    "    epsGF         = NULL;\n",
    "    entropyGF     = NULL;\n",
    "  }\n",
    "\n",
    "  list_of_functions_to_interpolate(cctkGH,cctk_lsh,cctk_nghostzones,\n",
    "                                   x,y,z,\n",
    "                                   invdx0,invdx1,invdx2,\n",
    "                                   InterpCounter,\n",
    "                                   rho_b,P,temperatureGF,Y_eGF,epsGF,entropyGF,\n",
    "                                   vx,vy,vz,\n",
    "                                   Bx,By,Bz,\n",
    "                                   gxx,gxy,gxz,gyy,gyz,gzz,\n",
    "                                   betax,betay,betaz,alp, interped_gf,\n",
    "                                   Ax_unstaggered,Ay_unstaggered,Az_unstaggered,\n",
    "                                   Psi0r,Psi0i,Psi1r,Psi1i,Psi2r,Psi2i,Psi3r,Psi3i,Psi4r,Psi4i,\n",
    "                                   curvIr,curvIi,curvJr,curvJi,\n",
    "                                   InterpHGF,InterpMU0GF,InterpMU1GF,InterpMU2GF);\n",
    "\n",
    "  // interped_gf will be interpolated across AMR boundaries, meaning that\n",
    "  //    it must be prolongated. Only gridfunctions with 3 timelevels stored\n",
    "  //    may be prolongated (provided time_interpolation_order is set to the\n",
    "  //    usual value of 2). We should only call this interpolation routine\n",
    "  //    at iterations in which all gridfunctions are on the same timelevel\n",
    "  //    (usually a power of 2), which will ensure that the following\n",
    "  //    \"filling of the timelevels\" is completely correct.\n",
    "#pragma omp parallel for\n",
    "  for(int i=0;i<cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2];i++) {\n",
    "    interped_gf_p[i]   = interped_gf[i];\n",
    "    interped_gf_p_p[i] = interped_gf[i];\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpygetgfname'></a>\n",
    "\n",
    "## Step 2.c: The `get_gf_name()` function \\[Back to [top](#toc)\\]\n",
    "\n",
    "$$\\label{nrpygetgfname}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.567119Z",
     "iopub.status.busy": "2021-03-07T17:07:27.566157Z",
     "iopub.status.idle": "2021-03-07T17:07:27.568582Z",
     "shell.execute_reply": "2021-03-07T17:07:27.569199Z"
    }
   },
   "outputs": [],
   "source": [
    "with open(os.path.join(Ccodesdir,\"src\",\"get_gf_name.h\"), \"w\") as file:\n",
    "    file.write(\"void get_gf_name(const int InterpCounter,char gf_name[100]) {\\n\")\n",
    "    for i in range(1,which_InterpCounter):\n",
    "        file.write(\"    if(InterpCounter==\"+str(i)+\") { snprintf(gf_name,100,\\\"\"+gf_interp_list[i].gf_description+\"\\\"); return; }\\n\")\n",
    "    file.write(\"    printf(\\\"Error. InterpCounter = %d unsupported. I should not be here.\\\\n\\\",InterpCounter); exit(1);\\n\")\n",
    "    file.write(\"}\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy_interp_counter'></a>\n",
    "\n",
    "## Step 2.d: C Code for Initializing and incrementing `InterpCounter` \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy_interp_counter}$$\n",
    "\n",
    "The gridfunctions are interpolated one at a time based on the current value of the index quantity `InterpCounter`. Here we write the C code needed for initializing and incrementing this variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.575270Z",
     "iopub.status.busy": "2021-03-07T17:07:27.574405Z",
     "iopub.status.idle": "2021-03-07T17:07:27.577526Z",
     "shell.execute_reply": "2021-03-07T17:07:27.578152Z"
    }
   },
   "outputs": [],
   "source": [
    "with open(os.path.join(Ccodesdir,\"src\",\"define_NumInterpFunctions.h\"), \"w\") as file:\n",
    "    file.write(\"#define NumInterpFunctions (\"+str(which_InterpCounter)+\"-4*(1-enable_nuc_eos))\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.585137Z",
     "iopub.status.busy": "2021-03-07T17:07:27.583938Z",
     "iopub.status.idle": "2021-03-07T17:07:27.589009Z",
     "shell.execute_reply": "2021-03-07T17:07:27.588248Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/src/interp_counter.cc\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/interp_counter.cc\n",
    "#include <assert.h>\n",
    "#include <stdio.h>\n",
    "#include <stdlib.h>\n",
    "#include <string.h>\n",
    "#include <math.h>\n",
    "#include <ctype.h>\n",
    "#include \"cctk.h\"\n",
    "#include \"cctk_Arguments.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "\n",
    "#include \"define_NumInterpFunctions.h\"\n",
    "\n",
    "void ArbGrid_InitializeInterpCounterToZero(CCTK_ARGUMENTS)\n",
    "{\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  *InterpCounter = 0;\n",
    "\n",
    "  if(verbose==2) printf(\"interp_arbgrid_MO_ETK: Just set InterpCounter to %d\\n\",*InterpCounter);\n",
    "}\n",
    "\n",
    "void ArbGrid_InitializeInterpCounter(CCTK_ARGUMENTS)\n",
    "{\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  if(cctk_iteration == interp_out_iteration) {\n",
    "    *InterpCounter = 1;\n",
    "    if(verbose==2) printf(\"interp_arbgrid_MO_ETK: Just set InterpCounter to %d ; ready to start looping over interpolated gridfunctions!\\n\",\n",
    "                          *InterpCounter);\n",
    "  }\n",
    "}\n",
    "\n",
    "// This function increments InterpCounter if we are at the interp_out_iteration until\n",
    "// it hits NumInterpFunctions. At this iteration, InterpCounter is set to zero, which\n",
    "// exits the loop.\n",
    "void ArbGrid_IncrementInterpCounter(CCTK_ARGUMENTS)\n",
    "{\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  if(*InterpCounter == NumInterpFunctions-1) {\n",
    "    *InterpCounter = 0;\n",
    "    if(verbose==2) printf(\"interp_arbgrid_MO_ETK: Finished! Just zeroed InterpCounter.\\n\");\n",
    "  } else {\n",
    "    (*InterpCounter)++;\n",
    "    if(verbose==2) printf(\"interp_arbgrid_MO_ETK: Just incremented InterpCounter to %d of %d\\n\",*InterpCounter,NumInterpFunctions-1);\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='validationagainstfm'></a>\n",
    "\n",
    "# Step 2.e: Validation of interpolated data against exact Fishbone-Moncrief data \\[Back to [top](#toc)\\]\n",
    "$$\\label{validationagainstfm}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.597921Z",
     "iopub.status.busy": "2021-03-07T17:07:27.597193Z",
     "iopub.status.idle": "2021-03-07T17:07:27.600155Z",
     "shell.execute_reply": "2021-03-07T17:07:27.599621Z"
    }
   },
   "outputs": [],
   "source": [
    "# # Step 1c: Call the FishboneMoncriefID() function from within the\n",
    "# #          FishboneMoncriefID/FishboneMoncriefID.py module.\n",
    "# import FishboneMoncriefID.FishboneMoncriefID as fmid\n",
    "\n",
    "# old_glb_gridfcs_list = gri.glb_gridfcs_list\n",
    "# # Step 1: Set up the Fishbone-Moncrief initial data. This sets all the ID gridfunctions.\n",
    "# gri.glb_gridfcs_list = [] # Reset list of gridfunctions\n",
    "# fmid.FishboneMoncriefID(\"Spherical\")\n",
    "# gammaDD = ixp.zerorank2()\n",
    "\n",
    "# DIM = 3\n",
    "# for i in range(DIM):\n",
    "#     for j in range(DIM):\n",
    "#         if i<=j:\n",
    "#             gammaDD[i][j] = fmid.IDgammaDD[i][j]\n",
    "#         else:\n",
    "#             gammaDD[i][j] = fmid.IDgammaDD[j][i]\n",
    "\n",
    "# # gamma_{ij} v^i_{(n)} v^j_{(n)}\n",
    "# Gammacontraction = sp.sympify(0)\n",
    "# for i in range(DIM):\n",
    "#     for j in range(DIM):\n",
    "#         Gammacontraction += gammaDD[i][j] * fmid.IDValencia3velocityU[i] * fmid.IDValencia3velocityU[j]\n",
    "\n",
    "# Gammafactor = sp.sqrt(1 / (1 - Gammacontraction))\n",
    "\n",
    "# # -={ F-M quantities: Generate C code from expressions and output to file }=-\n",
    "# FishboneMoncrief_to_print = [\\\n",
    "#                      lhrh(lhs=\"Gammafactor\",rhs=Gammafactor),\\\n",
    "#                      lhrh(lhs=\"Gamma_times_ValenciavU0\",rhs=Gammafactor*fmid.IDValencia3velocityU[0]),\\\n",
    "#                      lhrh(lhs=\"Gamma_times_ValenciavU1\",rhs=Gammafactor*fmid.IDValencia3velocityU[1]),\\\n",
    "#                      lhrh(lhs=\"Gamma_times_ValenciavU2\",rhs=Gammafactor*fmid.IDValencia3velocityU[2]),\\\n",
    "#                      ]\n",
    "# fin.FD_outputC(os.path.join(Ccodesdir,\"src\",\"FM_Gamma__Gamma_times_Valenciavs_sphbasis.h\"),FishboneMoncrief_to_print,\n",
    "#                params=\"outCverbose=False,CSE_enable=True\")\n",
    "\n",
    "# # Restore old gridfunctions list:\n",
    "# gri.glb_gridfcs_list = old_glb_gridfcs_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.607861Z",
     "iopub.status.busy": "2021-03-07T17:07:27.606893Z",
     "iopub.status.idle": "2021-03-07T17:07:27.610411Z",
     "shell.execute_reply": "2021-03-07T17:07:27.610882Z"
    }
   },
   "outputs": [],
   "source": [
    "# %%writefile $Ccodesdir/src/FM_validation.cc\n",
    "\n",
    "# #include <assert.h>\n",
    "# #include <stdio.h>\n",
    "# #include <stdlib.h>\n",
    "# #include <string.h>\n",
    "# #include <math.h>\n",
    "# #include <ctype.h>\n",
    "# // Needed for dealing with Cactus/ETK infrastructure\n",
    "# #include \"cctk.h\"\n",
    "# #include \"cctk_Arguments.h\"\n",
    "# #include \"cctk_Parameters.h\"\n",
    "# // Needed for low-level interpolation functions\n",
    "# #include \"util_Table.h\"\n",
    "# #include \"util_String.h\"\n",
    "\n",
    "# // C++ function prototypes:\n",
    "# extern void Interpolate_to_dest_grid(const cGH *cctkGH,const CCTK_INT interp_num_points, const CCTK_INT interp_order,\n",
    "#                                      const CCTK_REAL *point_x_temp,const CCTK_REAL *point_y_temp,const CCTK_REAL *point_z_temp,\n",
    "#                                      const CCTK_STRING input_array_names[1], CCTK_REAL *output_f[1]);\n",
    "# extern void get_gf_name(const int InterpCounter,char gf_name[100]);\n",
    "\n",
    "# #define FREE_2D_GENERIC(type,array,ni,nj) for(int cc = 0; cc < ni;cc++) free((void *)array[cc]); \\\n",
    "# /**/                                                                    free((void *)array);\n",
    "\n",
    "# void FM_validation(CCTK_ARGUMENTS)\n",
    "# {\n",
    "#     DECLARE_CCTK_ARGUMENTS;\n",
    "#     DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "#     const CCTK_INT sph_Nr  = 3200;\n",
    "#     const CCTK_INT sph_Nth = 1;\n",
    "#     const CCTK_INT sph_Nph = 160;\n",
    "#     const CCTK_REAL sph_rmin = 0.1;\n",
    "#     const CCTK_REAL sph_rmax = 50.0;\n",
    "#     const CCTK_REAL sph_thmin = M_PI/2.0;\n",
    "#     const CCTK_REAL sph_thmax = M_PI/2.0;\n",
    "#     const CCTK_REAL sph_phmin = 0;\n",
    "#     const CCTK_REAL sph_phmax = 2.0*M_PI;\n",
    "\n",
    "#     const CCTK_INT num_interp_points = sph_Nr*sph_Nth*sph_Nph;\n",
    "\n",
    "#   // STEP 1: IF GAMMA*VALENCIA,PROCEED. IF RHO_B, OUTPUT FM FOR VEL DATA. ELSE RETURN.\n",
    "#   // Perform interpolation only at iteration == interp_out_iteration:\n",
    "#   if(cctk_iteration != interp_out_iteration) return;\n",
    "#   char gf_name[100]; get_gf_name(*InterpCounter,gf_name);\n",
    "\n",
    "# //  if(strncmp(gf_name,\"Lorentz factor, times Valencia\",30) == 0) {\n",
    "#   if(0 == 0) {\n",
    "#       // Perform interpolation!\n",
    "#       // Process zero (CCTK_MyProc(cctkGH)==0) is responsible for directing the interpolation.\n",
    "#       //    All other processes must see the cctk_InterpGridArrays() within Interpolate_to_dest_grid(),\n",
    "#       //    so that the MPI calls work properly, but these nonzero processes can call\n",
    "#       //    Interpolate_to_dest_grid() with number of interpolated points set to zero, and\n",
    "#       //    without needing a malloc().\n",
    "\n",
    "#       if(CCTK_MyProc(cctkGH)==0) {\n",
    "#         CCTK_REAL *points_x,*points_y,*points_z,**output_f;\n",
    "#         // The name of the input gridfunction is always \"interp_arbgrid_MO_ETK::interped_gf\":\n",
    "#         const CCTK_STRING input_array_names[1] = { \"interp_arbgrid_MO_ETK::interped_gf\" };\n",
    "\n",
    "#         // STEP 1: Construct list of desired interpolation destination points:\n",
    "#         points_x = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*num_interp_points);\n",
    "#         points_y = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*num_interp_points);\n",
    "#         points_z = (CCTK_REAL *)malloc(sizeof(CCTK_REAL)*num_interp_points);\n",
    "#         output_f = (CCTK_REAL **)malloc(1 * sizeof(CCTK_REAL *));\n",
    "#         for(int cc = 0; cc < 1; cc++) output_f[cc]=(CCTK_REAL *)malloc(num_interp_points * sizeof(CCTK_REAL));\n",
    "\n",
    "#         // STEP 2: ALLOCATE INTERPOLATION ARRAYS SET INTERPOLATION POINT ARRAYS\n",
    "#         CCTK_INT pointcount = 0;\n",
    "#         for(int ir=0;ir<sph_Nr;ir++) for(int ith=0;ith<sph_Nth;ith++) for(int iph=0;iph<sph_Nph;iph++) {\n",
    "#             const CCTK_REAL r  = sph_rmin  + (CCTK_REAL)ir /((CCTK_REAL)sph_Nr ) * (sph_rmax - sph_rmin);\n",
    "#             const CCTK_REAL th = sph_thmin + (CCTK_REAL)ith/((CCTK_REAL)sph_Nth) * (sph_thmax - sph_thmin);\n",
    "#             const CCTK_REAL ph = sph_phmin + (CCTK_REAL)iph/((CCTK_REAL)sph_Nph) * (sph_phmax - sph_phmin);\n",
    "#             points_x[pointcount] = r*sin(th)*cos(ph);\n",
    "#             points_y[pointcount] = r*sin(th)*sin(ph);\n",
    "#             points_z[pointcount] = r*cos(th);\n",
    "#             pointcount++;\n",
    "#         } // END for(int ir=0;ir<sph_Nr;ir++) for...\n",
    "\n",
    "#         // STEP 3: Looping over interp order as desired, interpolate to destination points & output to file\n",
    "#         for(int order=1;order<=4;order*=2) {\n",
    "#             printf(\"ASCII FM Validation: %d pts; Interpolating\\033[1m %s \\033[0m... using interpolation order = %d\\n\",num_interp_points,gf_name,order);\n",
    "#       //Interpolate_to_dest_grid(cGH *cctkGH,CCTK_INT interp_num_points, CCTK_INT interp_order,\n",
    "#       //                        CCTK_REAL *point_x_temp,CCTK_REAL *point_y_temp,CCTK_REAL *point_z_temp,\n",
    "#       //                        const CCTK_STRING input_array_names[1], CCTK_REAL *output_f[1])\n",
    "#             Interpolate_to_dest_grid(cctkGH, num_interp_points, order,\n",
    "#                                      points_x,points_y,points_z, input_array_names, output_f);\n",
    "\n",
    "#             // Step 1.d.i: Sanity check -- check for bad point:\n",
    "#         #pragma omp parallel for\n",
    "#             for(int i=0;i<num_interp_points;i++) {\n",
    "#                 if(output_f[0][i] > 1e20) {\n",
    "#                     printf(\"BAD POINT: %s %d %e %e %e %e\\n\",gf_name,i,points_x[i],points_y[i],points_z[i], output_f[0][i]);\n",
    "#                     exit(1);\n",
    "#                 } // END if(output_f[0][i] > 1e20)\n",
    "#             } // END for(int i=0;i<num_interp_points;i++)\n",
    "#             char filename[500];\n",
    "#             sprintf (filename, \"%s/validation_points-%s-order%d.asc\", out_dir,gf_name,order);\n",
    "#             FILE *file;\n",
    "#             file = fopen (filename,\"w\");\n",
    "#             printf(\"WRITING to file %s\\n\",filename);\n",
    "#             if (! file) {\n",
    "#                 CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,\n",
    "#                             \"interp_dest_grid__ET_thorn: Cannot open ASCII output file '%s'\", filename);\n",
    "#                 exit(1);\n",
    "#             } // END if (! file)\n",
    "#             pointcount = 0;\n",
    "#             for(int ir=0;ir<sph_Nr;ir++) for(int ith=0;ith<sph_Nth;ith++) for(int iph=0;iph<sph_Nph;iph++) {\n",
    "#                 const CCTK_REAL xx = points_x[pointcount];\n",
    "#                 const CCTK_REAL yy = points_y[pointcount];\n",
    "#                 const CCTK_REAL zz = points_z[pointcount];\n",
    "#                 fprintf(file,\"%e %e %e %e\\n\",xx,yy,zz,output_f[0][pointcount]);\n",
    "#                 pointcount++;\n",
    "#             }\n",
    "#             fclose(file);\n",
    "#         } // END for(int order=1;order<=4;order*=2)\n",
    "\n",
    "#         // STEP 3: FREE THE MALLOCs for destination grids and interpolation output\n",
    "#         free(points_x);\n",
    "#         free(points_y);\n",
    "#         free(points_z);\n",
    "#         FREE_2D_GENERIC(CCTK_REAL,output_f,1,num_interp_points);\n",
    "\n",
    "#       } else if(CCTK_MyProc(cctkGH)!=0) {\n",
    "#         // On all MPI processes that are nonzero, only call the interpolation function\n",
    "\n",
    "#         CCTK_REAL *points_x,*points_y,*points_z,**output_f;\n",
    "#         // The name of the input gridfunction is always \"interp_arbgrid_MO_ETK::interped_gf\":\n",
    "#         const CCTK_STRING input_array_names[1] = { \"interp_arbgrid_MO_ETK::interped_gf\" };\n",
    "\n",
    "#         //    to ensure the MPI calls from the actual interpolation (driven by proc==0) are seen.\n",
    "#         for(int order=1;order<=4;order*=2) {\n",
    "#             Interpolate_to_dest_grid(cctkGH, 0, order,points_x,points_y,points_z, input_array_names, output_f);\n",
    "#         } // END for(int order=1;order<=4;order*=2)\n",
    "#       } // END if(CCTK_MyProc(cctkGH)...)\n",
    "\n",
    "#   }\n",
    "#     if(strncmp(gf_name,\"IGM density primitive\",21) == 0 && CCTK_MyProc(cctkGH)==0) {\n",
    "#       char filename[500];\n",
    "#       sprintf (filename, \"%s/FMvalidation_points.asc\", out_dir);\n",
    "#       FILE *file;\n",
    "#       file = fopen (filename,\"w\");\n",
    "#       printf(\"WRITING to file %s\\n\",filename);\n",
    "\n",
    "#       for(int ir=0;ir<sph_Nr;ir++) for(int ith=0;ith<sph_Nth;ith++) for(int iph=0;iph<sph_Nph;iph++) {\n",
    "#         const CCTK_REAL xx0 = sph_rmin  + (CCTK_REAL)ir /((CCTK_REAL)sph_Nr ) * (sph_rmax - sph_rmin);\n",
    "#         const CCTK_REAL xx1 = sph_thmin + (CCTK_REAL)ith/((CCTK_REAL)sph_Nth) * (sph_thmax - sph_thmin);\n",
    "#         const CCTK_REAL xx2 = sph_phmin + (CCTK_REAL)iph/((CCTK_REAL)sph_Nph) * (sph_phmax - sph_phmin);\n",
    "#         const CCTK_REAL xx = xx0*sin(xx1)*cos(xx2);\n",
    "#         const CCTK_REAL yy = xx0*sin(xx1)*sin(xx2);\n",
    "#         const CCTK_REAL zz = xx0*cos(xx1);\n",
    "#         if(xx0 < r_in) {\n",
    "#             fprintf(file,\"%e %e %e %e %e %e %e\\n\",xx,yy,zz,\n",
    "#                     0.0,0.0,0.0,1.0);\n",
    "#         } else {\n",
    "#             CCTK_REAL Gammafactor,Gamma_times_ValenciavU0,Gamma_times_ValenciavU1,Gamma_times_ValenciavU2;\n",
    "#     #include \"FM_Gamma__Gamma_times_Valenciavs_sphbasis.h\"\n",
    "#             fprintf(file,\"%e %e %e %e %e %e %e\\n\",xx,yy,zz,\n",
    "#                     Gamma_times_ValenciavU0,Gamma_times_ValenciavU1,Gamma_times_ValenciavU2,Gammafactor);\n",
    "#         } // END if(xx0 < r_in)\n",
    "#       } // END for(int ir=0;ir<sph_Nr;ir++) for...\n",
    "#       fclose(file);\n",
    "#       return;\n",
    "#   } // END if(strncmp(gf_name,\"Lorentz factor, times Valencia\",30) == 0)\n",
    "# } // END function\n",
    "# #undef FREE_2D_GENERIC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cclfiles'></a>\n",
    "\n",
    "# Step 3: Define how this module interacts and interfaces with the larger Einstein Toolkit infrastructure \\[Back to [top](#toc)\\]\n",
    "$$\\label{cclfiles}$$\n",
    "\n",
    "Writing a module (\"thorn\") within the Einstein Toolkit requires that three \"ccl\" files be constructed, all in the root directory of the thorn:\n",
    "\n",
    "1. `interface.ccl`: defines the gridfunction groups needed, and provides keywords denoting what this thorn provides and what it should inherit from other thorns.\n",
    "1. `param.ccl`: specifies free parameters within the thorn.\n",
    "1. `schedule.ccl`: allocates storage for gridfunctions, defines how the thorn's functions should be scheduled in a broader simulation, and specifies the regions of memory written to or read from gridfunctions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='makecodedefn'></a>\n",
    "\n",
    "## Step 3.a: `make.code.defn` \\[Back to [top](#toc)\\]\n",
    "$$\\label{makecodedefn}$$\n",
    "\n",
    "Before writing the \"ccl\" files, we first add Einstein Toolkit's equivalent of a Makefile, the `make.code.defn` file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.616070Z",
     "iopub.status.busy": "2021-03-07T17:07:27.615257Z",
     "iopub.status.idle": "2021-03-07T17:07:27.619057Z",
     "shell.execute_reply": "2021-03-07T17:07:27.618519Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/src/make.code.defn\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/make.code.defn\n",
    "# Main make.code.defn file for thorn interp_arbgrid_MO_ETK\n",
    "\n",
    "# Source files in this directory\n",
    "SRCS =  main_function.cc unstagger_A_fields.cc interp_counter.cc \\\n",
    "        construct_function_to_interpolate__store_to_interped_gf.cc # FM_validation.cc # <- For FishboneMoncriefID validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='interfaceccl'></a>\n",
    "\n",
    "## Step 3.b: `interface.ccl` \\[Back to [top](#toc)\\]\n",
    "$$\\label{interfaceccl}$$\n",
    "\n",
    "Let's now write `interface.ccl`. The [official Einstein Toolkit (Cactus) documentation](http://einsteintoolkit.org/usersguide/UsersGuide.html) defines what must/should be included in an `interface.ccl` file [**here**](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-178000D2.2). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.625492Z",
     "iopub.status.busy": "2021-03-07T17:07:27.624325Z",
     "iopub.status.idle": "2021-03-07T17:07:27.629892Z",
     "shell.execute_reply": "2021-03-07T17:07:27.629255Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/interface.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/interface.ccl\n",
    "\n",
    "# With \"implements\", we give our thorn its unique name.\n",
    "implements: interp_arbgrid_MO_ETK\n",
    "\n",
    "# By \"inheriting\" other thorns, we tell the Toolkit that we\n",
    "#   will rely on variables/function that exist within those\n",
    "#   functions.\n",
    "inherits:   admbase IllinoisGRMHD Grid\n",
    "inherits:   WeylScal4 # Needed for Weyl scalars psi4, psi3, psi..., and Weyl invariants I & J.\n",
    "# For FM ID comparisons:\n",
    "# inherits:   FishboneMoncriefID\n",
    "\n",
    "# Tell the Toolkit that we want \"interped_gf\" and \"InterpCounter\"\n",
    "#    and invariants to NOT be visible to other thorns, by using\n",
    "#    the keyword \"private\". Note that declaring these\n",
    "#    gridfunctions here *does not* allocate memory for them;\n",
    "#    that is done by the schedule.ccl file.\n",
    "private:\n",
    "CCTK_REAL interpolation_gf type=GF timelevels=3 tags='Checkpoint=\"no\"'\n",
    "{\n",
    "  interped_gf\n",
    "} \"Gridfunction containing output from interpolation.\"\n",
    "\n",
    "CCTK_REAL unstaggered_A_fields type=GF timelevels=3 tags='Checkpoint=\"no\"'\n",
    "{\n",
    "  Ax_unstaggered,Ay_unstaggered,Az_unstaggered\n",
    "} \"Unstaggered A-field components.\"\n",
    "\n",
    "\n",
    "int InterpCounterVar type = SCALAR tags='checkpoint=\"no\"'\n",
    "{\n",
    "  InterpCounter\n",
    "} \"Counter that keeps track of which function we are interpolating.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='paramccl'></a>\n",
    "\n",
    "## Step 3.c: `param.ccl` \\[Back to [top](#toc)\\]\n",
    "$$\\label{paramccl}$$\n",
    "\n",
    "We will now write the file `param.ccl`. This file allows the listed parameters to be set at runtime. We also give allowed ranges and default values for each parameter. More information on this file's syntax can be found in the [official Einstein Toolkit documentation](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-183000D2.3). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.635836Z",
     "iopub.status.busy": "2021-03-07T17:07:27.635008Z",
     "iopub.status.idle": "2021-03-07T17:07:27.639947Z",
     "shell.execute_reply": "2021-03-07T17:07:27.638955Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/param.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/param.ccl\n",
    "\n",
    "# Output the interpolated data to the IO::out_dir directory:\n",
    "shares: IO\n",
    "USES STRING out_dir\n",
    "\n",
    "# These parameters are used to output Hybrid EOS information\n",
    "shares: IllinoisGRMHD\n",
    "USES CCTK_INT neos\n",
    "USES CCTK_REAL Gamma_th\n",
    "\n",
    "# These parameters are used to output nuclear (tabulated) EOS information\n",
    "shares: EOS_Omni\n",
    "USES CCTK_STRING nuceos_table_name\n",
    "\n",
    "# For FM ID comparisons:\n",
    "# shares: FishboneMoncriefID\n",
    "# USES KEYWORD M\n",
    "# USES KEYWORD a\n",
    "# USES KEYWORD r_in\n",
    "# USES KEYWORD r_at_max_density\n",
    "\n",
    "restricted:\n",
    "\n",
    "########################################\n",
    "# BASIC THORN STEERING PARAMETERS\n",
    "CCTK_INT interp_out_iteration \"Which iteration to interpolate to destination grids?\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 960000\n",
    "\n",
    "## Interpolator information\n",
    "CCTK_INT verbose \"Set verbosity level: 1=useful info; 2=moderately annoying (though useful for debugging)\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:2 :: \"0 = no output; 1=useful info; 2=moderately annoying (though useful for debugging)\"\n",
    "} 2\n",
    "\n",
    "CCTK_INT A_fields_are_staggered \"Are A fields staggered? 1 = yes; 0 = no. Default to yes.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:1 :: \"\"\n",
    "} 1\n",
    "\n",
    "##########\n",
    "# Cartesian position of center of output grid (usually center of BH).\n",
    "CCTK_REAL x_center \"x-position of center.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  *:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "CCTK_REAL y_center \"y-position of center.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  *:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "CCTK_REAL z_center \"z-position of center.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  *:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "##########\n",
    "# Shift offset:\n",
    "CCTK_REAL beta_offsetU0 \"Offset to betax, to account for coordinate drift in x direction.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  *:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "CCTK_REAL beta_offsetU1 \"Offset to betay, to account for coordinate drift in y direction.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  *:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "CCTK_REAL beta_offsetU2 \"Offset to betaz, to account for coordinate drift in z direction.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  *:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "CCTK_REAL out_of_bounds_interp_xyz \"Do not interpolate points with fabs(xyz) > out_of_bounds_interp_xyz, where xyz are centered at xyz_center (usually center of BH). Fill dataset with NaN instead.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"Any positive number\"\n",
    "} 1E100\n",
    "\n",
    "###########\n",
    "# Nuclear EOS options\n",
    "CCTK_INT enable_nuc_eos \"Use nuclear (tabulated/microphysics) equation of state? 1 = yes; 0 = no. Default to no.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:1 :: \"\"\n",
    "} 0\n",
    "\n",
    "CCTK_STRING temperatureVarString \"Temperature GF name. Defaults to IllinoisGRMHD's.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  \"IllinoisGRMHD::igm_temperature\" :: \"IllinoisGRMHD's temperature gridfunction name\"\n",
    "  \"HydroBase::temperature\"         :: \"HydroBase's temperature gridfunction name\"\n",
    "  \".+\"                             :: \"Or use you can use your own thorn's temperature gridfunction name\"\n",
    "} \"IllinoisGRMHD::igm_temperature\"\n",
    "\n",
    "CCTK_STRING Y_eVarString \"Electron fraction GF name. Defaults to IllinoisGRMHD's.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  \"IllinoisGRMHD::igm_Ye\" :: \"IllinoisGRMHD's electron fraction gridfunction name\"\n",
    "  \"HydroBase::Y_e\"        :: \"HydroBase's  electron fraction gridfunction name\"\n",
    "  \".+\"                    :: \"Or use you can use your own thorn's  electron fraction gridfunction name\"\n",
    "} \"IllinoisGRMHD::igm_Ye\"\n",
    "\n",
    "CCTK_STRING epsVarString \"Specific internal energy GF name. Defaults to IllinoisGRMHD's.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  \"IllinoisGRMHD::igm_eps\" :: \"IllinoisGRMHD's specific internal energy gridfunction name\"\n",
    "  \"HydroBase::eps\"         :: \"HydroBase's specific internal energy gridfunction name\"\n",
    "  \".+\"                     :: \"Or use you can use your own thorn's specific internal energy gridfunction name\"\n",
    "} \"IllinoisGRMHD::igm_eps\"\n",
    "\n",
    "CCTK_STRING entropyVarString \"Entropy GF name. Defaults to IllinoisGRMHD's.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  \"IllinoisGRMHD::igm_entropy\" :: \"IllinoisGRMHD's entropy gridfunction name\"\n",
    "  \"HydroBase::entropy\"         :: \"HydroBase's entropy gridfunction name\"\n",
    "  \".+\"                         :: \"Or use you can use your own thorn's entropy gridfunction name\"\n",
    "} \"IllinoisGRMHD::igm_entropy\"\n",
    "\n",
    "CCTK_STRING HVarString \"Hamiltonian constraint GF name. Defaults to Baika's.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  \"Baikal::HGF\" :: \"Baikal's Hamiltonian constraint gridfunction name\"\n",
    "  \"ML_BSSN::H\"  :: \"ML_BSSN's Hamiltonian constraint gridfunction name\"\n",
    "  \".+\"                         :: \"Or use you can use your own thorn's entropy gridfunction name\"\n",
    "} \"Baikal::HGF\"\n",
    "\n",
    "CCTK_STRING M0VarString \"Momentum constraint (x-direction) Defaults to Baikal's.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  \"Baikal::MU0GF\" :: \"Baikal's Momentum constraint (x-direction) gridfunction name\"\n",
    "  \"ML_BSSN::M1\"   :: \"ML_BSSN's Momentum constraint (x-direction) gridfunction name\"\n",
    "  \".+\"                         :: \"Or use you can use your own thorn's entropy gridfunction name\"\n",
    "} \"Baikal::MU0GF\"\n",
    "\n",
    "CCTK_STRING M1VarString \"Momentum constraint (y-direction) GF name. Defaults to Baikal's.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  \"Baikal::MU1GF\" :: \"Baikal's Momentum constraint (y-direction) gridfunction name\"\n",
    "  \"ML_BSSN::M2\"   :: \"ML_BSSN's Momentum constraint (y-direction) gridfunction name\"\n",
    "  \".+\"                         :: \"Or use you can use your own thorn's entropy gridfunction name\"\n",
    "} \"Baikal::MU1GF\"\n",
    "\n",
    "CCTK_STRING M2VarString \"Momentum constraint (z-direction) GF name. Defaults to Baikal's.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  \"Baikal::MU2GF\" :: \"Baikal's Momentum constraint (z-direction) gridfunction name\"\n",
    "  \"ML_BSSN::M3\"   :: \"ML_BSSN's Momentum constraint (z-direction) gridfunction name\"\n",
    "  \".+\"                         :: \"Or use you can use your own thorn's entropy gridfunction name\"\n",
    "} \"Baikal::MU2GF\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='scheduleccl'></a>\n",
    "\n",
    "## Step 3.d: `schedule.ccl` \\[Back to [top](#toc)\\]\n",
    "$$\\label{scheduleccl}$$\n",
    "\n",
    "Finally, we will write the file `schedule.ccl`; its official documentation is found [here](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-186000D2.4). \n",
    "\n",
    "This file declares storage for variables declared in the `interface.ccl` file and specifies when the various parts of the thorn will be run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.647153Z",
     "iopub.status.busy": "2021-03-07T17:07:27.646158Z",
     "iopub.status.idle": "2021-03-07T17:07:27.650736Z",
     "shell.execute_reply": "2021-03-07T17:07:27.650034Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_arbgrid_MO_ETK/schedule.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/schedule.ccl\n",
    "\n",
    "STORAGE: interpolation_gf[3]\n",
    "STORAGE: unstaggered_A_fields[3]\n",
    "STORAGE: InterpCounterVar\n",
    "# STORAGE: interp_pointcoords_and_output_arrays\n",
    "\n",
    "#############################\n",
    "SCHEDULE ArbGrid_InitializeInterpCounterToZero AT CCTK_INITIAL\n",
    "{\n",
    "  LANG: C\n",
    "  OPTIONS: GLOBAL\n",
    "} \"Initialize InterpCounter variable to zero\"\n",
    "\n",
    "SCHEDULE ArbGrid_InitializeInterpCounterToZero AT CCTK_POST_RECOVER_VARIABLES\n",
    "{\n",
    "  LANG: C\n",
    "  OPTIONS: GLOBAL\n",
    "} \"Initialize InterpCounter variable to zero\"\n",
    "\n",
    "SCHEDULE ArbGrid_InitializeInterpCounter before ArbGrid_InterpGroup AT CCTK_ANALYSIS\n",
    "{\n",
    "  LANG: C\n",
    "  OPTIONS: GLOBAL\n",
    "} \"Initialize InterpCounter variable\"\n",
    "##################\n",
    "\n",
    "SCHEDULE GROUP ArbGrid_InterpGroup AT CCTK_ANALYSIS BEFORE CarpetLib_printtimestats BEFORE CarpetLib_printmemstats AFTER Convert_to_HydroBase WHILE interp_arbgrid_MO_ETK::InterpCounter\n",
    "{\n",
    "} \"Perform all interpolations. This group is only actually scheduled at cctk_iteration==interp_out_iteration.\"\n",
    "\n",
    "SCHEDULE unstagger_A_fields in ArbGrid_InterpGroup before construct_function_to_interpolate__store_to_interped_gf\n",
    "{\n",
    "  STORAGE: unstaggered_A_fields[3]\n",
    "  OPTIONS: GLOBAL,LOOP-LOCAL\n",
    "  SYNC:    unstaggered_A_fields\n",
    "  LANG: C\n",
    "} \"Unstagger A fields.\"\n",
    "\n",
    "SCHEDULE construct_function_to_interpolate__store_to_interped_gf in ArbGrid_InterpGroup before DoSum\n",
    "{\n",
    "  STORAGE: interpolation_gf[3],InterpCounterVar\n",
    "  OPTIONS: GLOBAL,LOOP-LOCAL\n",
    "  SYNC:    interpolation_gf\n",
    "  LANG: C\n",
    "} \"Construct the function to interpolate\"\n",
    "\n",
    "SCHEDULE Interpolate_to_dest_grid_main_function in ArbGrid_InterpGroup after construct_function_to_interpolate__store_to_interped_gf\n",
    "{\n",
    "  OPTIONS: GLOBAL\n",
    "  LANG: C\n",
    "} \"Perform interpolation and output result to file.\"\n",
    "\n",
    "# For FishboneMoncriefID validation only.\n",
    "# SCHEDULE FM_validation in ArbGrid_InterpGroup after Interpolate_to_dest_grid_main_function\n",
    "# {\n",
    "#   OPTIONS: GLOBAL\n",
    "#   LANG: C\n",
    "# } \"Perform interpolation and output result to 2D ASCII file.\"\n",
    "\n",
    "#######\n",
    "SCHEDULE ArbGrid_IncrementInterpCounter in ArbGrid_InterpGroup after Interpolate_to_dest_grid_main_function\n",
    "{\n",
    "  LANG: C\n",
    "  OPTIONS: GLOBAL\n",
    "} \"Increment InterpCounter variable, or set to zero once loop is complete.\"\n",
    "##################"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='readingoutputfile'></a>\n",
    "\n",
    "# Step 4: Python Script for Reading the Output File \\[Back to [top](#toc)\\]\n",
    "$$\\label{readingoutputfile}$$\n",
    "\n",
    "Here is a Python code for reading the output file generated by this thorn. It is based on a collection of Python scripts written by Bernard Kelly, available [here](https://bitbucket.org/zach_etienne/nrpy/src/master/mhd_diagnostics/). \n",
    "\n",
    "After generating the output file `interp_arbgrid_MO_ETK.dat` using the Einstein Toolkit thorn above, this script will read in all the data. Processing can then be done by straightforward modification of this script. Save the script as \"Interp_Arb_ReadIn.py\", and run it using the command\n",
    "\n",
    "**`python Interp_Arb_ReadIn.py interp_arbgrid_MO_ETK.dat 58 outfile`**\n",
    "\n",
    "Currently the last parameter \"outfile\" is required but not used.\n",
    "\n",
    "```python\n",
    "\"\"\"\n",
    "interp_arbgrid_MO_ETK.dat File Reader. Compatible with Python 2.7+ and 3.6+ at least.\n",
    "\n",
    "Zachariah B. Etienne\n",
    "\n",
    "Based on Python scripts written by Bernard Kelly:\n",
    "https://bitbucket.org/zach_etienne/nrpy/src/master/mhd_diagnostics/\n",
    "\n",
    "Find the latest version of this reader at the bottom of this Jupyter notebook:\n",
    "https://github.com/zachetienne/nrpytutorial/blob/master/Tutorial-ETK_thorn-Interpolation_to_Arbitrary_Grids_multi_order.ipynb\n",
    "\n",
    "Usage instructions:\n",
    "\n",
    "From the command-line, run via:\n",
    "python Interp_Arb_ReadIn.py interp_arbgrid_MO_ETK.dat [number of gridfunctions (58 or so)] [outfile]\n",
    "\n",
    "Currently the last parameter \"outfile\" is required but not actually used.\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "import struct\n",
    "import sys\n",
    "import argparse\n",
    "\n",
    "parser = argparse.ArgumentParser(description='Read file.')\n",
    "parser.add_argument(\"datafile\", help=\"main data file\")\n",
    "parser.add_argument(\"number_of_gridfunctions\", help=\"number of gridfunctions\")\n",
    "\n",
    "parser.add_argument(\"outfileroot\", help=\"root of output file names\")\n",
    "\n",
    "args = parser.parse_args()\n",
    "\n",
    "datafile = args.datafile\n",
    "outfileroot = args.outfileroot\n",
    "number_of_gridfunctions = int(args.number_of_gridfunctions)\n",
    "\n",
    "print(\"reading from \"+str(datafile))\n",
    "\n",
    "\"\"\"\n",
    "read_char_array():\n",
    "Reads a character array of size=\"size\"\n",
    "from a file (with file handle = \"filehandle\")\n",
    "and returns the character array as a proper \n",
    "Python string.\n",
    "\"\"\"\n",
    "def read_char_array(filehandle,size):\n",
    "    reached_end_of_string = False\n",
    "    chartmp = struct.unpack(str(size)+'s', filehandle.read(size))[0]\n",
    "\n",
    "    #https://docs.python.org/3/library/codecs.html#codecs.decode\n",
    "    char_array_orig = chartmp.decode('utf-8',errors='ignore')\n",
    "\n",
    "    char_array = \"\"\n",
    "    for i in range(len(char_array_orig)):\n",
    "        char = char_array_orig[i]\n",
    "        # C strings end in '\\0', which in Python-ese is '\\x00'.\n",
    "        #   As characters read after the end of the string will\n",
    "        #   generally be gibberish, we no inter append \n",
    "        #   to the output string after '\\0' is reached.\n",
    "        if   sys.version_info[0]==3 and bytes(char.encode('utf-8')) == b'\\x00':\n",
    "            reached_end_of_string = True\n",
    "        elif sys.version_info[0]==2 and char ==  '\\x00':\n",
    "            reached_end_of_string = True\n",
    "\n",
    "        if reached_end_of_string == False:\n",
    "            char_array += char\n",
    "        else:\n",
    "            pass # Continue until we've read 'size' bytes\n",
    "    return char_array\n",
    "\n",
    "\"\"\"\n",
    "read_header()\n",
    "Reads the header from a file.\n",
    "\"\"\"\n",
    "def read_header(filehandle):\n",
    "    # This function makes extensive use of Python's struct.unpack\n",
    "    # https://docs.python.org/3/library/struct.html\n",
    "    # First store gridfunction name and interpolation order used:\n",
    "    # fwrite(gf_name, 100*sizeof(char), 1, file);\n",
    "    gf_name = read_char_array(filehandle,100)\n",
    "    # fwrite(order, sizeof(CCTK_INT), 1, file);\n",
    "    order = struct.unpack('i',filehandle.read(4))[0]\n",
    "\n",
    "    # Then the number of interpolation points (stored as an int)\n",
    "    num_interp_points = struct.unpack('i',filehandle.read(4))[0]\n",
    "\n",
    "    magic_number_check = 1.130814081305130e-21\n",
    "    # fwrite( & magic_number, sizeof(CCTK_REAL), 1, file);\n",
    "    magic_number = struct.unpack('d', filehandle.read(8))[0]\n",
    "    if magic_number != magic_number_check:\n",
    "        print(\"Error: Possible file corruption: Magic number mismatch. Found magic number = \"+str(magic_number)+\" . Expected \"+str(magic_number_check))\n",
    "        exit(1)\n",
    "    # fwrite( & cctk_iteration, sizeof(CCTK_INT), 1, file);\n",
    "    cctk_iteration = struct.unpack('i', filehandle.read(4))[0]\n",
    "    # fwrite( & cctk_time, sizeof(CCTK_REAL), 1, file);\n",
    "    cctk_time      = struct.unpack('d', filehandle.read(8))[0]\n",
    "\n",
    "    return gf_name,order,num_interp_points,cctk_iteration,cctk_time\n",
    "\n",
    "# Now open the file and read all the data\n",
    "with open(datafile,\"rb\") as f:\n",
    "    # Main loop over all gridfunctions\n",
    "    for i in range(number_of_gridfunctions):\n",
    "        # Data are output in chunks, one gridfunction at a time, with metadata\n",
    "        #    for each gridfunction stored at the top of each chunk\n",
    "        # First read in the metadata:\n",
    "        gf_name, order, num_interp_points, cctk_iteration, cctk_time = read_header(f)\n",
    "        print(\"\\nReading gridfunction \"+gf_name+\", stored at interp order = \"+str(order))\n",
    "        data_chunk_size = num_interp_points*8 # 8 bytes per double-precision number\n",
    "        # Next read in the full gridfunction data\n",
    "        bytechunk = f.read(data_chunk_size)\n",
    "        # Process the data using NumPy's frombuffer() function:\n",
    "        #   https://docs.scipy.org/doc/numpy/reference/generated/numpy.frombuffer.html\n",
    "        buffer_res = np.frombuffer(bytechunk)\n",
    "        # Reshape the data into a 3D NumPy array:\n",
    "        #   https://docs.scipy.org/doc/numpy/reference/generated/numpy.reshape.html\n",
    "        #this_data = buffer_res.reshape(N0,N1,N2)\n",
    "\n",
    "        # Sanity check: Output data at all points:\n",
    "        with open(\"output-gf\"+str(i)+\".txt\",\"w\") as file:\n",
    "            for ii in range(num_interp_points):\n",
    "                file.write(str(ii) + \"\\t\" + str(buffer_res[ii])+\"\\n\")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: 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-ETK_thorn-Interpolation_to_Arbitrary_Grids_multi_order.pdf](Tutorial-ETK_thorn-Interpolation_to_Arbitrary_Grids_multi_order.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": 36,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:27.657108Z",
     "iopub.status.busy": "2021-03-07T17:07:27.655896Z",
     "iopub.status.idle": "2021-03-07T17:07:33.957574Z",
     "shell.execute_reply": "2021-03-07T17:07:33.956580Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ETK_thorn-\n",
      "    Interpolation_to_Arbitrary_Grids_multi_order.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-ETK_thorn-\n",
      "    Interpolation_to_Arbitrary_Grids_multi_order.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ETK_thorn-Interpolation_to_Arbitrary_Grids_multi_order\")"
   ]
  }
 ],
 "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}