{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Start-to-Finish Example: Validating `GiRaFFEfood_NRPy` against original, trusted `GiRaFFEfood`: \n",
    "\n",
    "## Author: Patrick Nelson\n",
    "\n",
    "## The notebook validates GiRaFFE_NRPy's initial data routines, known as GiRaFFEfood_NRPy, by contrasting them with the original GiRaFFEfood code, and assessing the level of numerical agreement.\n",
    "\n",
    "**Notebook Status:** <font color='green'><b>Validated</b></font>\n",
    "\n",
    "**Validation Notes:** This module validates all expressions used to set up initial data in \n",
    "* [Tutorial-GiRaFFEfood_NRPy-Exact_Wald](Tutorial-GiRaFFEfood_NRPy-Exact_Wald.ipynb)\n",
    "* [Tutorial-GiRaFFEfood_NRPy-Magnetospheric_Wald](Tutorial-GiRaFFEfood_NRPy-Magnetospheric_Wald.ipynb)\n",
    "* [Tutorial-GiRaFFEfood_NRPy-Split_Monopole](Tutorial-GiRaFFEfood_NRPy-Split_Monopole.ipynb)\n",
    "\n",
    "against the C-code implementation of these expressions found in the original (trusted) [`GiRaFFEfood` Einstein Toolkit thorn](link), and confirms roundoff-level agreement.\n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Exact_Wald.py](../../edit/in_progress/GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Exact_Wald.py) [\\[**tutorial**\\]](Tutorial-GiRaFFEfood_NRPy-Exact_Wald.ipynb) Generates Exact Wald initial data\n",
    "* [GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Magnetospheric_Wald.py](../../edit/in_progress/GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Magnetospheric_Wald.py) [\\[**tutorial**\\]](Tutorial-GiRaFFEfood_NRPy-Magnetospheric_Wald.ipynb) Generates Exact Wald initial data\n",
    "* [GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Split_Monopole.py](../../edit/in_progress/GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Split_Monopole.py) [\\[**tutorial**\\]](Tutorial-GiRaFFEfood_NRPy-Split_Monopole.ipynb) Generates Exact Wald initial data\n",
    "\n",
    "## Introduction:\n",
    "\n",
    "This notebook validates the initial data routines that we will use for `GiRaFFE_NRPy`, collectively referred to as `GiRaFFEfood_NRPy`. To do so, we will generate the initial data with both our code and the original `GiRaFFEfood` code. Then, we will directly compare the velocities and show round-off level agreement between the two. \n",
    "\n",
    "When this notebook is run, the significant digits of agreement between the old `GiRaFFE` and new `GiRaFFE_NRPy` versions of the algorithm will be evaluated. If the agreement falls below a thresold, the point, quantity, and level of agreement are reported [here](#compile_run).\n"
   ]
  },
  {
   "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](#setup): Set up core functions and parameters for unit testing the initial data algorithms\n",
    "    1. [Step 1.a](#spacetime) Generate the spacetime metric \n",
    "    1. [Step 1.b](#initial_data) Generate the initial data C function\n",
    "    1. [Step 1.c](#download) Download original `GiRaFFE` files\n",
    "    1. [Step 1.d](#free_params) Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h`\n",
    "    1. [Step 1.e](#interface) Create dummy files for the CCTK version of the code\n",
    "1. [Step 2](#mainc): `GiRaFFEfood_NRPy_unit_test.c`: The Main C Code\n",
    "    1. [Step 2.a](#compile_run): Compile and run the code to validate the output\n",
    "1. [Step 3](#drift_notes): [comment]: <> (TODO: drift notes)\n",
    "1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='setup'></a>\n",
    "\n",
    "# Step 1: Set up core functions and parameters for unit testing the initial data algorithms\" \\[Back to [top](#toc)\\]\n",
    "\n",
    "$$\\label{setup}$$\n",
    "\n",
    "We'll start by appending the relevant paths to `sys.path` so that we can access sympy modules in other places. Then, we'll import NRPy+ core functionality and set up a directory in which to carry out our test. We will also declare the gridfunctions that are needed for this portion of the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:08.747909Z",
     "iopub.status.busy": "2022-12-29T09:54:08.747613Z",
     "iopub.status.idle": "2022-12-29T09:54:08.756765Z",
     "shell.execute_reply": "2022-12-29T09:54:08.756249Z"
    }
   },
   "outputs": [],
   "source": [
    "# There are several initial data routines we need to test. We'll control which one we use with a string option\n",
    "initial_data = \"MagnetosphericWald\" # Valid options: \"AllTests\", \"ExactWald\", \"SplitMonopole\", \"MagnetosphericWald\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:08.758883Z",
     "iopub.status.busy": "2022-12-29T09:54:08.758658Z",
     "iopub.status.idle": "2022-12-29T09:54:09.301800Z",
     "shell.execute_reply": "2022-12-29T09:54:09.301262Z"
    }
   },
   "outputs": [],
   "source": [
    "import os, sys, shutil           # Standard Python modules for multiplatform OS-level functions\n",
    "# First, we'll add the parent directory to the list of directories Python will check for modules.\n",
    "nrpy_dir_path = os.path.join(\"..\")\n",
    "if nrpy_dir_path not in sys.path:\n",
    "    sys.path.append(nrpy_dir_path)\n",
    "nrpy_dir_path = os.path.join(\"..\", \"Deprecated\")\n",
    "if nrpy_dir_path not in sys.path:\n",
    "    sys.path.append(nrpy_dir_path)\n",
    "# nrpy_dir_path = os.path.join(\"..\",\"..\")\n",
    "# if nrpy_dir_path not in sys.path:\n",
    "#     sys.path.append(nrpy_dir_path)\n",
    "\n",
    "from outputC import outCfunction, lhrh # NRPy+: Core C code output module\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import finite_difference as fin  # NRPy+: Finite difference C code generation module\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import grid as gri               # NRPy+: Functions having to do with numerical grids\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface\n",
    "import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy_Common_Functions as gfcf # Some useful functions for GiRaFFE initial data.\n",
    "\n",
    "out_dir = os.path.join(\"Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Curved_Spacetime_Tests/\")\n",
    "# First remove C code output directory if it exists\n",
    "# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
    "shutil.rmtree(out_dir, ignore_errors=True)\n",
    "cmd.mkdir(out_dir)\n",
    "\n",
    "thismodule = \"Start_to_Finish_UnitTest-GiRaFFEfood_NRPy\"\n",
    "\n",
    "# Register the gridfunctions we need for this function\n",
    "AD = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"AD\")\n",
    "ValenciavU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"ValenciavU\")\n",
    "# gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gammaDD\",\"sym01\")\n",
    "betaU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"betaU\")\n",
    "betaSphU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"betaSphU\")\n",
    "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gammaDD\",\"sym01\")\n",
    "gammaSphDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gammaSphDD\",\"sym01\")\n",
    "alpha = gri.register_gridfunctions(\"AUXEVOL\",\"alpha\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='spacetime'></a>\n",
    "\n",
    "## Step 1.a: Generate the spacetime metric \\[Back to [top](#toc)\\]\n",
    "$$\\label{spacetime}$$\n",
    "\n",
    "While many of the initial data we will use assume a flat background spacetime, some will require a specific metric. We will set those up as needed here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:09.303960Z",
     "iopub.status.busy": "2022-12-29T09:54:09.303751Z",
     "iopub.status.idle": "2022-12-29T09:54:09.573032Z",
     "shell.execute_reply": "2022-12-29T09:54:09.572364Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function Shifted_Kerr_Schild_initial_metric_spherical() to file Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Curved_Spacetime_Tests/Shifted_Kerr_Schild_initial_metric_spherical.h\n",
      "Output C function Shifted_Kerr_Schild_initial_metric() to file Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Curved_Spacetime_Tests/Shifted_Kerr_Schild_initial_metric.h\n"
     ]
    }
   ],
   "source": [
    "# Exact Wald is more complicated. We'll need the Shifted Kerr Schild metric in Cartesian coordinates.\n",
    "import BSSN.ShiftedKerrSchild as sks\n",
    "sks.ShiftedKerrSchild()\n",
    "import reference_metric as rfm\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n",
    "rfm.reference_metric()\n",
    "# Use the Jacobian matrix to transform the vectors to Cartesian coordinates.\n",
    "gammaSphDD = ixp.zerorank2()\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        gammaSphDD[i][j] += sks.gammaDD[i][j].subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])\n",
    "\n",
    "unused_gammaUU,gammaSphDET = ixp.symm_matrix_inverter3x3(gammaSphDD)\n",
    "sqrtgammaSphDET = sp.sqrt(gammaSphDET)\n",
    "\n",
    "betaSphU = ixp.zerorank1()\n",
    "for i in range(3):\n",
    "    betaSphU[i] += sks.betaU[i].subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])\n",
    "alpha = sks.alpha.subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])\n",
    "\n",
    "# We only need to set alpha and betaU in C for the original Exact Wald\n",
    "name = \"Shifted_Kerr_Schild_initial_metric_spherical\"\n",
    "desc = \"Generate a spinning black hole with Shifted Kerr Schild metric in a spherical basis.\"\n",
    "values_to_print = [\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaSphDD00\"),rhs=gammaSphDD[0][0]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaSphDD01\"),rhs=gammaSphDD[0][1]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaSphDD02\"),rhs=gammaSphDD[0][2]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaSphDD11\"),rhs=gammaSphDD[1][1]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaSphDD12\"),rhs=gammaSphDD[1][2]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaSphDD22\"),rhs=gammaSphDD[2][2]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"betaSphU0\"),rhs=betaSphU[0]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"betaSphU1\"),rhs=betaSphU[1]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"betaSphU2\"),rhs=betaSphU[2]),\n",
    "                  ]\n",
    "\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(out_dir,name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs\",\n",
    "    body     = fin.FD_outputC(\"returnstring\",values_to_print,params=\"outCverbose=False\"),\n",
    "    loopopts =\"AllPoints,Read_xxs\")\n",
    "\n",
    "gammaDD = rfm.basis_transform_tensorDD_from_rfmbasis_to_Cartesian(gfcf.Jac_dUrfm_dDCartUD, gammaSphDD)\n",
    "\n",
    "unused_gammaUU,gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)\n",
    "sqrtgammaDET = sp.sqrt(gammaDET)\n",
    "\n",
    "betaU = rfm.basis_transform_vectorD_from_rfmbasis_to_Cartesian(gfcf.Jac_dUrfm_dDCartUD, betaSphU)\n",
    "\n",
    "# We only need to set alpha and betaU in C for the original Exact Wald\n",
    "name = \"Shifted_Kerr_Schild_initial_metric\"\n",
    "desc = \"Generate a spinning black hole with Shifted Kerr Schild metric in a Cartesian basis.\"\n",
    "values_to_print = [\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaDD00\"),rhs=gammaDD[0][0]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaDD01\"),rhs=gammaDD[0][1]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaDD02\"),rhs=gammaDD[0][2]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaDD11\"),rhs=gammaDD[1][1]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaDD12\"),rhs=gammaDD[1][2]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"gammaDD22\"),rhs=gammaDD[2][2]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"betaU0\"),rhs=betaU[0]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"betaU1\"),rhs=betaU[1]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"betaU2\"),rhs=betaU[2]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"alpha\"),rhs=sks.alpha.subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1]))\n",
    "                  ]\n",
    "\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(out_dir,name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs\",\n",
    "    body     = fin.FD_outputC(\"returnstring\",values_to_print,params=\"outCverbose=False\"),\n",
    "    loopopts =\"AllPoints,Read_xxs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initial_data'></a>\n",
    "\n",
    "## Step 1.b: Generate the initial data C function \\[Back to [top](#toc)\\]\n",
    "$$\\label{initial_data}$$\n",
    "\n",
    "First, we'll use NRPy+ to build the C function that will generate the initial data. There are several different cases here, one for each type of initial test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:09.602259Z",
     "iopub.status.busy": "2022-12-29T09:54:09.601747Z",
     "iopub.status.idle": "2022-12-29T09:54:09.928747Z",
     "shell.execute_reply": "2022-12-29T09:54:09.928230Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function GiRaFFE_NRPy_initial_data() to file Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Curved_Spacetime_Tests/GiRaFFE_NRPy_initial_data.h\n"
     ]
    }
   ],
   "source": [
    "import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy as gid\n",
    "if initial_data==\"ExactWald\":\n",
    "    gid.GiRaFFEfood_NRPy_generate_initial_data(ID_type = initial_data, stagger_enable = True,M=sks.M,KerrSchild_radial_shift=sks.r0,gammaDD=gammaDD,sqrtgammaDET=sqrtgammaDET)\n",
    "    desc = \"Generate exact Wald initial test data for GiRaFFEfood_NRPy.\"\n",
    "elif initial_data==\"MagnetosphericWald\":\n",
    "    import BSSN.ADMBSSN_tofrom_4metric as AB4m\n",
    "    AB4m.g4DD_ito_BSSN_or_ADM(\"ADM\",sks.gammaDD,sks.betaU,sks.alpha)\n",
    "    g4DD = AB4m.g4DD\n",
    "\n",
    "    for mu in range(4):\n",
    "        for nu in range(4):\n",
    "            g4DD[mu][nu] = g4DD[mu][nu].subs(sks.r,rfm.xxSph[0]).subs(sks.th,rfm.xxSph[1])\n",
    "\n",
    "    gid.GiRaFFEfood_NRPy_generate_initial_data(ID_type = initial_data, stagger_enable = True,a=sks.a,g4DD = g4DD)\n",
    "    desc = \"Generate exact Wald initial test data for GiRaFFEfood_NRPy.\"\n",
    "elif initial_data==\"SplitMonopole\":\n",
    "    gid.GiRaFFEfood_NRPy_generate_initial_data(ID_type = initial_data, stagger_enable = True,M=sks.M,a=sks.a,KerrSchild_radial_shift=sks.r0,alpha=alpha,betaU=betaSphU,gammaDD=gammaDD,sqrtgammaDET=sqrtgammaSphDET)\n",
    "    desc = \"Generate Split Monopole initial test data for GiRaFFEfood_NRPy.\"\n",
    "else:\n",
    "    print(\"Unsupported Initial Data string \"+initial_data+\"! Supported ID: AllTests, ExactWald, SplitMonopole, or Magnetospheric Wald\")\n",
    "\n",
    "name = \"GiRaFFE_NRPy_initial_data\"\n",
    "\n",
    "values_to_print = [\n",
    "                   lhrh(lhs=gri.gfaccess(\"out_gfs\",\"AD0\"),rhs=gid.AD[0]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"out_gfs\",\"AD1\"),rhs=gid.AD[1]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"out_gfs\",\"AD2\"),rhs=gid.AD[2]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"ValenciavU0\"),rhs=gid.ValenciavU[0]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"ValenciavU1\"),rhs=gid.ValenciavU[1]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"ValenciavU2\"),rhs=gid.ValenciavU[2])\n",
    "                  ]\n",
    "\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(out_dir,name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs,REAL *out_gfs\",\n",
    "    body     = fin.FD_outputC(\"returnstring\",values_to_print,params=\"outCverbose=False\"),\n",
    "    loopopts =\"AllPoints,Read_xxs\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='download'></a>\n",
    "\n",
    "## Step 1.b: Download original `GiRaFFE` files \\[Back to [top](#toc)\\]\n",
    "\n",
    "$$\\label{download}$$\n",
    "\n",
    "Here, we download the relevant portion of the original `GiRaFFE` code from Bitbucket. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:09.931262Z",
     "iopub.status.busy": "2022-12-29T09:54:09.931129Z",
     "iopub.status.idle": "2022-12-29T09:54:10.738188Z",
     "shell.execute_reply": "2022-12-29T09:54:10.737553Z"
    }
   },
   "outputs": [],
   "source": [
    "# First download the original GiRaFFE source code\n",
    "import urllib\n",
    "\n",
    "original_file_url  = [\n",
    "                      \"https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFEfood/src/ExactWald.cc\",\n",
    "                      \"https://bitbucket.org/zach_etienne/wvuthorns/raw/8a2f5662114af154d1915ee8460e91769518e05d/GiRaFFEfood/src/MagnetoWald.cc\",\n",
    "                      \"https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFEfood/src/SplitMonopole.cc\"\n",
    "                     ]\n",
    "original_file_name = [\n",
    "                      \"ExactWald.cc\",\n",
    "                      \"MagnetoWald.cc\",\n",
    "                      \"SplitMonopole.cc\"\n",
    "                     ]\n",
    "\n",
    "for i in range(len(original_file_url)):\n",
    "    original_file_path = os.path.join(out_dir,original_file_name[i])\n",
    "\n",
    "    # Then download the original GiRaFFE source code\n",
    "    # We try it here in a couple of ways in an attempt to keep\n",
    "    # the code more portable\n",
    "    try:\n",
    "        original_file_code = urllib.request.urlopen(original_file_url[i]).read().decode('utf-8')\n",
    "    except:\n",
    "        original_file_code = urllib.urlopen(original_file_url[i]).read().decode('utf-8')\n",
    "\n",
    "    # Write down the file the original GiRaFFE source code\n",
    "    with open(original_file_path,\"w\") as file:\n",
    "        file.write(original_file_code)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='free_params'></a>\n",
    "\n",
    "## Step 1.c: Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h` \\[Back to [top](#toc)\\]\n",
    "\n",
    "$$\\label{free_params}$$\n",
    "\n",
    "Based on declared NRPy+ Cparameters, first we generate `declare_Cparameters_struct.h`, `set_Cparameters_default.h`, and `set_Cparameters[-SIMD].h`.\n",
    "\n",
    "Then we output `free_parameters.h`, which sets some basic grid parameters as well as the speed limit parameter we need for this function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:10.741483Z",
     "iopub.status.busy": "2022-12-29T09:54:10.740604Z",
     "iopub.status.idle": "2022-12-29T09:54:10.752862Z",
     "shell.execute_reply": "2022-12-29T09:54:10.751859Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.d\n",
    "# Step 3.d.ii: Set free_parameters.h\n",
    "with open(os.path.join(out_dir,\"free_parameters.h\"),\"w\") as file:\n",
    "    file.write(\"\"\"\n",
    "// Set free-parameter values.\n",
    "\n",
    "const int NGHOSTS = 3;\n",
    "\n",
    "// Set free-parameter values for the initial data.\n",
    "// Override parameter defaults with values based on command line arguments and NGHOSTS.\n",
    "const int Nx0x1x2 = 5;\n",
    "params.Nxx0 = Nx0x1x2;\n",
    "params.Nxx1 = Nx0x1x2;\n",
    "params.Nxx2 = Nx0x1x2;\n",
    "params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*NGHOSTS;\n",
    "params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*NGHOSTS;\n",
    "params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*NGHOSTS;\n",
    "// Step 0d: Set up space and time coordinates\n",
    "// Step 0d.i: Declare \\Delta x^i=dxx{0,1,2} and invdxx{0,1,2}, as well as xxmin[3] and xxmax[3]:\n",
    "const REAL xxmin[3] = {-1.0,-1.0,-1.0};\n",
    "const REAL xxmax[3] = { 1.0, 1.0, 1.0};\n",
    "\n",
    "params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx0);\n",
    "params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx1);\n",
    "params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx2);\n",
    "params.invdx0 = 1.0 / params.dxx0;\n",
    "params.invdx1 = 1.0 / params.dxx1;\n",
    "params.invdx2 = 1.0 / params.dxx2;\n",
    "\\n\"\"\")\n",
    "\n",
    "if initial_data==\"ExactWald\" or initial_data==\"AllTests\":\n",
    "    with open(os.path.join(out_dir,\"free_parameters.h\"),\"a\") as file:\n",
    "        file.write(\"\"\"params.r0 = 0.4;\n",
    "params.a = 0.0;\n",
    "\"\"\")\n",
    "if initial_data==\"SplitMonopole\" or initial_data==\"AllTests\":\n",
    "    with open(os.path.join(out_dir,\"free_parameters.h\"),\"a\") as file:\n",
    "        file.write(\"\"\"params.r0 = 0.4;\n",
    "params.a = 0.1;\n",
    "\"\"\")\n",
    "if initial_data==\"MagnetosphericWald\" or initial_data==\"AllTests\":\n",
    "    with open(os.path.join(out_dir,\"free_parameters.h\"),\"a\") as file:\n",
    "        file.write(\"\"\"params.r0 = 0.4359;\n",
    "params.a = 0.0;\n",
    "\"\"\")\n",
    "\n",
    "\n",
    "# Generates declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n",
    "import deprecated_NRPy_param_funcs as evil_par\n",
    "evil_par.generate_Cparameters_Ccodes(os.path.join(out_dir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='interface'></a>\n",
    "\n",
    "## Step 1.d: Create dummy files for the CCTK version of the code \\[Back to [top](#toc)\\]\n",
    "\n",
    "$$\\label{interface}$$\n",
    "\n",
    "The original `GiRaFFE` code depends on some functionalities of the CCTK. Since we only care about this one small function, we can get around this by creating some nearly-empty, non-functional files that can be included to satisfy the pre-processor without changing functionality. We will later replace what little functionality we need with some basic global variables and macros."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:10.756296Z",
     "iopub.status.busy": "2022-12-29T09:54:10.756155Z",
     "iopub.status.idle": "2022-12-29T09:54:10.760512Z",
     "shell.execute_reply": "2022-12-29T09:54:10.759572Z"
    }
   },
   "outputs": [],
   "source": [
    "#incldue \"cctk.h\"\n",
    "#include \"cctk_Arguments.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "#include \"Symmetry.h\"\n",
    "with open(os.path.join(out_dir,\"cctk.h\"),\"w\") as file:\n",
    "    file.write(\"\"\"//\"\"\")\n",
    "\n",
    "with open(os.path.join(out_dir,\"cctk_Arguments.h\"),\"w\") as file:\n",
    "    file.write(\"\"\"#define DECLARE_CCTK_ARGUMENTS //\n",
    "#define CCTK_ARGUMENTS void\n",
    "\"\"\")\n",
    "\n",
    "with open(os.path.join(out_dir,\"cctk_Parameters.h\"),\"w\") as file:\n",
    "    file.write(\"\"\"#define DECLARE_CCTK_PARAMETERS //\n",
    "\"\"\")\n",
    "\n",
    "with open(os.path.join(out_dir,\"Symmetry.h\"),\"w\") as file:\n",
    "    file.write(\"\"\"//\"\"\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mainc'></a>\n",
    "\n",
    "# Step 2: `GiRaFFEfood_NRPy_unit_test.c`: The Main C Code \\[Back to [top](#toc)\\]\n",
    "\n",
    "$$\\label{mainc}$$\n",
    "\n",
    "Now that we have our vector potential and analytic magnetic field to compare against, we will start writing our unit test. We'll also import common C functionality, define `REAL`, the number of ghost zones, and the faces, and set the standard macros for NRPy+ style memory access."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:10.765456Z",
     "iopub.status.busy": "2022-12-29T09:54:10.765029Z",
     "iopub.status.idle": "2022-12-29T09:54:10.774928Z",
     "shell.execute_reply": "2022-12-29T09:54:10.774415Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Curved_Spacetime_Tests//GiRaFFEfood_NRPy_unit_test.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile $out_dir/GiRaFFEfood_NRPy_unit_test.C\n",
    "\n",
    "// These are common packages that we are likely to need.\n",
    "#include \"stdio.h\"\n",
    "#include \"stdlib.h\"\n",
    "#include \"math.h\"\n",
    "#include <string> // Needed for strncmp, etc.\n",
    "#include \"stdint.h\" // Needed for Windows GCC 6.x compatibility\n",
    "#include <time.h>   // Needed to set a random seed.\n",
    "#include \"gsl/gsl_sf_dilog.h\"\n",
    "\n",
    "#define REAL double\n",
    "#include \"declare_Cparameters_struct.h\"\n",
    "\n",
    "// Standard NRPy+ memory access:\n",
    "#define IDX4S(g,i,j,k) \\\n",
    "( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )\n",
    "\n",
    "// Standard formula to calculate significant digits of agreement:\n",
    "//#define SDA(a,b) 1.0-log10(2.0*fabs(a-b)/(fabs(a)+fabs(b)))\n",
    "REAL SDA(REAL a, REAL b) {\n",
    "    return (a<1.0e-14 && b<1.0e-14) ? 16.0 : 1.0-log10(2.0*fabs(a-b)/(fabs(a)+fabs(b)));\n",
    "}\n",
    "\n",
    "// Memory access definitions for NRPy+\n",
    "#define BU0GF 0\n",
    "#define BU1GF 1\n",
    "#define BU2GF 2\n",
    "#define VALENCIAVU0GF 3\n",
    "#define VALENCIAVU1GF 4\n",
    "#define VALENCIAVU2GF 5\n",
    "#define BETAU0GF 6\n",
    "#define BETAU1GF 7\n",
    "#define BETAU2GF 8\n",
    "#define ALPHAGF 9\n",
    "#define GAMMADD00GF 10\n",
    "#define GAMMADD01GF 11\n",
    "#define GAMMADD02GF 12\n",
    "#define GAMMADD11GF 13\n",
    "#define GAMMADD12GF 14\n",
    "#define GAMMADD22GF 15\n",
    "#define BETASPHU0GF 16\n",
    "#define BETASPHU1GF 17\n",
    "#define BETASPHU2GF 18\n",
    "#define GAMMASPHDD00GF 19\n",
    "#define GAMMASPHDD01GF 20\n",
    "#define GAMMASPHDD02GF 21\n",
    "#define GAMMASPHDD11GF 22\n",
    "#define GAMMASPHDD12GF 23\n",
    "#define GAMMASPHDD22GF 24\n",
    "#define NUM_AUXEVOL_GFS 25\n",
    "\n",
    "#define AD0GF 0\n",
    "#define AD1GF 1\n",
    "#define AD2GF 2\n",
    "#define PSI6PHIGF 3\n",
    "#define NUM_EVOL_GFS 4\n",
    "\n",
    "// Include the functions that we want to test:\n",
    "#include \"GiRaFFE_NRPy_initial_data.h\"\n",
    "#include \"Shifted_Kerr_Schild_initial_metric.h\"\n",
    "#include \"Shifted_Kerr_Schild_initial_metric_spherical.h\"\n",
    "\n",
    "// Define CCTK macros\n",
    "#define CCTK_REAL double\n",
    "#define CCTK_INT int\n",
    "struct cGH{};\n",
    "const cGH* cctkGH;\n",
    "\n",
    "// GiRaFFE parameters in ETK\n",
    "const CCTK_REAL min_radius_inside_of_which_conserv_to_prims_FFE_and_FFE_evolution_is_DISABLED = -1;\n",
    "const int current_sheet_null_v = 1;\n",
    "\n",
    "// More definitions to interface with ETK code:\n",
    "const int cctk_lsh[3] = {11,11,11};\n",
    "const int grid_size = cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2];\n",
    "const int grid_size_3_comp = grid_size*3;\n",
    "CCTK_REAL *Avec = (REAL *)malloc(sizeof(REAL) * grid_size_3_comp);\n",
    "CCTK_REAL *vel = (REAL *)malloc(sizeof(REAL) * grid_size_3_comp);\n",
    "CCTK_REAL *Aphi = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *psi6phi = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *Ax = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *Ay = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *Az = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *vx = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *vy = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *vz = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *Bx = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *By = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *Bz = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *x = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *y = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *z = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *r = (REAL *)malloc(sizeof(REAL) * grid_size);;\n",
    "CCTK_REAL *alp;\n",
    "CCTK_REAL *betax;\n",
    "CCTK_REAL *betay;\n",
    "CCTK_REAL *betaz;\n",
    "CCTK_REAL *gxx;\n",
    "CCTK_REAL *gxy;\n",
    "CCTK_REAL *gxz;\n",
    "CCTK_REAL *gyy;\n",
    "CCTK_REAL *gyz;\n",
    "CCTK_REAL *gzz;\n",
    "CCTK_REAL *SKSbetar;\n",
    "CCTK_REAL *SKSbetath;\n",
    "CCTK_REAL *SKSbetaph;\n",
    "CCTK_REAL *SKSgrr;\n",
    "CCTK_REAL *SKSgrth;\n",
    "CCTK_REAL *SKSgrph;\n",
    "CCTK_REAL *SKSgthth;\n",
    "CCTK_REAL *SKSgthph;\n",
    "CCTK_REAL *SKSgphph;\n",
    "\n",
    "// We need to declare these to compile functions we won't call:\n",
    "int Compute_Exact_Every;\n",
    "int cctk_iteration;\n",
    "CCTK_REAL *delpsi6phi;\n",
    "CCTK_REAL *delAx;\n",
    "CCTK_REAL *delAy;\n",
    "CCTK_REAL *delAz;\n",
    "CCTK_REAL *exactBx;\n",
    "CCTK_REAL *exactBy;\n",
    "CCTK_REAL *exactBz;\n",
    "CCTK_REAL *delBx;\n",
    "CCTK_REAL *delBy;\n",
    "CCTK_REAL *delBz;\n",
    "CCTK_REAL *exactVx;\n",
    "CCTK_REAL *exactVy;\n",
    "CCTK_REAL *exactVz;\n",
    "CCTK_REAL *delvx;\n",
    "CCTK_REAL *delvy;\n",
    "CCTK_REAL *delvz;\n",
    "CCTK_REAL cctk_time = 0.0;\n",
    "CCTK_REAL *exactBx_ThreeWaves;\n",
    "CCTK_REAL *exactBy_ThreeWaves;\n",
    "CCTK_REAL *exactBz_ThreeWaves;\n",
    "CCTK_REAL *delBx_ThreeWaves;\n",
    "CCTK_REAL *delBy_ThreeWaves;\n",
    "CCTK_REAL *delBz_ThreeWaves;\n",
    "CCTK_REAL *mhd_st_x;\n",
    "CCTK_REAL *mhd_st_y;\n",
    "CCTK_REAL *mhd_st_z;\n",
    "CCTK_REAL *Ex;\n",
    "CCTK_REAL *Ey;\n",
    "CCTK_REAL *Ez;\n",
    "CCTK_REAL *B2mE2;\n",
    "\n",
    "// Set constants to default for comparison\n",
    "CCTK_REAL wave_speed = -0.5;\n",
    "CCTK_REAL Omega_aligned_rotator = 1e3;\n",
    "CCTK_REAL R_NS_aligned_rotator = 1.0;\n",
    "CCTK_REAL B_p_aligned_rotator = 1e-5;\n",
    "CCTK_REAL Wald_B0 = 1.0;\n",
    "CCTK_REAL KerrSchild_radial_shift = 0.4;\n",
    "CCTK_REAL BH_mass = 1.0;\n",
    "CCTK_REAL BH_spin = 0.0;\n",
    "CCTK_REAL split_C = 1.0;\n",
    "const int drop_fr_SplitM = 1; // We don't want to set f(r)=0\n",
    "CCTK_REAL GAMMA_SPEED_LIMIT = 2000.0; // Set to default--only used for diagnostics.\n",
    "\n",
    "// Define dz in CCTK\n",
    "CCTK_REAL cactus_dxx[3];\n",
    "#define CCTK_DELTA_SPACE(i) cactus_dxx[i]\n",
    "\n",
    "// Dummy ETK function:\n",
    "#define CCTK_GFINDEX3D(cctkGH,i,j,k) (i) + cctk_lsh[0] * ( (j) + cctk_lsh[1] * (k) )\n",
    "#define CCTK_GFINDEX4D(cctkGH,i,j,k,g) \\\n",
    "( (i) + cctk_lsh[0] * ( (j) + cctk_lsh[1] * ( (k) + cctk_lsh[2] * (g) ) ) )\n",
    "#define CCTK_VInfo(...) //\n",
    "#define CCTK_VWarn(...) //\n",
    "#define CCTK_VError(...) //\n",
    "\n",
    "#include \"ExactWald.cc\"\n",
    "#include \"MagnetoWald.cc\"\n",
    "#include \"SplitMonopole.cc\"\n",
    "\n",
    "int main(int argc, char** argv) {\n",
    "    paramstruct params;\n",
    "#include \"set_Cparameters_default.h\"\n",
    "\n",
    "    // Step 0c: Set free parameters, overwriting Cparameters defaults\n",
    "    //          by hand or with command-line input, as desired.\n",
    "#include \"free_parameters.h\"\n",
    "#include \"set_Cparameters-nopointer.h\"\n",
    "\n",
    "    // Now that we've calculated dxx2,  we can define a cactus equivalent\n",
    "    cactus_dxx[0] = dxx0;\n",
    "    cactus_dxx[1] = dxx1;\n",
    "    cactus_dxx[2] = dxx2;\n",
    "\n",
    "    BH_spin = a;\n",
    "    Wald_B0 = C0;\n",
    "    KerrSchild_radial_shift = r0;\n",
    "\n",
    "    // Step 0d.ii: Set up uniform coordinate grids\n",
    "    REAL *xx[3];\n",
    "    xx[0] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS0);\n",
    "    xx[1] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS1);\n",
    "    xx[2] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS2);\n",
    "    for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] = xxmin[0] + (j-NGHOSTS)*dxx0;\n",
    "    for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + (j-NGHOSTS)*dxx1;\n",
    "    for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + (j-NGHOSTS)*dxx2;\n",
    "\n",
    "    for(int k=0;k<Nxx_plus_2NGHOSTS2;k++)\n",
    "        for(int j=0;j<Nxx_plus_2NGHOSTS1;j++)\n",
    "            for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) {\n",
    "                int index = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "                x[index] = xx[0][i];\n",
    "                y[index] = xx[1][j];\n",
    "                z[index] = xx[2][k];\n",
    "                r[index] = sqrt(x[index]*x[index] + y[index]*y[index] + z[index]*z[index]);\n",
    "    }\n",
    "\n",
    "    //for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) printf(\"x[%d] = %.5e\\n\",j,xx[0][j]);\n",
    "\n",
    "    // This is the array to which we'll write the NRPy+ variables.\n",
    "    REAL *auxevol_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);\n",
    "    REAL *evol_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);\n",
    "\n",
    "    // Memory access for metric gridfunctions for Exact Wald:\n",
    "    alp = auxevol_gfs + (grid_size*ALPHAGF);\n",
    "    betax = auxevol_gfs + (grid_size*BETAU0GF);\n",
    "    betay = auxevol_gfs + (grid_size*BETAU1GF);\n",
    "    betaz = auxevol_gfs + (grid_size*BETAU2GF);\n",
    "    gxx = auxevol_gfs + (grid_size*GAMMADD00GF);\n",
    "    gxy = auxevol_gfs + (grid_size*GAMMADD01GF);\n",
    "    gxz = auxevol_gfs + (grid_size*GAMMADD02GF);\n",
    "    gyy = auxevol_gfs + (grid_size*GAMMADD11GF);\n",
    "    gyz = auxevol_gfs + (grid_size*GAMMADD12GF);\n",
    "    gzz = auxevol_gfs + (grid_size*GAMMADD22GF);\n",
    "    SKSbetar = auxevol_gfs + (grid_size*BETASPHU0GF);\n",
    "    SKSbetath = auxevol_gfs + (grid_size*BETASPHU1GF);\n",
    "    SKSbetaph = auxevol_gfs + (grid_size*BETASPHU2GF);\n",
    "    SKSgrr = auxevol_gfs + (grid_size*GAMMASPHDD00GF);\n",
    "    SKSgrth = auxevol_gfs + (grid_size*GAMMASPHDD01GF);\n",
    "    SKSgrph = auxevol_gfs + (grid_size*GAMMASPHDD02GF);\n",
    "    SKSgthth = auxevol_gfs + (grid_size*GAMMASPHDD11GF);\n",
    "    SKSgthph = auxevol_gfs + (grid_size*GAMMASPHDD12GF);\n",
    "    SKSgphph = auxevol_gfs + (grid_size*GAMMASPHDD22GF);\n",
    "\n",
    "    Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);\n",
    "    Shifted_Kerr_Schild_initial_metric_spherical(&params,xx,auxevol_gfs);\n",
    "    GiRaFFE_NRPy_initial_data(&params,xx,auxevol_gfs,evol_gfs);\n",
    "    if(atoi(argv[1])==0) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);GiRaFFEfood_ExactWald();}\n",
    "    else if(atoi(argv[1])==1) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);Shifted_Kerr_Schild_initial_metric_spherical(&params,xx,auxevol_gfs);GiRaFFEfood_SplitMonopole();}\n",
    "    else if(atoi(argv[1])==2) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);Shifted_Kerr_Schild_initial_metric_spherical(&params,xx,auxevol_gfs);GiRaFFEfood_MagnetoWald();}\n",
    "\n",
    "    int all_agree = 1;\n",
    "\n",
    "    for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++){\n",
    "        for(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++){\n",
    "            for(int i2=0;i2<Nxx_plus_2NGHOSTS2;i2++){\n",
    "                if(SDA(auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)])<10.0){\n",
    "                    printf(\"Quantity ValenciavU0 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\\n\",\n",
    "                           SDA(auxevol_gfs[IDX4S(ALPHAGF, i0,i1,i2)]*auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)]-auxevol_gfs[IDX4S(BETAU0GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)]),i0,i1,i2);\n",
    "                    all_agree=0;\n",
    "                }\n",
    "                if(SDA(auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)])<10.0){\n",
    "                    printf(\"Quantity ValenciavU1 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\\n\",\n",
    "                           SDA(auxevol_gfs[IDX4S(ALPHAGF, i0,i1,i2)]*auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)]-auxevol_gfs[IDX4S(BETAU1GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)]),i0,i1,i2);\n",
    "                    all_agree=0;\n",
    "                }\n",
    "                if(SDA(auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)])<10.0){\n",
    "                    printf(\"Quantity ValenciavU2 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\\n\",\n",
    "                           SDA(auxevol_gfs[IDX4S(ALPHAGF, i0,i1,i2)]*auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]-auxevol_gfs[IDX4S(BETAU2GF, i0,i1,i2)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)]),i0,i1,i2);\n",
    "                    all_agree=0;\n",
    "                }\n",
    "                //printf(\"NRPy: %.15e,%.15e,%.15e\\n\",auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)],auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)],auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]);\n",
    "                //printf(\"CCTK: %.15e,%.15e,%.15e\\n\",vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)],vel[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)]);\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "\n",
    "    // Shift the grid to compare A_x\n",
    "    //for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] += 0.5*dxx1;\n",
    "    //for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] += 0.5*dxx2;\n",
    "    //GiRaFFE_NRPy_initial_data(&params,xx,auxevol_gfs,evol_gfs);\n",
    "    //if(atoi(argv[1])==0) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);GiRaFFEfood_ExactWald();}\n",
    "    //else if(atoi(argv[1])==1) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);Shifted_Kerr_Schild_initial_metric_spherical(&params,xx,auxevol_gfs);GiRaFFEfood_SplitMonopole();}\n",
    "    //else if(atoi(argv[1])==2) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);Shifted_Kerr_Schild_initial_metric_spherical(&params,xx,auxevol_gfs);GiRaFFEfood_MagnetoWald();}\n",
    "\n",
    "    for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++){\n",
    "        for(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++){\n",
    "            for(int i2=0;i2<Nxx_plus_2NGHOSTS2;i2++){\n",
    "                if(SDA(evol_gfs[IDX4S(AD0GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)])<10.0){\n",
    "                    printf(\"Quantity AD0 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\\n\",\n",
    "                           SDA(evol_gfs[IDX4S(AD0GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)]),i0,i1,i2);\n",
    "                    all_agree=0;\n",
    "                }\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "\n",
    "    // Shift the grid to compare A_y\n",
    "    //for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] += 0.5*dxx0;\n",
    "    //for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] -= 0.5*dxx1;\n",
    "    //GiRaFFE_NRPy_initial_data(&params,xx,auxevol_gfs,evol_gfs);\n",
    "    //if(atoi(argv[1])==0) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);GiRaFFEfood_ExactWald();}\n",
    "    //else if(atoi(argv[1])==1) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);Shifted_Kerr_Schild_initial_metric_spherical(&params,xx,auxevol_gfs);GiRaFFEfood_SplitMonopole();}\n",
    "    //else if(atoi(argv[1])==2) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);Shifted_Kerr_Schild_initial_metric_spherical(&params,xx,auxevol_gfs);GiRaFFEfood_MagnetoWald();}\n",
    "\n",
    "    for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++){\n",
    "        for(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++){\n",
    "            for(int i2=0;i2<Nxx_plus_2NGHOSTS2;i2++){\n",
    "                if(SDA(evol_gfs[IDX4S(AD1GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)])<10.0){\n",
    "                    printf(\"Quantity AD1 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\\n\",\n",
    "                           SDA(evol_gfs[IDX4S(AD1GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)]),i0,i1,i2);\n",
    "                    printf(\"NRPy: %.15e\\n\",evol_gfs[IDX4S(AD1GF, i0,i1,i2)]);\n",
    "                    printf(\"CCTK: %.15e\\n\\n\",Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)]);\n",
    "                    all_agree=0;\n",
    "                }\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "\n",
    "    // Shift the grid to compare A_z\n",
    "    //for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] += 0.5*dxx1;\n",
    "    //for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] -= 0.5*dxx2;\n",
    "    //GiRaFFE_NRPy_initial_data(&params,xx,auxevol_gfs,evol_gfs);\n",
    "    //if(atoi(argv[1])==0) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);GiRaFFEfood_ExactWald();}\n",
    "    //else if(atoi(argv[1])==1) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);Shifted_Kerr_Schild_initial_metric_spherical(&params,xx,auxevol_gfs);GiRaFFEfood_SplitMonopole();}\n",
    "    //else if(atoi(argv[1])==2) {Shifted_Kerr_Schild_initial_metric(&params,xx,auxevol_gfs);Shifted_Kerr_Schild_initial_metric_spherical(&params,xx,auxevol_gfs);GiRaFFEfood_MagnetoWald();}\n",
    "\n",
    "    for(int i0=0;i0<Nxx_plus_2NGHOSTS0;i0++){\n",
    "        for(int i1=0;i1<Nxx_plus_2NGHOSTS1;i1++){\n",
    "            for(int i2=0;i2<Nxx_plus_2NGHOSTS2;i2++){\n",
    "                if(SDA(evol_gfs[IDX4S(AD2GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)])<10.0){\n",
    "                    printf(\"Quantity AD2 only agrees with the original GiRaFFE to %.2f digits at i0,i1,i2=%d,%d,%d!\\n\",\n",
    "                           SDA(evol_gfs[IDX4S(AD2GF, i0,i1,i2)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)]),i0,i1,i2);\n",
    "                    all_agree=0;\n",
    "                }\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "\n",
    "    //printf(\"NRPy: %.15e,%.15e,%.15e; i,j,k=%d,%d,%d\\n\",evol_gfs[IDX4S(AD0GF, i0,i1,i2)],evol_gfs[IDX4S(AD1GF, i0,i1,i2)],evol_gfs[IDX4S(AD2GF, i0,i1,i2)],i0,i1,i2);\n",
    "    //printf(\"CCTK: %.15e,%.15e,%.15e\\n\",Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,0)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,1)],Avec[CCTK_GFINDEX4D(cctkGH,i0,i1,i2,2)]);\n",
    "\n",
    "    if(all_agree) printf(\"All quantities agree at all points!\\n\");\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='compile_run'></a>\n",
    "\n",
    "## Step 2.a: Compile and run the code to validate the output \\[Back to [top](#toc)\\]\n",
    "\n",
    "$$\\label{compile_run}$$\n",
    "\n",
    "Finally, we can compile and run the code we have written. Once run, this code will output the level of agreement between the two codes and some information to help interpret those numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:10.778433Z",
     "iopub.status.busy": "2022-12-29T09:54:10.778106Z",
     "iopub.status.idle": "2022-12-29T09:54:12.142404Z",
     "shell.execute_reply": "2022-12-29T09:54:12.141659Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Now compiling, should take ~2 seconds...\n",
      "\n",
      "Finished in 1.1283152103424072 seconds.\n",
      "\n",
      "\n",
      "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./GiRaFFEfood_NRPy_unit_test 2`...\n",
      "(BENCH): Finished executing in 0.20 seconds.\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "print(\"Now compiling, should take ~2 seconds...\\n\")\n",
    "start = time.time()\n",
    "# cmd.C_compile(os.path.join(out_dir,\"GiRaFFEfood_NRPy_unit_test.C\"), os.path.join(out_dir,\"GiRaFFEfood_NRPy_unit_test\")additional_libraries=\"-lgsl -lgslcblas\")\n",
    "!g++ -Ofast -fopenmp -march=native -funroll-loops Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Curved_Spacetime_Tests/GiRaFFEfood_NRPy_unit_test.C -o Validation-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Curved_Spacetime_Tests/GiRaFFEfood_NRPy_unit_test -lstdc++ -lgsl -lgslcblas\n",
    "end = time.time()\n",
    "print(\"Finished in \"+str(end-start)+\" seconds.\\n\\n\")\n",
    "\n",
    "results_file = \"out_GiRaFFEfood_NRPy_test.txt\"\n",
    "\n",
    "# os.chdir(out_dir)\n",
    "os.chdir(out_dir)\n",
    "# cmd.Execute(os.path.join(\"GiRaFFEfood_NRPy_unit_test\"))\n",
    "if initial_data==\"ExactWald\":\n",
    "    cmd.Execute(\"GiRaFFEfood_NRPy_unit_test\",\"0\",results_file)\n",
    "elif initial_data==\"SplitMonopole\":\n",
    "    cmd.Execute(\"GiRaFFEfood_NRPy_unit_test\",\"1\",results_file)\n",
    "elif initial_data==\"MagnetosphericWald\":\n",
    "    cmd.Execute(\"GiRaFFEfood_NRPy_unit_test\",\"2\",results_file)\n",
    "os.chdir(os.path.join(\"../\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we add some emergency brakes so that if the output from the test isn't good, we throw an error to stop the notebook dead in its tracks. This way, our automatic testing infrastructure can let us know if something goes wrong. We will also print the output from the test for convenience's sake."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:12.144090Z",
     "iopub.status.busy": "2022-12-29T09:54:12.143970Z",
     "iopub.status.idle": "2022-12-29T09:54:12.148022Z",
     "shell.execute_reply": "2022-12-29T09:54:12.147342Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All quantities agree at all points!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "with open(os.path.join(out_dir,results_file),\"r\") as file:\n",
    "    output = file.readline()\n",
    "    print(output)\n",
    "    if output!=\"All quantities agree at all points!\\n\": # If this isn't the first line of this file, something went wrong!\n",
    "        sys.exit(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 4: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy.pdf](Tutorial-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy.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": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-12-29T09:54:12.150573Z",
     "iopub.status.busy": "2022-12-29T09:54:12.150348Z",
     "iopub.status.idle": "2022-12-29T09:54:16.172548Z",
     "shell.execute_reply": "2022-12-29T09:54:16.171997Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-\n",
      "    Curved_Spacetime_Tests.tex, and compiled LaTeX file to PDF file\n",
      "    Tutorial-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-\n",
      "    Curved_Spacetime_Tests.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy-Curved_Spacetime_Tests\")"
   ]
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}