{
 "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",
    "# Tutorial-IllinoisGRMHD: ID_converter_ILGRMHD ETKThorn\n",
    "\n",
    "## Authors: Leo Werneck & Zach Etienne\n",
    "\n",
    "<font color='red'>**This module is currently under development**</font>\n",
    "\n",
    "## This tutorial module generates the `ID_converter_ILGRMHD` ETK thorn files, compatible with the latest implementation of IllinoisGRMHD.\n",
    "\n",
    "### Required and recommended citations:\n",
    "\n",
    "* **(Required)** Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. ([arxiv:1501.07276](http://arxiv.org/abs/1501.07276)).\n",
    "* **(Required)** Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) ([astro-ph/0512420](https://arxiv.org/abs/astro-ph/0512420)).\n",
    "* **(Recommended)** Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 ([astro-ph/0210618](https://arxiv.org/abs/astro-ph/0210618))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This module is organized as follows\n",
    "\n",
    "0. [Step 0](#src_dir): **Source directory creation**\n",
    "1. [Step 1](#introduction): **Introduction**\n",
    "1. [Step 2](#convert_to_hydrobase__src): **`set_IllinoisGRMHD_metric_GRMHD_variables_based_on_HydroBase_and_ADMBase_variables.C`**\n",
    "1. [Step 3](#convert_to_hydrobase__param): **`param.ccl`**\n",
    "1. [Step 4](#convert_to_hydrobase__interface): **`interface.ccl`**\n",
    "1. [Step 5](#convert_to_hydrobase__schedule): **`schedule.ccl`**\n",
    "1. [Step 6](#convert_to_hydrobase__make): **`make.code.defn`**\n",
    "1. [Step n-1](#code_validation): **Code validation**\n",
    "1. [Step n](#latex_pdf_output): **Output this notebook to $\\LaTeX$-formatted PDF file**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='src_dir'></a>\n",
    "\n",
    "# Step 0: Source directory creation \\[Back to [top](#toc)\\]\n",
    "$$\\label{src_dir}$$\n",
    "\n",
    "We will now use the [cmdline_helper.py NRPy+ module](Tutorial-Tutorial-cmdline_helper.ipynb) to create the source directory within the `IllinoisGRMHD` NRPy+ directory if it does not exist yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 0: Creation of the IllinoisGRMHD source directory\n",
    "# Step 0a: Load up cmdline_helper and create the directory\n",
    "import os,sys\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",
    "import cmdline_helper as cmd\n",
    "IDcIGM_dir_path = os.path.join(\"..\",\"ID_converter_ILGRMHD\")\n",
    "cmd.mkdir(IDcIGM_dir_path)\n",
    "IDcIGM_src_dir_path = os.path.join(IDcIGM_dir_path,\"src\")\n",
    "cmd.mkdir(IDcIGM_src_dir_path)\n",
    "\n",
    "# Step 0b: Create the output file path\n",
    "outfile_path__ID_converter_ILGRMHD__source    = os.path.join(IDcIGM_src_dir_path,\"set_IllinoisGRMHD_metric_GRMHD_variables_based_on_HydroBase_and_ADMBase_variables.C\")\n",
    "outfile_path__ID_converter_ILGRMHD__make      = os.path.join(IDcIGM_src_dir_path,\"make.code.defn\")\n",
    "outfile_path__ID_converter_ILGRMHD__param     = os.path.join(IDcIGM_dir_path,\"param.ccl\")\n",
    "outfile_path__ID_converter_ILGRMHD__interface = os.path.join(IDcIGM_dir_path,\"interface.ccl\")\n",
    "outfile_path__ID_converter_ILGRMHD__schedule  = os.path.join(IDcIGM_dir_path,\"schedule.ccl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='introduction'></a>\n",
    "\n",
    "# Step 1: Introduction \\[Back to [top](#toc)\\]\n",
    "$$\\label{introduction}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convert_to_hydrobase__src'></a>\n",
    "\n",
    "# Step 2: `set_IllinoisGRMHD_metric_GRMHD_variables _based_on_HydroBase_and_ADMBase_variables.C` \\[Back to [top](#toc)\\]\n",
    "$$\\label{convert_to_hydrobase__src}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ../ID_converter_ILGRMHD/src/set_IllinoisGRMHD_metric_GRMHD_variables_based_on_HydroBase_and_ADMBase_variables.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outfile_path__ID_converter_ILGRMHD__source\n",
    "/********************************\n",
    " * CONVERT ET ID TO IllinoisGRMHD\n",
    " *\n",
    " * Written in 2014 by Zachariah B. Etienne\n",
    " *\n",
    " * Sets metric & MHD variables needed\n",
    " * by IllinoisGRMHD, converting from\n",
    " * HydroBase and ADMBase.\n",
    " ********************************/\n",
    "\n",
    "#include <stdio.h>\n",
    "#include <stdlib.h>\n",
    "#include <math.h>\n",
    "#include <sys/time.h>\n",
    "#include \"cctk.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "#include \"cctk_Arguments.h\"\n",
    "#include \"IllinoisGRMHD_headers.h\"\n",
    "\n",
    "extern \"C\" void set_IllinoisGRMHD_metric_GRMHD_variables_based_on_HydroBase_and_ADMBase_variables(CCTK_ARGUMENTS) {\n",
    "\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  if(rho_b_atm > 1e199) {\n",
    "    CCTK_VError(VERR_DEF_PARAMS, \"You MUST set rho_b_atm to some reasonable value in your param.ccl file.\\n\");\n",
    "  }\n",
    "\n",
    "  // Convert ADM variables (from ADMBase) to the BSSN-based variables expected by this routine.\n",
    "  IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(cctkGH,cctk_lsh,  gxx,gxy,gxz,gyy,gyz,gzz,alp,\n",
    "                                                                              gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,\n",
    "                                                                              gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,\n",
    "                                                                              phi_bssn,psi_bssn,lapm1);\n",
    "\n",
    "  /***************\n",
    "   * PPEOS Patch *\n",
    "   ***************/\n",
    "  eos_struct eos;\n",
    "  initialize_EOS_struct_from_input(eos);\n",
    "\n",
    "  if(pure_hydro_run) {\n",
    "#pragma omp parallel for\n",
    "    for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {\n",
    "          int index=CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "          Avec[CCTK_GFINDEX4D(cctkGH,i,j,k,0)]=0;\n",
    "          Avec[CCTK_GFINDEX4D(cctkGH,i,j,k,1)]=0;\n",
    "          Avec[CCTK_GFINDEX4D(cctkGH,i,j,k,2)]=0;\n",
    "          Aphi[index]=0;\n",
    "        }\n",
    "  }\n",
    "\n",
    "#pragma omp parallel for\n",
    "  for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) {\n",
    "        int index=CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "\n",
    "        rho_b[index] = rho[index];\n",
    "        P[index] = press[index];\n",
    "\n",
    "        /***************\n",
    "         * PPEOS Patch *\n",
    "         ***************\n",
    "         * We now verify that the initial data\n",
    "         * provided by the user is indeed \"cold\",\n",
    "         * i.e. it contains no Thermal part and\n",
    "         * P = P_cold.\n",
    "         */\n",
    "        /* Compute P_cold */\n",
    "        const int polytropic_index = find_polytropic_K_and_Gamma_index(eos, rho_b[index]);\n",
    "        const double K_poly     = eos.K_ppoly_tab[polytropic_index];\n",
    "        const double Gamma_poly = eos.Gamma_ppoly_tab[polytropic_index];\n",
    "        const double P_cold     = K_poly*pow(rho_b[index],Gamma_poly);\n",
    "\n",
    "        /* Compare P and P_cold */\n",
    "        double P_rel_error = fabs(P[index] - P_cold)/P[index];\n",
    "        if( rho_b[index] > rho_b_atm && P_rel_error > 1e-2 ) {\n",
    "          const double Gamma_poly_local = log(P[index]/K_poly) / log(rho_b[index]);\n",
    "          /* Determine the value of Gamma_poly_local associated with P[index] */\n",
    "          CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,\n",
    "\"Expected a PP EOS with local Gamma_poly = %.15e, but found a point such that Gamma_poly_local = %.15e.\\n\",\n",
    "                     Gamma_poly, Gamma_poly_local);\n",
    "          CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,\n",
    "\"{rho_b; rho_b_atm; P; P_cold; P_rel_Error} = %.10e %e %.10e %.10e %e\\n\",\n",
    "                      rho_b[index], rho_b_atm, P[index],P_cold,P_rel_error);\n",
    "        }\n",
    "\n",
    "        Ax[index] = Avec[CCTK_GFINDEX4D(cctkGH,i,j,k,0)];\n",
    "        Ay[index] = Avec[CCTK_GFINDEX4D(cctkGH,i,j,k,1)];\n",
    "        Az[index] = Avec[CCTK_GFINDEX4D(cctkGH,i,j,k,2)];\n",
    "        psi6phi[index] = Aphi[index];\n",
    "\n",
    "        double ETvx = vel[CCTK_GFINDEX4D(cctkGH,i,j,k,0)];\n",
    "        double ETvy = vel[CCTK_GFINDEX4D(cctkGH,i,j,k,1)];\n",
    "        double ETvz = vel[CCTK_GFINDEX4D(cctkGH,i,j,k,2)];\n",
    "\n",
    "        // IllinoisGRMHD defines v^i = u^i/u^0.\n",
    "\n",
    "        // Meanwhile, the ET/HydroBase formalism, called the Valencia\n",
    "        // formalism, splits the 4 velocity into a purely spatial part\n",
    "        // and a part that is normal to the spatial hypersurface:\n",
    "        // u^a = G (n^a + U^a), (Eq. 14 of arXiv:1304.5544; G=W, U^a=v^a)\n",
    "        // where n^a is the unit normal vector to the spatial hypersurface,\n",
    "        // n_a = {-\\alpha,0,0,0}, and U^a is the purely spatial part, which\n",
    "        // is defined in HydroBase as the vel[] vector gridfunction.\n",
    "        // Then u^a n_a = - \\alpha u^0 = G n^a n_a = -G, and\n",
    "        // of course \\alpha u^0 = 1/sqrt(1+γ^ij u_i u_j) = \\Gamma,\n",
    "        // the standard Lorentz factor.\n",
    "\n",
    "        // Note that n^i = - \\beta^i / \\alpha, so\n",
    "        // u^a = \\Gamma (n^a + U^a)\n",
    "        // -> u^i = \\Gamma ( U^i - \\beta^i / \\alpha )\n",
    "        // which implies\n",
    "        // v^i = u^i/u^0\n",
    "        //     = \\Gamma/u^0 ( U^i - \\beta^i / \\alpha ) <- \\Gamma = \\alpha u^0\n",
    "        //     = \\alpha ( U^i - \\beta^i / \\alpha )\n",
    "        //     = \\alpha U^i - \\beta^i\n",
    "\n",
    "        vx[index] = alp[index]*ETvx - betax[index];\n",
    "        vy[index] = alp[index]*ETvy - betay[index];\n",
    "        vz[index] = alp[index]*ETvz - betaz[index];\n",
    "\n",
    "      }\n",
    "\n",
    "  // Neat feature for debugging: Add a roundoff-error perturbation\n",
    "  //    to the initial data.\n",
    "  // Set random_pert variable to ~1e-14 for a random 15th digit\n",
    "  //    perturbation.\n",
    "  srand(random_seed); // Use srand() as rand() is thread-safe.\n",
    "  for(int k=0;k<cctk_lsh[2];k++)\n",
    "    for(int j=0;j<cctk_lsh[1];j++)\n",
    "      for(int i=0;i<cctk_lsh[0];i++) {\n",
    "        int index=CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "        double pert = (random_pert*(double)rand() / RAND_MAX);\n",
    "        double one_plus_pert=(1.0+pert);\n",
    "        rho_b[index]*=one_plus_pert;\n",
    "        vx[index]*=one_plus_pert;\n",
    "        vy[index]*=one_plus_pert;\n",
    "        vz[index]*=one_plus_pert;\n",
    "\n",
    "        psi6phi[index]*=one_plus_pert;\n",
    "        Ax[index]*=one_plus_pert;\n",
    "        Ay[index]*=one_plus_pert;\n",
    "        Az[index]*=one_plus_pert;\n",
    "      }\n",
    "\n",
    "  // Next compute B & B_stagger from A_i. Note that this routine also depends on\n",
    "  //   the psi_bssn[] gridfunction being set to exp(phi).\n",
    "\n",
    "  double dxi = 1.0/CCTK_DELTA_SPACE(0);\n",
    "  double dyi = 1.0/CCTK_DELTA_SPACE(1);\n",
    "  double dzi = 1.0/CCTK_DELTA_SPACE(2);\n",
    "\n",
    "#pragma omp parallel for\n",
    "  for(int k=0;k<cctk_lsh[2];k++)\n",
    "    for(int j=0;j<cctk_lsh[1];j++)\n",
    "      for(int i=0;i<cctk_lsh[0];i++) {\n",
    "        // Look Mom, no if() statements!\n",
    "        int shiftedim1 = (i-1)*(i!=0); // This way, i=0 yields shiftedim1=0 and shiftedi=1, used below for our COPY boundary condition.\n",
    "        int shiftedi   = shiftedim1+1;\n",
    "\n",
    "        int shiftedjm1 = (j-1)*(j!=0);\n",
    "        int shiftedj   = shiftedjm1+1;\n",
    "\n",
    "        int shiftedkm1 = (k-1)*(k!=0);\n",
    "        int shiftedk   = shiftedkm1+1;\n",
    "\n",
    "        int index,indexim1,indexjm1,indexkm1;\n",
    "\n",
    "        int actual_index = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "\n",
    "        double Psi = psi_bssn[actual_index];\n",
    "        double Psim3 = 1.0/(Psi*Psi*Psi);\n",
    "\n",
    "        // For the lower boundaries, the following applies a \"copy\"\n",
    "        //    boundary condition on Bi_stagger where needed.\n",
    "        //    E.g., Bx_stagger(i,jmin,k) = Bx_stagger(i,jmin+1,k)\n",
    "        //    We find the copy BC works better than extrapolation.\n",
    "        // For the upper boundaries, we do the following copy:\n",
    "        //    E.g., Psi(imax+1,j,k)=Psi(imax,j,k)\n",
    "        /**************/\n",
    "        /* Bx_stagger */\n",
    "        /**************/\n",
    "\n",
    "        index    = CCTK_GFINDEX3D(cctkGH,i,shiftedj,shiftedk);\n",
    "        indexjm1 = CCTK_GFINDEX3D(cctkGH,i,shiftedjm1,shiftedk);\n",
    "        indexkm1 = CCTK_GFINDEX3D(cctkGH,i,shiftedj,shiftedkm1);\n",
    "        // Set Bx_stagger = \\partial_y A_z - partial_z A_y\n",
    "        // \"Grid\" Ax(i,j,k) is actually Ax(i,j+1/2,k+1/2)\n",
    "        // \"Grid\" Ay(i,j,k) is actually Ay(i+1/2,j,k+1/2)\n",
    "        // \"Grid\" Az(i,j,k) is actually Ay(i+1/2,j+1/2,k)\n",
    "        // Therefore, the 2nd order derivative \\partial_z A_y at (i+1/2,j,k) is:\n",
    "        //          [\"Grid\" Ay(i,j,k) - \"Grid\" Ay(i,j,k-1)]/dZ\n",
    "        Bx_stagger[actual_index] = (Az[index]-Az[indexjm1])*dyi - (Ay[index]-Ay[indexkm1])*dzi;\n",
    "\n",
    "        // Now multiply Bx and Bx_stagger by 1/sqrt(gamma(i+1/2,j,k)]) = 1/sqrt(1/2 [gamma + gamma_ip1]) = exp(-6 x 1/2 [phi + phi_ip1] )\n",
    "        int imax_minus_i = (cctk_lsh[0]-1)-i;\n",
    "        int indexip1jk = CCTK_GFINDEX3D(cctkGH,i + ( (imax_minus_i > 0) - (0 > imax_minus_i) ),j,k);\n",
    "        double Psi_ip1 = psi_bssn[indexip1jk];\n",
    "        Bx_stagger[actual_index] *= Psim3/(Psi_ip1*Psi_ip1*Psi_ip1);\n",
    "\n",
    "        /**************/\n",
    "        /* By_stagger */\n",
    "        /**************/\n",
    "\n",
    "        index    = CCTK_GFINDEX3D(cctkGH,shiftedi,j,shiftedk);\n",
    "        indexim1 = CCTK_GFINDEX3D(cctkGH,shiftedim1,j,shiftedk);\n",
    "        indexkm1 = CCTK_GFINDEX3D(cctkGH,shiftedi,j,shiftedkm1);\n",
    "        // Set By_stagger = \\partial_z A_x - \\partial_x A_z\n",
    "        By_stagger[actual_index] = (Ax[index]-Ax[indexkm1])*dzi - (Az[index]-Az[indexim1])*dxi;\n",
    "\n",
    "        // Now multiply By and By_stagger by 1/sqrt(gamma(i,j+1/2,k)]) = 1/sqrt(1/2 [gamma + gamma_jp1]) = exp(-6 x 1/2 [phi + phi_jp1] )\n",
    "        int jmax_minus_j = (cctk_lsh[1]-1)-j;\n",
    "        int indexijp1k = CCTK_GFINDEX3D(cctkGH,i,j + ( (jmax_minus_j > 0) - (0 > jmax_minus_j) ),k);\n",
    "        double Psi_jp1 = psi_bssn[indexijp1k];\n",
    "        By_stagger[actual_index] *= Psim3/(Psi_jp1*Psi_jp1*Psi_jp1);\n",
    "\n",
    "\n",
    "        /**************/\n",
    "        /* Bz_stagger */\n",
    "        /**************/\n",
    "\n",
    "        index    = CCTK_GFINDEX3D(cctkGH,shiftedi,shiftedj,k);\n",
    "        indexim1 = CCTK_GFINDEX3D(cctkGH,shiftedim1,shiftedj,k);\n",
    "        indexjm1 = CCTK_GFINDEX3D(cctkGH,shiftedi,shiftedjm1,k);\n",
    "        // Set Bz_stagger = \\partial_x A_y - \\partial_y A_x\n",
    "        Bz_stagger[actual_index] = (Ay[index]-Ay[indexim1])*dxi - (Ax[index]-Ax[indexjm1])*dyi;\n",
    "\n",
    "        // Now multiply Bz_stagger by 1/sqrt(gamma(i,j,k+1/2)]) = 1/sqrt(1/2 [gamma + gamma_kp1]) = exp(-6 x 1/2 [phi + phi_kp1] )\n",
    "        int kmax_minus_k = (cctk_lsh[2]-1)-k;\n",
    "        int indexijkp1 = CCTK_GFINDEX3D(cctkGH,i,j,k + ( (kmax_minus_k > 0) - (0 > kmax_minus_k) ));\n",
    "        double Psi_kp1 = psi_bssn[indexijkp1];\n",
    "        Bz_stagger[actual_index] *= Psim3/(Psi_kp1*Psi_kp1*Psi_kp1);\n",
    "\n",
    "      }\n",
    "\n",
    "#pragma omp parallel for\n",
    "  for(int k=0;k<cctk_lsh[2];k++)\n",
    "    for(int j=0;j<cctk_lsh[1];j++)\n",
    "      for(int i=0;i<cctk_lsh[0];i++) {\n",
    "        // Look Mom, no if() statements!\n",
    "        int shiftedim1 = (i-1)*(i!=0); // This way, i=0 yields shiftedim1=0 and shiftedi=1, used below for our COPY boundary condition.\n",
    "        int shiftedi   = shiftedim1+1;\n",
    "\n",
    "        int shiftedjm1 = (j-1)*(j!=0);\n",
    "        int shiftedj   = shiftedjm1+1;\n",
    "\n",
    "        int shiftedkm1 = (k-1)*(k!=0);\n",
    "        int shiftedk   = shiftedkm1+1;\n",
    "\n",
    "        int index,indexim1,indexjm1,indexkm1;\n",
    "\n",
    "        int actual_index = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "\n",
    "        // For the lower boundaries, the following applies a \"copy\"\n",
    "        //    boundary condition on Bi and Bi_stagger where needed.\n",
    "        //    E.g., Bx(imin,j,k) = Bx(imin+1,j,k)\n",
    "        //    We find the copy BC works better than extrapolation.\n",
    "        /******/\n",
    "        /* Bx */\n",
    "        /******/\n",
    "        index = CCTK_GFINDEX3D(cctkGH,shiftedi,j,k);\n",
    "        indexim1 = CCTK_GFINDEX3D(cctkGH,shiftedim1,j,k);\n",
    "        // Set Bx = 0.5 ( Bx_stagger + Bx_stagger_im1 )\n",
    "        // \"Grid\" Bx_stagger(i,j,k) is actually Bx_stagger(i+1/2,j,k)\n",
    "        Bx[actual_index] = 0.5 * ( Bx_stagger[index] + Bx_stagger[indexim1] );\n",
    "\n",
    "        /******/\n",
    "        /* By */\n",
    "        /******/\n",
    "        index = CCTK_GFINDEX3D(cctkGH,i,shiftedj,k);\n",
    "        indexjm1 = CCTK_GFINDEX3D(cctkGH,i,shiftedjm1,k);\n",
    "        // Set By = 0.5 ( By_stagger + By_stagger_im1 )\n",
    "        // \"Grid\" By_stagger(i,j,k) is actually By_stagger(i,j+1/2,k)\n",
    "        By[actual_index] = 0.5 * ( By_stagger[index] + By_stagger[indexjm1] );\n",
    "\n",
    "        /******/\n",
    "        /* Bz */\n",
    "        /******/\n",
    "        index = CCTK_GFINDEX3D(cctkGH,i,j,shiftedk);\n",
    "        indexkm1 = CCTK_GFINDEX3D(cctkGH,i,j,shiftedkm1);\n",
    "        // Set Bz = 0.5 ( Bz_stagger + Bz_stagger_im1 )\n",
    "        // \"Grid\" Bz_stagger(i,j,k) is actually Bz_stagger(i,j+1/2,k)\n",
    "        Bz[actual_index] = 0.5 * ( Bz_stagger[index] + Bz_stagger[indexkm1] );\n",
    "      }\n",
    "\n",
    "  // Finally, enforce limits on primitives & compute conservative variables.\n",
    "#pragma omp parallel for\n",
    "  for(int k=0;k<cctk_lsh[2];k++)\n",
    "    for(int j=0;j<cctk_lsh[1];j++)\n",
    "      for(int i=0;i<cctk_lsh[0];i++) {\n",
    "        static const int zero_int=0;\n",
    "        int index = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "\n",
    "        int ww;\n",
    "\n",
    "        double PRIMS[MAXNUMVARS];\n",
    "        ww=0;\n",
    "        PRIMS[ww] = rho_b[index]; ww++;\n",
    "        PRIMS[ww] = P[index];     ww++;\n",
    "        PRIMS[ww] = vx[index];    ww++;\n",
    "        PRIMS[ww] = vy[index];    ww++;\n",
    "        PRIMS[ww] = vz[index];    ww++;\n",
    "        PRIMS[ww] = Bx[index];    ww++;\n",
    "        PRIMS[ww] = By[index];    ww++;\n",
    "        PRIMS[ww] = Bz[index];    ww++;\n",
    "\n",
    "        double METRIC[NUMVARS_FOR_METRIC],dummy=0;\n",
    "        ww=0;\n",
    "        // FIXME: NECESSARY?\n",
    "        //psi_bssn[index] = exp(phi[index]);\n",
    "        METRIC[ww] = phi_bssn[index];ww++;\n",
    "        METRIC[ww] = dummy;          ww++; // Don't need to set psi.\n",
    "        METRIC[ww] = gtxx[index];    ww++;\n",
    "        METRIC[ww] = gtxy[index];    ww++;\n",
    "        METRIC[ww] = gtxz[index];    ww++;\n",
    "        METRIC[ww] = gtyy[index];    ww++;\n",
    "        METRIC[ww] = gtyz[index];    ww++;\n",
    "        METRIC[ww] = gtzz[index];    ww++;\n",
    "        METRIC[ww] = lapm1[index];   ww++;\n",
    "        METRIC[ww] = betax[index];   ww++;\n",
    "        METRIC[ww] = betay[index];   ww++;\n",
    "        METRIC[ww] = betaz[index];   ww++;\n",
    "        METRIC[ww] = gtupxx[index];  ww++;\n",
    "        METRIC[ww] = gtupyy[index];  ww++;\n",
    "        METRIC[ww] = gtupzz[index];  ww++;\n",
    "        METRIC[ww] = gtupxy[index];  ww++;\n",
    "        METRIC[ww] = gtupxz[index];  ww++;\n",
    "        METRIC[ww] = gtupyz[index];  ww++;\n",
    "\n",
    "        double CONSERVS[NUM_CONSERVS] = {0,0,0,0,0};\n",
    "        double g4dn[4][4];\n",
    "        double g4up[4][4];\n",
    "        double TUPMUNU[10],TDNMUNU[10];\n",
    "\n",
    "        struct output_stats stats; stats.failure_checker=0;\n",
    "        IllinoisGRMHD_enforce_limits_on_primitives_and_recompute_conservs(zero_int,PRIMS,stats,eos,\n",
    "                                                                          METRIC,g4dn,g4up,TUPMUNU,TDNMUNU,CONSERVS);\n",
    "        rho_b[index] = PRIMS[RHOB];\n",
    "        P[index]     = PRIMS[PRESSURE];\n",
    "        vx[index]    = PRIMS[VX];\n",
    "        vy[index]    = PRIMS[VY];\n",
    "        vz[index]    = PRIMS[VZ];\n",
    "\n",
    "        rho_star[index] = CONSERVS[RHOSTAR];\n",
    "        mhd_st_x[index] = CONSERVS[STILDEX];\n",
    "        mhd_st_y[index] = CONSERVS[STILDEY];\n",
    "        mhd_st_z[index] = CONSERVS[STILDEZ];\n",
    "        tau[index]      = CONSERVS[TAUENERGY];\n",
    "\n",
    "        if(update_Tmunu) {\n",
    "          ww=0;\n",
    "          eTtt[index] = TDNMUNU[ww]; ww++;\n",
    "          eTtx[index] = TDNMUNU[ww]; ww++;\n",
    "          eTty[index] = TDNMUNU[ww]; ww++;\n",
    "          eTtz[index] = TDNMUNU[ww]; ww++;\n",
    "          eTxx[index] = TDNMUNU[ww]; ww++;\n",
    "          eTxy[index] = TDNMUNU[ww]; ww++;\n",
    "          eTxz[index] = TDNMUNU[ww]; ww++;\n",
    "          eTyy[index] = TDNMUNU[ww]; ww++;\n",
    "          eTyz[index] = TDNMUNU[ww]; ww++;\n",
    "          eTzz[index] = TDNMUNU[ww];\n",
    "        }\n",
    "      }\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convert_to_hydrobase__param'></a>\n",
    "\n",
    "# Step 3: `param.ccl` \\[Back to [top](#toc)\\]\n",
    "$$\\label{convert_to_hydrobase__param}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ../ID_converter_ILGRMHD/param.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outfile_path__ID_converter_ILGRMHD__param\n",
    "# Parameter definitions for thorn ID_converter_ILGRMHD\n",
    "\n",
    "shares: IllinoisGRMHD\n",
    "USES INT neos\n",
    "USES REAL K_ppoly_tab0\n",
    "USES REAL rho_ppoly_tab_in[10]\n",
    "USES REAL Gamma_ppoly_tab_in[10]\n",
    "USES REAL rho_b_max\n",
    "USES REAL rho_b_atm\n",
    "USES REAL tau_atm\n",
    "USES REAL Sym_Bz\n",
    "USES REAL GAMMA_SPEED_LIMIT\n",
    "USES REAL Psi6threshold\n",
    "USES BOOLEAN update_Tmunu\n",
    "\n",
    "private:\n",
    "\n",
    "INT random_seed \"Random seed for random, generally roundoff-level perturbation on initial data. Seeds srand(), and rand() is used for the RNG.\"\n",
    "{\n",
    " 0:99999999 :: \"Anything unsigned goes.\"\n",
    "} 0\n",
    "\n",
    "REAL random_pert \"Random perturbation atop data\"\n",
    "{\n",
    " *:* :: \"Anything goes.\"\n",
    "} 0\n",
    "\n",
    "BOOLEAN pure_hydro_run \"Set the vector potential and corresponding EM gauge quantity to zero\"\n",
    "{\n",
    "} \"no\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convert_to_hydrobase__interface'></a>\n",
    "\n",
    "# Step 4: `interface.ccl` \\[Back to [top](#toc)\\]\n",
    "$$\\label{convert_to_hydrobase__interface}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ../ID_converter_ILGRMHD/interface.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outfile_path__ID_converter_ILGRMHD__interface\n",
    "# Interface definition for thorn ID_converter_ILGRMHD\n",
    "\n",
    "implements: ID_converter_ILGRMHD\n",
    "\n",
    "inherits: ADMBase, Boundary, SpaceMask, Tmunubase, HydroBase, grid, IllinoisGRMHD\n",
    "\n",
    "uses include header: IllinoisGRMHD_headers.h\n",
    "USES INCLUDE: Symmetry.h\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convert_to_hydrobase__schedule'></a>\n",
    "\n",
    "# Step 5: `schedule.ccl` \\[Back to [top](#toc)\\]\n",
    "$$\\label{convert_to_hydrobase__schedule}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ../ID_converter_ILGRMHD/schedule.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outfile_path__ID_converter_ILGRMHD__schedule\n",
    "# Schedule definitions for thorn ID_converter_ILGRMHD\n",
    "\n",
    "schedule group IllinoisGRMHD_ID_Converter at CCTK_INITIAL after HydroBase_Initial before Convert_to_HydroBase\n",
    "{\n",
    "} \"Translate ET-generated, HydroBase-compatible initial data and convert into variables used by IllinoisGRMHD\"\n",
    "\n",
    "schedule set_IllinoisGRMHD_metric_GRMHD_variables_based_on_HydroBase_and_ADMBase_variables IN IllinoisGRMHD_ID_Converter as first_initialdata before TOV_Initial_Data\n",
    "{\n",
    "   LANG:       C\n",
    "   OPTIONS:    LOCAL\n",
    "   # What the heck, let's synchronize everything!\n",
    "   SYNC: IllinoisGRMHD::grmhd_primitives_Bi, IllinoisGRMHD::grmhd_primitives_Bi_stagger, IllinoisGRMHD::grmhd_primitives_allbutBi, IllinoisGRMHD::em_Ax,IllinoisGRMHD::em_Ay,IllinoisGRMHD::em_Az,IllinoisGRMHD::em_psi6phi,IllinoisGRMHD::grmhd_conservatives,IllinoisGRMHD::BSSN_quantities,ADMBase::metric,ADMBase::lapse,ADMBase::shift,ADMBase::curv\n",
    "} \"Convert HydroBase initial data (ID) to ID that IllinoisGRMHD can read.\"\n",
    "\n",
    "schedule IllinoisGRMHD_InitSymBound IN IllinoisGRMHD_ID_Converter as third_initialdata after second_initialdata\n",
    "{\n",
    "  SYNC: IllinoisGRMHD::grmhd_conservatives,IllinoisGRMHD::em_Ax,IllinoisGRMHD::em_Ay,IllinoisGRMHD::em_Az,IllinoisGRMHD::em_psi6phi\n",
    "  LANG: C\n",
    "} \"Schedule symmetries -- Actually just a placeholder function to ensure prolongation / processor syncs are done BEFORE the primitives solver.\"\n",
    "\n",
    "schedule IllinoisGRMHD_compute_B_and_Bstagger_from_A IN IllinoisGRMHD_ID_Converter as fourth_initialdata after third_initialdata\n",
    "{\n",
    "  SYNC: IllinoisGRMHD::grmhd_primitives_Bi, IllinoisGRMHD::grmhd_primitives_Bi_stagger\n",
    "  LANG: C\n",
    "} \"Compute B and B_stagger from A\"\n",
    "\n",
    "schedule IllinoisGRMHD_conserv_to_prims IN IllinoisGRMHD_ID_Converter as fifth_initialdata after fourth_initialdata\n",
    "{\n",
    "  LANG: C\n",
    "} \"Compute primitive variables from conservatives. This is non-trivial, requiring a Newton-Raphson root-finder.\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convert_to_hydrobase__make'></a>\n",
    "\n",
    "# Step 6: `make.code.defn` \\[Back to [top](#toc)\\]\n",
    "$$\\label{convert_to_hydrobase__make}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting ../ID_converter_ILGRMHD/src/make.code.defn\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outfile_path__ID_converter_ILGRMHD__make\n",
    "# Main make.code.defn file for thorn ID_converter_ILGRMHD\n",
    "\n",
    "# Source files in this directory\n",
    "SRCS =  set_IllinoisGRMHD_metric_GRMHD_variables_based_on_HydroBase_and_ADMBase_variables.C\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step n-1: Code validation \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "First, we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Verify if the code generated by this tutorial module\n",
    "# # matches the original IllinoisGRMHD source code\n",
    "\n",
    "# # First download the original IllinoisGRMHD source code\n",
    "# import urllib\n",
    "# from os import path\n",
    "\n",
    "# original_IGM_file_url  = \"https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/A_i_rhs_no_gauge_terms.C\"\n",
    "# original_IGM_file_name = \"A_i_rhs_no_gauge_terms-original.C\"\n",
    "# original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)\n",
    "\n",
    "# # Then download the original IllinoisGRMHD 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_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
    "#     # Write down the file the original IllinoisGRMHD source code\n",
    "#     with open(original_IGM_file_path,\"w\") as file:\n",
    "#         file.write(original_IGM_file_code)\n",
    "# except:\n",
    "#     try:\n",
    "#         original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
    "#         # Write down the file the original IllinoisGRMHD source code\n",
    "#         with open(original_IGM_file_path,\"w\") as file:\n",
    "#             file.write(original_IGM_file_code)\n",
    "#     except:\n",
    "#         # If all else fails, hope wget does the job\n",
    "#         !wget -O $original_IGM_file_path $original_IGM_file_url\n",
    "\n",
    "# # Perform validation\n",
    "# Validation__A_i_rhs_no_gauge_terms__C  = !diff $original_IGM_file_path $outfile_path__A_i_rhs_no_gauge_terms__C\n",
    "\n",
    "# if Validation__A_i_rhs_no_gauge_terms__C == []:\n",
    "#     # If the validation passes, we do not need to store the original IGM source code file\n",
    "#     !rm $original_IGM_file_path\n",
    "#     print(\"Validation test for A_i_rhs_no_gauge_terms.C: PASSED!\")\n",
    "# else:\n",
    "#     # If the validation fails, we keep the original IGM source code file\n",
    "#     print(\"Validation test for A_i_rhs_no_gauge_terms.C: FAILED!\")\n",
    "#     # We also print out the difference between the code generated\n",
    "#     # in this tutorial module and the original IGM source code\n",
    "#     print(\"Diff:\")\n",
    "#     for diff_line in Validation__A_i_rhs_no_gauge_terms__C:\n",
    "#         print(diff_line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step n: 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-IllinoisGRMHD__A_i_rhs_no_gauge_terms.pdf](Tutorial-IllinoisGRMHD__A_i_rhs_no_gauge_terms.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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "latex_nrpy_style_path = os.path.join(nrpy_dir_path,\"latex_nrpy_style.tplx\")\n",
    "#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__A_i_rhs_no_gauge_terms.ipynb\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__A_i_rhs_no_gauge_terms.tex\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__A_i_rhs_no_gauge_terms.tex\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__A_i_rhs_no_gauge_terms.tex\n",
    "!rm -f Tut*.out Tut*.aux Tut*.log"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}