Tutorial-ETK_thorn-Interpolation_to_Spherical_Grids.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_sph_grids_ETK`: An Einstein Toolkit Module for Interpolation to Spherical Grids\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 spherical sampling.\n",
    "\n",
    "**Notebook Status:** <font color='red'><b> In progress </b></font>\n",
    "\n",
    "**Validation Notes:** This module has not yet undergone validation testing.\n",
    "\n",
    "## Introduction:\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 spherical 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 spherical 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",
    "This notebook is organized as follows\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 ETK Interpolation Function\n",
    "    1. [Step 1.b](#sphericalgridnotes): Setting up the Spherical Grids\n",
    "    1. [Step 1.c](#fileformat): Outputting to File\n",
    "    1. [Step 1.d](#maininterpolator) The Main Interpolation Driver Function\n",
    "1. [Step 2](#nrpy): Use NRPy+ C Output to Set All Output 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): GRMHD quantities (***IN PROGRESS***)\n",
    "        1. [Step 2.a.ii](#nrpy4metric): Compute all 10 components of the 4-metric $g_{\\mu\\nu}$\n",
    "        1. [Step 2.a.iii](#nrpy4christoffels): Compute all 40 4-Christoffels $\\Gamma^{\\mu}_{\\nu\\delta}$\n",
    "    1. [Step 2.b](#nrpy_c_callingfunction): C code calling function for the NRPy+ C output\n",
    "    1. [Step 2.c](#nrpygetgfinterporder): `get_gf_name()` and `get_interp_order()` functions\n",
    "    1. [Step 2.d](#nrpy_interp_counter): C Code for Initializing and incrementing `InterpCounter`\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:36.497543Z",
     "iopub.status.busy": "2021-03-07T17:07:36.496425Z",
     "iopub.status.idle": "2021-03-07T17:07:36.503810Z",
     "shell.execute_reply": "2021-03-07T17:07:36.504728Z"
    }
   },
   "outputs": [],
   "source": [
    "import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface\n",
    "import shutil, os, sys, time     # Standard Python modules for multiplatform OS-level functions, benchmarking\n",
    "\n",
    "# Create C code output directory:\n",
    "Ccodesdir = \"interp_sph_grids_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/\"))"
   ]
  },
  {
   "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_sph_grid()`**, which  to file. \n",
    "\n",
    "**`Interpolate_to_sph_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_sph_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:36.513538Z",
     "iopub.status.busy": "2021-03-07T17:07:36.512634Z",
     "iopub.status.idle": "2021-03-07T17:07:36.518224Z",
     "shell.execute_reply": "2021-03-07T17:07:36.518727Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_sph_grids_ETK/src/Interpolate_to_sph_grid.h\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/Interpolate_to_sph_grid.h\n",
    "\n",
    "void Interpolate_to_sph_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",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  CCTK_INT ierr;\n",
    "\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",
    "\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]\n",
    "    = { (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]\n",
    "    = { (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",
    "  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": "markdown",
   "metadata": {},
   "source": [
    "<a id='sphericalgridnotes'></a>\n",
    "\n",
    "## Step 1.b: Setting up the Spherical Grids \\[Back to [top](#toc)\\]\n",
    "$$\\label{sphericalgridnotes}$$\n",
    "\n",
    "+ By default, we set logarithmic radial coordinates: $r(x_{0,i}) = R_0 + e^{x_{0,i}}$, where\n",
    "\n",
    "    + $x_{0,i} = x_{0, \\mathrm{beg}} + \\left(i+\\frac{1}{2}\\right) \\Delta x_0$\n",
    "    + $x_{0, {\\mathrm{beg}}} = \\log\\left( R_{\\mathrm{in}} - R_0 \\right)$\n",
    "    + $\\Delta x_0 = \\frac{1}{N_0}\\log\\left(\\frac{R_\\mathrm{out} - R_0}{R_\\mathrm{in} - R_0}\\right)$\n",
    "\n",
    "\n",
    "+ As for the polar angle $\\theta$, there are two options:\n",
    "    + **Option 1**: \n",
    "  $$ \\theta(x_{1,j})  \\, = \\, \\theta_c \\, + \\, \\left( \\pi - 2 \\theta_c \\right) x_{1,j} \\, + \\, \\xi \\, \\sin\\left(2 \\pi x_{1,j} \\right),$$\n",
    "      where\n",
    "        + $x_{1,j} = x_{1, \\mathrm{beg}} + \\left(j+\\frac{1}{2}\\right) \\Delta x_1$\n",
    "        + $\\Delta x_1 = \\frac{1}{N_1}$\n",
    "\n",
    "    + **Option 2**: \n",
    "  $$ \\theta(x_{1,j}) = \\frac{\\pi}{2} \\left[  1  + \\left(1-\\xi \\right) \\left(2 x_{1,j} - 1 \\right) + \\left( \\xi - \\frac{2 \\theta_c}{\\pi} \\right) \\left( 2 x_{1,j} - 1 \\right)^n \\right],$$\n",
    "      where\n",
    "        + $n$ is odd\n",
    "        + $x_{1,j} = x_{1, \\mathrm{beg}} + \\left(j+\\frac{1}{2}\\right) \\Delta x_1$\n",
    "        + $\\Delta x_1 = \\frac{1}{N_1}$\n",
    "\n",
    "\n",
    "+ The azimuthal angle $\\phi$ is uniform, so that $\\phi(x_{2,k}) = x_{2,k}$:\n",
    "  \n",
    "    + $x_{2,k} \\in [0,2\\pi]$\n",
    "    + $x_{2,k} = x_{2, \\mathrm{beg}} + \\left(k+\\frac{1}{2}\\right)\\Delta x_{2}$\n",
    "    + $\\Delta x_{2} = \\frac{ 2 \\pi }{N_2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:36.525134Z",
     "iopub.status.busy": "2021-03-07T17:07:36.524011Z",
     "iopub.status.idle": "2021-03-07T17:07:36.528212Z",
     "shell.execute_reply": "2021-03-07T17:07:36.527529Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_sph_grids_ETK/src/Set_up_interp_points_on_sph_grid.h\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/Set_up_interp_points_on_sph_grid.h\n",
    "\n",
    "void sph_grid_Interpolate_many_pts__set_interp_pts(CCTK_ARGUMENTS) {\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  CCTK_REAL dx0 = log( (Rout - R0) / (Rin - R0) ) / ((CCTK_REAL)N0);\n",
    "  CCTK_REAL dx1 =      1.0 / ((CCTK_REAL)N1);\n",
    "  CCTK_REAL dx2 = 2.0*M_PI / ((CCTK_REAL)N2);\n",
    "  CCTK_REAL x0_beg = log( Rin - R0 );\n",
    "  CCTK_INT which_pt = 0;\n",
    "  for(CCTK_INT k=0;k<N2;k++) for(CCTK_INT j=0;j<N1;j++) for(CCTK_INT i=0;i<N0;i++) {\n",
    "    CCTK_REAL x0_i = x0_beg + ((CCTK_REAL)i + 0.5)*dx0;\n",
    "    CCTK_REAL rr = R0 + exp(x0_i);\n",
    "\n",
    "    CCTK_REAL x1_j = x1_beg + ((CCTK_REAL)j + 0.5)*dx1;\n",
    "    CCTK_REAL th = -1e300;\n",
    "    if(theta_option == 1) {\n",
    "       th = th_c + (M_PI - 2.0*th_c)*x1_j + xi*sin(2.0*M_PI*x1_j);\n",
    "    } else if (theta_option == 2) {\n",
    "       th = M_PI/2.0 * ( 1.0 + (1.0 - xi)*(2.0*x1_j - 1.0) + (xi - 2.0*th_c/M_PI)*pow(2.0*x1_j - 1.0 ,th_n) );\n",
    "    } else {\n",
    "       printf(\"Error: theta_option = %d NOT SUPPORTED.\",theta_option);\n",
    "       exit(1);\n",
    "    }\n",
    "\n",
    "    CCTK_REAL x2_k = x2_beg + ((CCTK_REAL)k + 0.5)*dx2;\n",
    "    CCTK_REAL ph = x2_k;\n",
    "\n",
    "    points_x[which_pt] = rr*sin(th)*cos(ph);\n",
    "    points_y[which_pt] = rr*sin(th)*sin(ph);\n",
    "    points_z[which_pt] = rr*cos(th);\n",
    "    which_pt++;\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='fileformat'></a>\n",
    "\n",
    "## Step 1.c: 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:36.539566Z",
     "iopub.status.busy": "2021-03-07T17:07:36.536346Z",
     "iopub.status.idle": "2021-03-07T17:07:36.541945Z",
     "shell.execute_reply": "2021-03-07T17:07:36.542453Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_sph_grids_ETK/src/output_to_file.h\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/output_to_file.h\n",
    "\n",
    "#include \"define_NumInterpFunctions.h\"\n",
    "\n",
    "void output_to_file(CCTK_ARGUMENTS,char gf_name[100],int *order,CCTK_REAL *output_f[1]) {\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  char filename[100];\n",
    "  sprintf (filename, \"%s/interp_sph_grids.dat\", out_dir);\n",
    "  FILE *file;\n",
    "  if(*InterpCounter == 1) {\n",
    "    file = fopen (filename,\"w\");\n",
    "  } else {\n",
    "    file = fopen (filename,\"a+\");\n",
    "  }\n",
    "  if (! file) {\n",
    "    CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,\n",
    "                \"interp_sph_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",
    "\n",
    "  fwrite(&N0, sizeof(CCTK_INT), 1, file);\n",
    "  fwrite(&R0, sizeof(CCTK_REAL), 1, file);\n",
    "  fwrite(&Rin, sizeof(CCTK_REAL), 1, file);\n",
    "  fwrite(&Rout, sizeof(CCTK_REAL), 1, file);\n",
    "\n",
    "  fwrite(&N1, sizeof(CCTK_INT), 1, file);\n",
    "  fwrite(&x1_beg, sizeof(CCTK_REAL), 1, file);\n",
    "  fwrite(&theta_option, sizeof(CCTK_INT), 1, file);\n",
    "  fwrite(&th_c, sizeof(CCTK_REAL), 1, file);\n",
    "  fwrite(&xi, sizeof(CCTK_REAL), 1, file);\n",
    "  fwrite(&th_n, sizeof(CCTK_INT), 1, file);\n",
    "\n",
    "  fwrite(&N2, sizeof(CCTK_INT), 1, file);\n",
    "  fwrite(&x2_beg, sizeof(CCTK_REAL), 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)*N0*N1*N2, 1, file);\n",
    "  }\n",
    "\n",
    "  fclose(file);\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='maininterpolator'></a>\n",
    "\n",
    "## Step 1.d: The Main Interpolation Driver Function \\[Back to [top](#toc)\\]\n",
    "$$\\label{maininterpolator}$$\n",
    "\n",
    "The **`Interpolate_to_sph_grid_main_function()`** function calls the above functions as follows:\n",
    "\n",
    "1. **`sph_grid_Interpolate_many_pts__set_interp_pts()`**: First set up the spherical grids\n",
    "1. **`Interpolate_to_sph_grid()`**: Output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:36.548255Z",
     "iopub.status.busy": "2021-03-07T17:07:36.547359Z",
     "iopub.status.idle": "2021-03-07T17:07:36.550400Z",
     "shell.execute_reply": "2021-03-07T17:07:36.550889Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_sph_grids_ETK/src/main_function.cc\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/main_function.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",
    "// 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 \"Set_up_interp_points_on_sph_grid.h\"\n",
    "#include \"Interpolate_to_sph_grid.h\"\n",
    "#include \"output_to_file.h\"\n",
    "#include \"get_interp_order.h\"\n",
    "#include \"get_gf_name.h\"\n",
    "\n",
    "void Interpolate_to_sph_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",
    "\n",
    "  // Set up spherically sampled interpolation grid arrays points_x,points_y,points_z:\n",
    "  sph_grid_Interpolate_many_pts__set_interp_pts(CCTK_PASS_CTOC);\n",
    "\n",
    "  // Set up output array:\n",
    "  CCTK_REAL *output_f[1];\n",
    "  output_f[0] = output_interped;\n",
    "  // The name of the input gridfunction is always \"interp_sph_grids_ETK::interped_gf\":\n",
    "  const CCTK_STRING input_array_names[1] = { \"interp_sph_grids_ETK::interped_gf\" };\n",
    "\n",
    "  // Perform interpolation!\n",
    "  int order = get_interp_order(*InterpCounter);\n",
    "  char gf_name[100];\n",
    "  get_gf_name(*InterpCounter,gf_name);\n",
    "  printf(\"Interpolating\\033[1m %s \\033[0m... using interpolation order = %d\\n\",gf_name,order);\n",
    "  Interpolate_to_sph_grid(cctkGH, N0*N1*N2, order,\n",
    "                             points_x,points_y,points_z, input_array_names, output_f);\n",
    "\n",
    "  if(CCTK_MyProc(cctkGH)==0) {\n",
    "    for(int i=0;i<N0*N1*N2;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(0);\n",
    "        }\n",
    "    }\n",
    "    output_to_file(CCTK_PASS_CTOC,gf_name,&order,output_f);\n",
    "    printf(\"Interpolate_to_sph_grid_main_function(): Just output to file at iteration %d\\n\",cctk_iteration);\n",
    "  } else {\n",
    "    printf(\"Interpolate_to_sph_grid_main_function(): Process !=0 waiting for file output at iteration %d\\n\",cctk_iteration);\n",
    "  }\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}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:36.561278Z",
     "iopub.status.busy": "2021-03-07T17:07:36.560170Z",
     "iopub.status.idle": "2021-03-07T17:07:36.908600Z",
     "shell.execute_reply": "2021-03-07T17:07:36.907995Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2 P0: Import needed NRPy+ modules/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', 'description order')\n",
    "gf_interp_list = []\n",
    "gf_interp_list.append(gf_interp(\"dummy -- used because this is a 1-offset array\",order=\"-100\"))\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",
    "        output_type=\"w\"\n",
    "\n",
    "    with open(filename, output_type) as file:\n",
    "        file.write(\"if(*InterpCounter == \"+str(which_InterpCounter)+\") {\\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\",\"\",\"\"],\"   \",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": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:36.913245Z",
     "iopub.status.busy": "2021-03-07T17:07:36.912418Z",
     "iopub.status.idle": "2021-03-07T17:07:36.914922Z",
     "shell.execute_reply": "2021-03-07T17:07:36.915446Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.a:  Set up NRPy-based \"list_of_functions_to_interpolate.h\"\n",
    "\n",
    "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 (*IN PROGRESS*) \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpygrmhd}$$\n",
    "\n",
    "These include\n",
    "\n",
    "* $\\rho_b$, the baryonic density (i.e., the HydroBase variable $\\verb|rho|$)\n",
    "* $P$, the total gas pressure (i.e., the HydroBase variable $\\verb|press|$)\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",
    "$$\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",
    "\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": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:36.928424Z",
     "iopub.status.busy": "2021-03-07T17:07:36.927510Z",
     "iopub.status.idle": "2021-03-07T17:07:36.950180Z",
     "shell.execute_reply": "2021-03-07T17:07:36.950664Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.a.i: GRMHD quantities (*IN PROGRESS*)\n",
    "# INPUT GRIDFUNCTIONS: The AUX or EVOL designation is *not* used in diagnostic modules.\n",
    "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUX\",\"gammaDD\", \"sym01\")\n",
    "betaU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"betaU\")\n",
    "alpha = gri.register_gridfunctions(\"AUX\",\"alpha\")\n",
    "\n",
    "DIM=3\n",
    "\n",
    "gf_interp_list.append(gf_interp(\"IGM density primitive\",order=\"1\"))\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\",order=\"1\"))\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",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:36.960945Z",
     "iopub.status.busy": "2021-03-07T17:07:36.959822Z",
     "iopub.status.idle": "2021-03-07T17:07:37.612516Z",
     "shell.execute_reply": "2021-03-07T17:07:37.611989Z"
    }
   },
   "outputs": [],
   "source": [
    "IGMvU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"IGMvU\")\n",
    "Valenciav = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    Valenciav[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]*Valenciav[i]*Valenciav[j]\n",
    "\n",
    "Gamma_times_ValenciavU = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    Gamma_times_ValenciavU[i] = sp.sqrt(1/(1 - v_dot_v))*Valenciav[i]\n",
    "    gf_interp_list.append(gf_interp(\"Lorentz factor, times Valencia vU\"+str(i),order=\"1\"))\n",
    "    interp_expr = Gamma_times_ValenciavU[i]\n",
    "    which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "# For testing:\n",
    "# gf_interp_list.append(gf_interp(\"Lorentz factor\",order=\"1\"))\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),order=\"1\"))\n",
    "#     interp_expr = Valenciav[i]\n",
    "#     which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)\n",
    "\n",
    "BU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"BU\")\n",
    "for i in range(DIM):\n",
    "    gf_interp_list.append(gf_interp(\"IGM magnetic field component B\"+str(i),order=\"1\"))\n",
    "    interp_expr = BU[i]\n",
    "    which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy4metric'></a>\n",
    "\n",
    "### Step 2.a.ii: 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": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:37.645821Z",
     "iopub.status.busy": "2021-03-07T17:07:37.635895Z",
     "iopub.status.idle": "2021-03-07T17:07:37.672082Z",
     "shell.execute_reply": "2021-03-07T17:07:37.672597Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.a.ii: Compute all 10 components of the 4-metric g_{\\mu\\nu}\n",
    "# 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),order=\"4\"))\n",
    "        interp_expr = g4DD[mu][nu]\n",
    "        which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy4christoffels'></a>\n",
    "\n",
    "### Step 2.a.iii: Compute all 40 4-Christoffels $\\Gamma^{\\mu}_{\\nu\\delta}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy4christoffels}$$\n",
    "\n",
    "By definition,\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": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:37.701275Z",
     "iopub.status.busy": "2021-03-07T17:07:37.685742Z",
     "iopub.status.idle": "2021-03-07T17:07:37.725098Z",
     "shell.execute_reply": "2021-03-07T17:07:37.724355Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.a.iii: Compute all 40 4-Christoffels \\Gamma^{\\mu}_{\\nu\\delta}\n",
    "\n",
    "betaDdD = ixp.zerorank2()\n",
    "gammaDD_dD = ixp.declarerank3(\"gammaDD_dD\",\"sym01\")\n",
    "betaU_dD   = ixp.declarerank2(\"betaU_dD\",\"nosym\")\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",
    "alpha_dD   = ixp.declarerank1(\"alpha_dD\")\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": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:37.759881Z",
     "iopub.status.busy": "2021-03-07T17:07:37.743394Z",
     "iopub.status.idle": "2021-03-07T17:07:50.639945Z",
     "shell.execute_reply": "2021-03-07T17:07:50.639318Z"
    }
   },
   "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),order=\"4\"))\n",
    "            interp_expr = Gamma4UDD[mu][nu][delta]\n",
    "            which_InterpCounter = interp_fileout(which_InterpCounter,interp_expr,NRPyoutfilename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy_c_callingfunction'></a>\n",
    "\n",
    "## Step 2.b: C code calling function for the NRPy+ C output \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy_c_callingfunction}$$\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": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:50.645834Z",
     "iopub.status.busy": "2021-03-07T17:07:50.644970Z",
     "iopub.status.idle": "2021-03-07T17:07:50.648759Z",
     "shell.execute_reply": "2021-03-07T17:07:50.648218Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_sph_grids_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 *cctkGH,const CCTK_INT *cctk_lsh,const CCTK_INT *cctk_nghostzones,\n",
    "                                     const CCTK_REAL invdx0,const CCTK_REAL invdx1,const CCTK_REAL invdx2,\n",
    "                                     const CCTK_INT *InterpCounter,\n",
    "                                     const CCTK_REAL *rho_bGF,const CCTK_REAL *PGF,\n",
    "                                     const CCTK_REAL *IGMvU0GF,const CCTK_REAL *IGMvU1GF,const CCTK_REAL *IGMvU2GF,\n",
    "                                     const CCTK_REAL *BU0GF,const CCTK_REAL *BU1GF,const CCTK_REAL *BU2GF,\n",
    "                                     const CCTK_REAL *gammaDD00GF,const CCTK_REAL *gammaDD01GF,const CCTK_REAL *gammaDD02GF,\n",
    "                                     const CCTK_REAL *gammaDD11GF,const CCTK_REAL *gammaDD12GF,const CCTK_REAL *gammaDD22GF,\n",
    "                                     const CCTK_REAL *betaU0GF,const CCTK_REAL *betaU1GF,const CCTK_REAL *betaU2GF,\n",
    "                                     const CCTK_REAL *alphaGF,   CCTK_REAL *interped_gfGF) {\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",
    "  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",
    "  list_of_functions_to_interpolate(cctkGH,cctk_lsh,cctk_nghostzones,invdx0,invdx1,invdx2,\n",
    "                                   InterpCounter,\n",
    "                                   rho_b,P,\n",
    "                                   vx,vy,vz,\n",
    "                                   Bx,By,Bz,\n",
    "                                   gxx,gxy,gxz,gyy,gyz,gzz,\n",
    "                                   betax,betay,betaz,alp, interped_gf);\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='nrpygetgfinterporder'></a>\n",
    "\n",
    "## Step 2.c: `get_gf_name()` and `get_interp_order()` functions \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpygetgfinterporder}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:50.657251Z",
     "iopub.status.busy": "2021-03-07T17:07:50.656396Z",
     "iopub.status.idle": "2021-03-07T17:07:50.659393Z",
     "shell.execute_reply": "2021-03-07T17:07:50.659926Z"
    }
   },
   "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].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\")\n",
    "\n",
    "with open(os.path.join(Ccodesdir,\"src\",\"get_interp_order.h\"), \"w\") as file:\n",
    "    file.write(\"int get_interp_order(const int InterpCounter) {\\n\")\n",
    "    for i in range(1,which_InterpCounter):\n",
    "        file.write(\"    if(InterpCounter==\"+str(i)+\") return \"+gf_interp_list[i].order+\";\\n\")\n",
    "    file.write(\"    return -1;\\n}\\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": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:50.665134Z",
     "iopub.status.busy": "2021-03-07T17:07:50.664435Z",
     "iopub.status.idle": "2021-03-07T17:07:50.666906Z",
     "shell.execute_reply": "2021-03-07T17:07:50.667495Z"
    }
   },
   "outputs": [],
   "source": [
    "with open(os.path.join(Ccodesdir,\"src\",\"define_NumInterpFunctions.h\"), \"w\") as file:\n",
    "    file.write(\"#define NumInterpFunctions \"+str(which_InterpCounter)+\"\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:50.673874Z",
     "iopub.status.busy": "2021-03-07T17:07:50.672595Z",
     "iopub.status.idle": "2021-03-07T17:07:50.677655Z",
     "shell.execute_reply": "2021-03-07T17:07:50.676910Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_sph_grids_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 SphGrid_InitializeInterpCounterToZero(CCTK_ARGUMENTS)\n",
    "{\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  *InterpCounter = 0;\n",
    "\n",
    "  if(verbose==2) printf(\"interp_sph_grids_ETK: Just set InterpCounter to %d\\n\",*InterpCounter);\n",
    "}\n",
    "\n",
    "void SphGrid_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_sph_grids_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 SphGrid_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_sph_grids_ETK: Finished! Just zeroed InterpCounter.\\n\");\n",
    "    } else {\n",
    "        (*InterpCounter)++;\n",
    "        if(verbose==2) printf(\"interp_sph_grids_ETK: Just incremented InterpCounter to %d of %d\\n\",*InterpCounter,NumInterpFunctions-1);\n",
    "    }\n",
    "}"
   ]
  },
  {
   "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": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:50.683843Z",
     "iopub.status.busy": "2021-03-07T17:07:50.682889Z",
     "iopub.status.idle": "2021-03-07T17:07:50.686974Z",
     "shell.execute_reply": "2021-03-07T17:07:50.687726Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_sph_grids_ETK/src/make.code.defn\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/src/make.code.defn\n",
    "# Main make.code.defn file for thorn interp_sph_grids_ETK\n",
    "\n",
    "# Source files in this directory\n",
    "SRCS =  main_function.cc interp_counter.cc construct_function_to_interpolate__store_to_interped_gf.cc"
   ]
  },
  {
   "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": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:50.694976Z",
     "iopub.status.busy": "2021-03-07T17:07:50.693890Z",
     "iopub.status.idle": "2021-03-07T17:07:50.697718Z",
     "shell.execute_reply": "2021-03-07T17:07:50.698394Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_sph_grids_ETK/interface.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/interface.ccl\n",
    "\n",
    "# With \"implements\", we give our thorn its unique name.\n",
    "implements: interp_sph_grids_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",
    "\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",
    "int InterpCounterVar type = SCALAR tags='checkpoint=\"no\"'\n",
    "{\n",
    "  InterpCounter\n",
    "} \"Counter that keeps track of which function we are interpolating.\"\n",
    "\n",
    "CCTK_REAL interp_pointcoords_and_output_arrays TYPE=ARRAY DISTRIB=CONSTANT DIM=1 SIZE=N0*N1*N2 tags='checkpoint=\"no\"'\n",
    "{\n",
    "  points_x,points_y,points_z,\n",
    "  output_interped\n",
    "}"
   ]
  },
  {
   "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": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:50.705784Z",
     "iopub.status.busy": "2021-03-07T17:07:50.704632Z",
     "iopub.status.idle": "2021-03-07T17:07:50.708578Z",
     "shell.execute_reply": "2021-03-07T17:07:50.709407Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_sph_grids_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",
    "restricted:\n",
    "\n",
    "########################################\n",
    "# BASIC THORN STEERING PARAMETERS\n",
    "CCTK_INT interp_out_iteration \"Which iteration to interpolate to spherical grids?\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 960000\n",
    "\n",
    "## Interpolator information\n",
    "CCTK_STRING interpolator_name \"Which interpolator to use?\" STEERABLE=ALWAYS\n",
    "{\n",
    "  \".+\" :: \"Any nonempty string; an unsupported value will throw an error.\"\n",
    "} \"Lagrange polynomial interpolation\"\n",
    "\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",
    "# SPHERICAL COORDINATE SYSTEM PARAMETERS\n",
    "CCTK_INT N0 \"Number of points in r direction\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 96\n",
    "\n",
    "CCTK_INT N1 \"Number of points in theta direction\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 96\n",
    "\n",
    "CCTK_INT N2 \"Number of points in phi direction\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 96\n",
    "\n",
    "##########\n",
    "# Cartesian position of center of spherical grid (usually center of BH) -- CURRENTLY UNSUPPORTED!\n",
    "CCTK_REAL x_center \"x-position of center.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "CCTK_REAL y_center \"y-position of center.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "CCTK_REAL z_center \"z-position of center.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "##########\n",
    "# Radial parameters:\n",
    "CCTK_REAL R0 \"Radial offset: r(x0) = R_0 + exp(x0). Probably should keep it set to zero.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "CCTK_REAL Rin \"x0 offset: x0 = log(Rin-R0) + (i + 0.5)Dx0.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 1.08986052555408\n",
    "\n",
    "CCTK_REAL Rout \"Dx0 = log( (Rout-R0) / (Rin-R0) )/N0\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 80.0\n",
    "\n",
    "##########\n",
    "# Theta parameters:\n",
    "CCTK_REAL x1_beg \"x1 offset: x1 = x1_beg + (j + 0.5)Dx1. Probably should keep it set to zero.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 0.0\n",
    "\n",
    "CCTK_INT theta_option \"Which prescription for theta should be used? 1 or 2?\" STEERABLE=ALWAYS\n",
    "{\n",
    "  1:2 :: \"\"\n",
    "} 1\n",
    "\n",
    "CCTK_REAL th_c \"theta_c: Angular cutout size for theta = 0 and pi\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 0.053407075111026485 # 0.017*pi\n",
    "\n",
    "CCTK_REAL xi \"Amplitude of nonlinear part of the theta distribution.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 0.25\n",
    "\n",
    "CCTK_INT th_n \"Power of nonlinear part of theta distribution. Only for theta_option=2\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 9\n",
    "\n",
    "##########\n",
    "# Phi parameters:\n",
    "CCTK_REAL x2_beg \"x2 offset: x2 = x2_beg + (k + 0.5)Dx2. Probably should keep it set to zero.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  0:* :: \"\"\n",
    "} 0.0\n",
    "########################################"
   ]
  },
  {
   "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": 20,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:50.715829Z",
     "iopub.status.busy": "2021-03-07T17:07:50.714889Z",
     "iopub.status.idle": "2021-03-07T17:07:50.718375Z",
     "shell.execute_reply": "2021-03-07T17:07:50.718902Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing interp_sph_grids_ETK/schedule.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/schedule.ccl\n",
    "\n",
    "STORAGE: interpolation_gf[3]\n",
    "STORAGE: InterpCounterVar\n",
    "STORAGE: interp_pointcoords_and_output_arrays\n",
    "\n",
    "#############################\n",
    "SCHEDULE SphGrid_InitializeInterpCounterToZero AT CCTK_INITIAL\n",
    "{\n",
    "  LANG: C\n",
    "  OPTIONS: GLOBAL\n",
    "} \"Initialize InterpCounter variable to zero\"\n",
    "\n",
    "SCHEDULE SphGrid_InitializeInterpCounterToZero AT CCTK_POST_RECOVER_VARIABLES\n",
    "{\n",
    "  LANG: C\n",
    "  OPTIONS: GLOBAL\n",
    "} \"Initialize InterpCounter variable to zero\"\n",
    "\n",
    "SCHEDULE SphGrid_InitializeInterpCounter before SphGrid_InterpGroup AT CCTK_ANALYSIS\n",
    "{\n",
    "  LANG: C\n",
    "  OPTIONS: GLOBAL\n",
    "} \"Initialize InterpCounter variable\"\n",
    "##################\n",
    "\n",
    "SCHEDULE GROUP SphGrid_InterpGroup AT CCTK_ANALYSIS BEFORE CarpetLib_printtimestats BEFORE CarpetLib_printmemstats AFTER Convert_to_HydroBase WHILE interp_sph_grids_ETK::InterpCounter\n",
    "{\n",
    "} \"Perform all spherical interpolations. This group is only actually scheduled at cctk_iteration==interp_out_iteration.\"\n",
    "\n",
    "SCHEDULE construct_function_to_interpolate__store_to_interped_gf in SphGrid_InterpGroup before DoSum\n",
    "{\n",
    "  STORAGE: interpolation_gf[3],InterpCounterVar,interp_pointcoords_and_output_arrays\n",
    "  OPTIONS: GLOBAL,LOOP-LOCAL\n",
    "  SYNC: interpolation_gf\n",
    "  LANG: C\n",
    "} \"Construct the function to interpolate\"\n",
    "\n",
    "SCHEDULE Interpolate_to_sph_grid_main_function in SphGrid_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",
    "SCHEDULE SphGrid_IncrementInterpCounter in SphGrid_InterpGroup after Interpolate_to_sph_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_sph_grids.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_Sph_ReadIn.py\", and run it using the command\n",
    "\n",
    "**`python Interp_Sph_ReadIn.py interp_sph_grids.dat 58 outfile`**\n",
    "\n",
    "Currently the last parameter \"outfile\" is required but not used.\n",
    "\n",
    "```python\n",
    "\"\"\"\n",
    "interp_sph_grids.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_Spherical_Grids.ipynb\n",
    "\n",
    "Usage instructions:\n",
    "\n",
    "From the command-line, run via:\n",
    "python Interp_Sph_ReadIn.py interp_sph_grids.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 longer 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 radial grid parameters:\n",
    "    # fwrite( & N0, sizeof(CCTK_INT), 1, file);\n",
    "    N0 =    struct.unpack('i',filehandle.read(4))[0]\n",
    "    # fwrite( & R0, sizeof(CCTK_REAL), 1, file);\n",
    "    R0 =    struct.unpack('d',filehandle.read(8))[0]\n",
    "    # fwrite( & Rin, sizeof(CCTK_REAL), 1, file);\n",
    "    Rin =   struct.unpack('d',filehandle.read(8))[0]\n",
    "    # fwrite( & Rout, sizeof(CCTK_REAL), 1, file);\n",
    "    Rout =  struct.unpack('d',filehandle.read(8))[0]\n",
    "\n",
    "    # Then the grid parameters related to the theta coordinate:\n",
    "    # fwrite( & N1, sizeof(CCTK_INT), 1, file);\n",
    "    N1           = struct.unpack('i', filehandle.read(4))[0]\n",
    "    # fwrite( & x1_beg, sizeof(CCTK_REAL), 1, file);\n",
    "    x1_beg       = struct.unpack('d', filehandle.read(8))[0]\n",
    "    # fwrite( & theta_option, sizeof(CCTK_INT), 1, file);\n",
    "    theta_option = struct.unpack('i', filehandle.read(4))[0]\n",
    "    # fwrite( & th_c, sizeof(CCTK_REAL), 1, file);\n",
    "    th_c         = struct.unpack('d', filehandle.read(8))[0]\n",
    "    # fwrite( & xi, sizeof(CCTK_REAL), 1, file);\n",
    "    xi           = struct.unpack('d', filehandle.read(8))[0]\n",
    "    # fwrite( & th_n, sizeof(CCTK_INT), 1, file);\n",
    "    th_n         = struct.unpack('i', filehandle.read(4))[0]\n",
    "\n",
    "    # Then the grid parameters related to the phi coordinate:\n",
    "    # fwrite( & N2, sizeof(CCTK_INT), 1, file);\n",
    "    N2     = struct.unpack('i', filehandle.read(4))[0]\n",
    "    # fwrite( & x2_beg, sizeof(CCTK_REAL), 1, file);\n",
    "    x2_beg = struct.unpack('d', filehandle.read(8))[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,N0,R0,Rin,Rout,N1,x1_beg,theta_option,th_c,xi,th_n,N2,x2_beg,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, N0, R0, Rin, Rout, N1, x1_beg, theta_option, th_c, xi, th_n, N2, x2_beg, cctk_iteration, cctk_time = read_header(f)\n",
    "        print(\"\\n\\nReading gridfunction \"+gf_name)\n",
    "        data_chunk_size = N0*N1*N2*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: Make sure the output in the \"middle\" of the grid looks reasonable.\n",
    "        ii = int(N0/2)\n",
    "        jj = int(N1/2)\n",
    "        kk = int(N2/2)\n",
    "        with open(\"output-gf\"+str(i)+\".txt\",\"w\") as file:\n",
    "            for ii in range(N0):\n",
    "                for kk in range(N2):\n",
    "                    r  = ii*1.0/N0\n",
    "                    th = (jj*1.0)*np.pi/N1\n",
    "                    ph = (kk*1.0)*2.0*np.pi/N2\n",
    "                    xx = r*np.sin(th)*np.cos(ph)\n",
    "                    yy = r*np.sin(th)*np.sin(ph)\n",
    "                    zz = r*np.cos(th)\n",
    "                    file.write(str(xx)+\" \"+str(yy)+\" \"+str(zz)+\" \"+str(this_data[kk,jj,ii])+\"\\n\")\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_Spherical_Grids.pdf](Tutorial-ETK_thorn-Interpolation_to_Spherical_Grids.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": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:07:50.725876Z",
     "iopub.status.busy": "2021-03-07T17:07:50.723007Z",
     "iopub.status.idle": "2021-03-07T17:07:55.810834Z",
     "shell.execute_reply": "2021-03-07T17:07:55.809875Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ETK_thorn-Interpolation_to_Spherical_Grids.tex, and\n",
      "    compiled LaTeX file to PDF file Tutorial-ETK_thorn-\n",
      "    Interpolation_to_Spherical_Grids.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_Spherical_Grids\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}