{
 "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: Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C\n",
    "\n",
    "## Authors: Leo Werneck & Zach Etienne\n",
    "\n",
    "<font color='red'>**This module is currently under development**</font>\n",
    "\n",
    "## This tutorial module constructs the right-hand side of the evolution equations of $\\left[\\sqrt{\\gamma}\\Phi\\right]$ and adds gauge terms to the right-hand side of the evolution equations of $A_{i}$.\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](#lorenz_psi6phi_rhs__add_gauge_terms_to_a_i_rhs__c): **`Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C`**\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: Add NRPy's directory to the path\n",
    "# https://stackoverflow.com/questions/16780014/import-file-from-parent-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",
    "# Step 0b: Load up cmdline_helper and create the directory\n",
    "import cmdline_helper as cmd\n",
    "IGM_src_dir_path = os.path.join(\"..\",\"src\")\n",
    "cmd.mkdir(IGM_src_dir_path)\n",
    "\n",
    "# Step 0c: Create the output file path\n",
    "outfile_path__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs__C = os.path.join(IGM_src_dir_path,\"Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C\")"
   ]
  },
  {
   "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='lorenz_psi6phi_rhs__add_gauge_terms_to_a_i_rhs__c'></a>\n",
    "\n",
    "# Step 2: `Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C` \\[Back to [top](#toc)\\]\n",
    "$$\\label{lorenz_psi6phi_rhs__add_gauge_terms_to_a_i_rhs__c}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing ../src/Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outfile_path__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs__C\n",
    "static inline CCTK_REAL avg(CCTK_REAL f[PLUS2+1][PLUS2+1][PLUS2+1],int imin,int imax, int jmin,int jmax, int kmin,int kmax);\n",
    "\n",
    "static void Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,CCTK_REAL *dX,CCTK_REAL **in_vars,CCTK_REAL *psi6phi,\n",
    "                                                           CCTK_REAL *shiftx_iphjphkph,CCTK_REAL *shifty_iphjphkph,CCTK_REAL *shiftz_iphjphkph,\n",
    "                                                           CCTK_REAL *alpha_iphjphkph,CCTK_REAL *alpha_Phi_minus_betaj_A_j_iphjphkph,CCTK_REAL *alpha_sqrtg_Ax_interp,\n",
    "                                                           CCTK_REAL *alpha_sqrtg_Ay_interp,CCTK_REAL *alpha_sqrtg_Az_interp,\n",
    "                                                           CCTK_REAL *psi6phi_rhs,CCTK_REAL *Ax_rhs,CCTK_REAL *Ay_rhs,CCTK_REAL *Az_rhs) {\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  /* Compute \\partial_t psi6phi = -\\partial_i (  \\alpha psi^6 A^i - psi6phi \\beta^i)\n",
    "   *    (Eq 13 of http://arxiv.org/pdf/1110.4633.pdf), using Lorenz gauge.\n",
    "   * Note that the RHS consists of a shift advection term on psi6phi and\n",
    "   *    a term depending on the vector potential.\n",
    "   * psi6phi is defined at (i+1/2,j+1/2,k+1/2), but instead of reconstructing\n",
    "   *    to compute the RHS of \\partial_t psi6phi, we instead use standard\n",
    "   *    interpolations.\n",
    "   */\n",
    "\n",
    "  CCTK_REAL dXm1=1.0/dX[0];\n",
    "  CCTK_REAL dYm1=1.0/dX[1];\n",
    "  CCTK_REAL dZm1=1.0/dX[2];\n",
    "\n",
    "  // The stencil here is {-1,1},{-1,1},{-1,1} for x,y,z directions, respectively.\n",
    "  //     Note that ALL input variables are defined at ALL gridpoints, so no\n",
    "  //     worries about ghostzones.\n",
    "#pragma omp parallel for\n",
    "  for(int k=1;k<cctk_lsh[2]-1;k++) for(int j=1;j<cctk_lsh[1]-1;j++) for(int i=1;i<cctk_lsh[0]-1;i++) {\n",
    "        int index=CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "        CCTK_REAL INTERP_VARS[MAXNUMINTERP][PLUS2+1][PLUS2+1][PLUS2+1];\n",
    "\n",
    "        // First compute \\partial_j \\alpha \\sqrt{\\gamma} A^j (RHS of \\partial_i psi6phi)\n",
    "        // FIXME: Would be much cheaper & easier to unstagger A_i, raise, then interpolate A^i.\n",
    "        //        However, we keep it this way to be completely compatible with the original\n",
    "        //        Illinois GRMHD thorn, called mhd_evolve.\n",
    "        //\n",
    "        //Step 1) j=x: Need to raise A_i, but to do that, we must have all variables at the same gridpoints:\n",
    "        // The goal is to compute \\partial_j (\\alpha \\sqrt{\\gamma} A^j) at (i+1/2,j+1/2,k+1/2)\n",
    "        //    We do this by first interpolating (RHS1x) = (\\alpha \\sqrt{\\gamma} A^x) at\n",
    "        //    (i,j+1/2,k+1/2)and (i+1,j+1/2,k+1/2), then taking \\partial_x (RHS1x) =\n",
    "        //    [ RHS1x(i+1,j+1/2,k+1/2) - RHS1x(i,j+1/2,k+1/2) ]/dX.\n",
    "        // First bring gup's, psi, and alpha to (i,j+1/2,k+1/2):\n",
    "        int num_vars_to_interp;\n",
    "        int vars_to_interpolate[MAXNUMINTERP] = {GUPXXI,GUPXYI,GUPXZI,GUPYYI,GUPYZI,GUPZZI,LAPM1I,PSII,SHIFTXI,SHIFTYI,SHIFTZI};\n",
    "        num_vars_to_interp = 11;\n",
    "        // We may set interp_limits to be more general than we need.\n",
    "        int interp_limits[6] = {-1,1,-1,1,-1,1}; SET_INDEX_ARRAYS_3DBLOCK(interp_limits);\n",
    "        //SET_INDEX_ARRAYS_3DBLOCK(interp_limits);\n",
    "        for(int ww=0;ww<num_vars_to_interp;ww++) {\n",
    "          int whichvar=vars_to_interpolate[ww];\n",
    "          // Read in variable at interp. stencil points from main memory, store in INTERP_VARS.\n",
    "          for(int kk=PLUS0;kk<=PLUS1;kk++) for(int jj=PLUS0;jj<=PLUS1;jj++) for(int ii=PLUS0;ii<=PLUS1;ii++) {\n",
    "                INTERP_VARS[whichvar][kk][jj][ii] = in_vars[whichvar][index_arr_3DB[kk][jj][ii]]; }\n",
    "        }\n",
    "\n",
    "        // Next set \\alpha at (i+1/2,j+1/2,k+1/2). Will come in handy when computing damping term later.\n",
    "        alpha_iphjphkph[index] = avg(INTERP_VARS[LAPM1I] , PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS1)+1.0;\n",
    "\n",
    "        //A_x needs a stencil s.t. interp_limits={0,1,-1,1,-1,1}:\n",
    "        for(int kk=MINUS1;kk<=PLUS1;kk++) for(int jj=MINUS1;jj<=PLUS1;jj++) for(int ii=PLUS0;ii<=PLUS1;ii++) {\n",
    "              INTERP_VARS[A_XI][kk][jj][ii] = in_vars[A_XI][index_arr_3DB[kk][jj][ii]]; }\n",
    "        //A_y needs a stencil s.t. interp_limits={-1,1,0,1,-1,1}:\n",
    "        for(int kk=MINUS1;kk<=PLUS1;kk++) for(int jj=PLUS0;jj<=PLUS1;jj++) for(int ii=MINUS1;ii<=PLUS1;ii++) {\n",
    "              INTERP_VARS[A_YI][kk][jj][ii] = in_vars[A_YI][index_arr_3DB[kk][jj][ii]]; }\n",
    "        //A_z needs a stencil s.t. interp_limits={-1,1,-1,1,0,1}:\n",
    "        for(int kk=PLUS0;kk<=PLUS1;kk++) for(int jj=MINUS1;jj<=PLUS1;jj++) for(int ii=MINUS1;ii<=PLUS1;ii++) {\n",
    "              INTERP_VARS[A_ZI][kk][jj][ii] = in_vars[A_ZI][index_arr_3DB[kk][jj][ii]]; }\n",
    "\n",
    "        // FIRST DO A^X TERM (interpolate to (i,j+1/2,k+1/2) )\n",
    "        // \\alpha \\sqrt{\\gamma} A^x = \\alpha psi^6 A^x (RHS of \\partial_i psi6phi)\n",
    "        // Note that gupij is \\tilde{\\gamma}^{ij}, so we need to multiply by \\psi^{-4}.\n",
    "        CCTK_REAL gupxx_jphkph = avg(INTERP_VARS[GUPXXI], PLUS0,PLUS0, PLUS0,PLUS1, PLUS0,PLUS1);\n",
    "        CCTK_REAL gupxy_jphkph = avg(INTERP_VARS[GUPXYI], PLUS0,PLUS0, PLUS0,PLUS1, PLUS0,PLUS1);\n",
    "        CCTK_REAL gupxz_jphkph = avg(INTERP_VARS[GUPXZI], PLUS0,PLUS0, PLUS0,PLUS1, PLUS0,PLUS1);\n",
    "\n",
    "        for(int kk=PLUS0;kk<=PLUS1;kk++) for(int jj=PLUS0;jj<=PLUS1;jj++) for(int ii=PLUS0;ii<=PLUS1;ii++) {\n",
    "              CCTK_REAL Psi2 = INTERP_VARS[PSII][kk][jj][ii]*INTERP_VARS[PSII][kk][jj][ii];\n",
    "              CCTK_REAL alpha = INTERP_VARS[LAPM1I][kk][jj][ii]+1.0;\n",
    "              INTERP_VARS[LAPSE_PSI2I][kk][jj][ii]=alpha*Psi2;\n",
    "              INTERP_VARS[LAPSE_OVER_PSI6I][kk][jj][ii]=alpha/(Psi2*Psi2*Psi2);\n",
    "            }\n",
    "\n",
    "        CCTK_REAL lapse_Psi2_jphkph = avg(INTERP_VARS[LAPSE_PSI2I], PLUS0,PLUS0, PLUS0,PLUS1, PLUS0,PLUS1);\n",
    "\n",
    "        CCTK_REAL A_x_jphkph   = avg(INTERP_VARS[A_XI], PLUS0,PLUS0, PLUS0,PLUS0, PLUS0,PLUS0); // @ (i,j+1/2,k+1/2)\n",
    "        CCTK_REAL A_y_jphkph   = avg(INTERP_VARS[A_YI],MINUS1,PLUS0, PLUS0,PLUS1, PLUS0,PLUS0); // @ (i+1/2,j,k+1/2)\n",
    "        CCTK_REAL A_z_jphkph   = avg(INTERP_VARS[A_ZI],MINUS1,PLUS0, PLUS0,PLUS0, PLUS0,PLUS1); // @ (i+1/2,j+1/2,k)\n",
    "\n",
    "        alpha_sqrtg_Ax_interp[index] = lapse_Psi2_jphkph*\n",
    "          ( gupxx_jphkph*A_x_jphkph + gupxy_jphkph*A_y_jphkph + gupxz_jphkph*A_z_jphkph );\n",
    "\n",
    "\n",
    "        // DO A^Y TERM (interpolate to (i+1/2,j,k+1/2) )\n",
    "        // \\alpha \\sqrt{\\gamma} A^y = \\alpha psi^6 A^y (RHS of \\partial_i psi6phi)\n",
    "        // Note that gupij is \\tilde{\\gamma}^{ij}, so we need to multiply by \\psi^{-4}.\n",
    "        CCTK_REAL gupxy_iphkph = avg(INTERP_VARS[GUPXYI], PLUS0,PLUS1, PLUS0,PLUS0, PLUS0,PLUS1);\n",
    "        CCTK_REAL gupyy_iphkph = avg(INTERP_VARS[GUPYYI], PLUS0,PLUS1, PLUS0,PLUS0, PLUS0,PLUS1);\n",
    "        CCTK_REAL gupyz_iphkph = avg(INTERP_VARS[GUPYZI], PLUS0,PLUS1, PLUS0,PLUS0, PLUS0,PLUS1);\n",
    "\n",
    "        CCTK_REAL lapse_Psi2_iphkph = avg(INTERP_VARS[LAPSE_PSI2I], PLUS0,PLUS1, PLUS0,PLUS0, PLUS0,PLUS1);\n",
    "        //CCTK_REAL lapse_iphkph = avg(INTERP_VARS[LAPM1I], PLUS0,PLUS1, PLUS0,PLUS0, PLUS0,PLUS1)+1.0;\n",
    "        //CCTK_REAL psi_iphkph   = avg(INTERP_VARS[PSII  ], PLUS0,PLUS1, PLUS0,PLUS0, PLUS0,PLUS1);\n",
    "\n",
    "        CCTK_REAL A_x_iphkph   = avg(INTERP_VARS[A_XI], PLUS0,PLUS1,MINUS1,PLUS0, PLUS0,PLUS0); // @ (i,j+1/2,k+1/2)\n",
    "        CCTK_REAL A_y_iphkph   = avg(INTERP_VARS[A_YI], PLUS0,PLUS0, PLUS0,PLUS0, PLUS0,PLUS0); // @ (i+1/2,j,k+1/2)\n",
    "        CCTK_REAL A_z_iphkph   = avg(INTERP_VARS[A_ZI], PLUS0,PLUS0,MINUS1,PLUS0, PLUS0,PLUS1); // @ (i+1/2,j+1/2,k)\n",
    "\n",
    "        alpha_sqrtg_Ay_interp[index] = lapse_Psi2_iphkph*\n",
    "          ( gupxy_iphkph*A_x_iphkph + gupyy_iphkph*A_y_iphkph + gupyz_iphkph*A_z_iphkph );\n",
    "\n",
    "        // DO A^Z TERM (interpolate to (i+1/2,j+1/2,k) )\n",
    "        // \\alpha \\sqrt{\\gamma} A^z = \\alpha psi^6 A^z (RHS of \\partial_i psi6phi)\n",
    "        // Note that gupij is \\tilde{\\gamma}^{ij}, so we need to multiply by \\psi^{-4}.\n",
    "        CCTK_REAL gupxz_iphjph = avg(INTERP_VARS[GUPXZI], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS0);\n",
    "        CCTK_REAL gupyz_iphjph = avg(INTERP_VARS[GUPYZI], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS0);\n",
    "        CCTK_REAL gupzz_iphjph = avg(INTERP_VARS[GUPZZI], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS0);\n",
    "        //CCTK_REAL lapse_iphjph = avg(INTERP_VARS[LAPM1I], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS0)+1.0;\n",
    "        //CCTK_REAL psi_iphjph   = avg(INTERP_VARS[PSII  ], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS0);\n",
    "\n",
    "        CCTK_REAL lapse_Psi2_iphjph = avg(INTERP_VARS[LAPSE_PSI2I], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS0);\n",
    "\n",
    "        CCTK_REAL A_x_iphjph   = avg(INTERP_VARS[A_XI], PLUS0,PLUS1, PLUS0,PLUS0,MINUS1,PLUS0); // @ (i,j+1/2,k+1/2)\n",
    "        CCTK_REAL A_y_iphjph   = avg(INTERP_VARS[A_YI], PLUS0,PLUS0, PLUS0,PLUS1,MINUS1,PLUS0); // @ (i+1/2,j,k+1/2)\n",
    "        CCTK_REAL A_z_iphjph   = avg(INTERP_VARS[A_ZI], PLUS0,PLUS0, PLUS0,PLUS0, PLUS0,PLUS0); // @ (i+1/2,j+1/2,k)\n",
    "\n",
    "        alpha_sqrtg_Az_interp[index] = lapse_Psi2_iphjph*\n",
    "          ( gupxz_iphjph*A_x_iphjph + gupyz_iphjph*A_y_iphjph + gupzz_iphjph*A_z_iphjph );\n",
    "\n",
    "\n",
    "        // Next set \\alpha \\Phi - \\beta^j A_j at (i+1/2,j+1/2,k+1/2):\n",
    "        //   We add a \"L\" suffix to shifti_iphjphkph to denote \"LOCAL\", as we set\n",
    "        //      shifti_iphjphkph[] gridfunction below.\n",
    "        CCTK_REAL shiftx_iphjphkphL = avg(INTERP_VARS[SHIFTXI], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS1);\n",
    "        CCTK_REAL shifty_iphjphkphL = avg(INTERP_VARS[SHIFTYI], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS1);\n",
    "        CCTK_REAL shiftz_iphjphkphL = avg(INTERP_VARS[SHIFTZI], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS1);\n",
    "        CCTK_REAL lapse_over_Psi6_iphjphkphL = avg(INTERP_VARS[LAPSE_OVER_PSI6I], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS1);\n",
    "        //CCTK_REAL psi_iphjphkph = avg(INTERP_VARS[PSII  ], PLUS0,PLUS1, PLUS0,PLUS1, PLUS0,PLUS1);\n",
    "        //CCTK_REAL psi2_iphjphkph= psi_iphjphkph*psi_iphjphkph;\n",
    "        //CCTK_REAL psi6_iphjphkph= psi2_iphjphkph*psi2_iphjphkph*psi2_iphjphkph;\n",
    "        CCTK_REAL A_x_iphjphkph = avg(INTERP_VARS[A_XI], PLUS0,PLUS1, PLUS0,PLUS0, PLUS0,PLUS0); // @ (i,j+1/2,k+1/2)\n",
    "        CCTK_REAL A_y_iphjphkph = avg(INTERP_VARS[A_YI], PLUS0,PLUS0, PLUS0,PLUS1, PLUS0,PLUS0); // @ (i+1/2,j,k+1/2)\n",
    "        CCTK_REAL A_z_iphjphkph = avg(INTERP_VARS[A_ZI], PLUS0,PLUS0, PLUS0,PLUS0, PLUS0,PLUS1); // @ (i+1/2,j+1/2,k)\n",
    "\n",
    "        alpha_Phi_minus_betaj_A_j_iphjphkph[index] = psi6phi[index]*lapse_over_Psi6_iphjphkphL\n",
    "          - (shiftx_iphjphkphL*A_x_iphjphkph + shifty_iphjphkphL*A_y_iphjphkph + shiftz_iphjphkphL*A_z_iphjphkph);\n",
    "\n",
    "        // Finally, save shifti_iphjphkph, for \\partial_j \\beta^j psi6phi\n",
    "        shiftx_iphjphkph[index]=shiftx_iphjphkphL;\n",
    "        shifty_iphjphkph[index]=shifty_iphjphkphL;\n",
    "        shiftz_iphjphkph[index]=shiftz_iphjphkphL;\n",
    "      }\n",
    "\n",
    "  // This loop requires two additional ghostzones in every direction. Hence the following loop definition:\n",
    "#pragma omp parallel for\n",
    "  for(int k=cctk_nghostzones[2];k<cctk_lsh[2]-cctk_nghostzones[2];k++) for(int j=cctk_nghostzones[1];j<cctk_lsh[1]-cctk_nghostzones[1];j++) for(int i=cctk_nghostzones[0];i<cctk_lsh[0]-cctk_nghostzones[0];i++) {\n",
    "        int index = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "\n",
    "        // \\partial_t A_i = [reconstructed stuff] + [gauge stuff],\n",
    "        //    where [gauge stuff] = -\\partial_i (\\alpha \\Phi - \\beta^j A_j)\n",
    "        CCTK_REAL alpha_Phi_minus_betaj_A_j_iphjphkphL = alpha_Phi_minus_betaj_A_j_iphjphkph[index];\n",
    "        // - partial_i -> - (A_{i} - A_{i-1})/dX = (A_{i-1} - A_{i})/dX, for Ax\n",
    "        Ax_rhs[index] += dXm1*(alpha_Phi_minus_betaj_A_j_iphjphkph[CCTK_GFINDEX3D(cctkGH,i-1,j,k)] - alpha_Phi_minus_betaj_A_j_iphjphkphL);\n",
    "        Ay_rhs[index] += dYm1*(alpha_Phi_minus_betaj_A_j_iphjphkph[CCTK_GFINDEX3D(cctkGH,i,j-1,k)] - alpha_Phi_minus_betaj_A_j_iphjphkphL);\n",
    "        Az_rhs[index] += dZm1*(alpha_Phi_minus_betaj_A_j_iphjphkph[CCTK_GFINDEX3D(cctkGH,i,j,k-1)] - alpha_Phi_minus_betaj_A_j_iphjphkphL);\n",
    "\n",
    "        // \\partial_t psi6phi = [shift advection term] + \\partial_j (\\alpha \\sqrt{\\gamma} A^j)\n",
    "        // Here we compute [shift advection term] = \\partial_j (\\beta^j psi6phi)\n",
    "        // Cache misses are likely more expensive than branch mispredictions here,\n",
    "        //       which is why we use if() statements and array lookups inside the if()'s.\n",
    "        CCTK_REAL psi6phi_rhsL=0.0;\n",
    "        CCTK_REAL psi6phiL=psi6phi[index];\n",
    "        CCTK_REAL shiftx_iphjphkphL=shiftx_iphjphkph[index];\n",
    "        CCTK_REAL shifty_iphjphkphL=shifty_iphjphkph[index];\n",
    "        CCTK_REAL shiftz_iphjphkphL=shiftz_iphjphkph[index];\n",
    "\n",
    "        // \\partial_x (\\beta^x psi6phi) :\n",
    "        if(shiftx_iphjphkphL < 0.0) {\n",
    "          psi6phi_rhsL+=0.5*dXm1*(+    shiftx_iphjphkph[CCTK_GFINDEX3D(cctkGH,i-2,j,k)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i-2,j,k)]\n",
    "                                  -4.0*shiftx_iphjphkph[CCTK_GFINDEX3D(cctkGH,i-1,j,k)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i-1,j,k)]\n",
    "                                  +3.0*shiftx_iphjphkphL*                               psi6phiL);\n",
    "        } else {\n",
    "          psi6phi_rhsL+=0.5*dXm1*(-    shiftx_iphjphkph[CCTK_GFINDEX3D(cctkGH,i+2,j,k)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i+2,j,k)]\n",
    "                                  +4.0*shiftx_iphjphkph[CCTK_GFINDEX3D(cctkGH,i+1,j,k)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i+1,j,k)]\n",
    "                                  -3.0*shiftx_iphjphkphL*                               psi6phiL);\n",
    "        }\n",
    "\n",
    "        // \\partial_y (\\beta^y psi6phi) :\n",
    "        if(shifty_iphjphkphL < 0.0) {\n",
    "          psi6phi_rhsL+=0.5*dYm1*(+    shifty_iphjphkph[CCTK_GFINDEX3D(cctkGH,i,j-2,k)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i,j-2,k)]\n",
    "                                  -4.0*shifty_iphjphkph[CCTK_GFINDEX3D(cctkGH,i,j-1,k)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i,j-1,k)]\n",
    "                                  +3.0*shifty_iphjphkphL*                               psi6phiL);\n",
    "        } else {\n",
    "          psi6phi_rhsL+=0.5*dYm1*(-    shifty_iphjphkph[CCTK_GFINDEX3D(cctkGH,i,j+2,k)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i,j+2,k)]\n",
    "                                  +4.0*shifty_iphjphkph[CCTK_GFINDEX3D(cctkGH,i,j+1,k)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i,j+1,k)]\n",
    "                                  -3.0*shifty_iphjphkphL*                               psi6phiL);\n",
    "        }\n",
    "\n",
    "        // \\partial_z (\\beta^z psi6phi) :\n",
    "        if(shiftz_iphjphkphL < 0.0) {\n",
    "          psi6phi_rhsL+=0.5*dZm1*(+    shiftz_iphjphkph[CCTK_GFINDEX3D(cctkGH,i,j,k-2)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i,j,k-2)]\n",
    "                                  -4.0*shiftz_iphjphkph[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i,j,k-1)]\n",
    "                                  +3.0*shiftz_iphjphkphL*                               psi6phiL);\n",
    "        } else {\n",
    "          psi6phi_rhsL+=0.5*dZm1*(-    shiftz_iphjphkph[CCTK_GFINDEX3D(cctkGH,i,j,k+2)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i,j,k+2)]\n",
    "                                  +4.0*shiftz_iphjphkph[CCTK_GFINDEX3D(cctkGH,i,j,k+1)]*psi6phi[CCTK_GFINDEX3D(cctkGH,i,j,k+1)]\n",
    "                                  -3.0*shiftz_iphjphkphL*                               psi6phiL);\n",
    "        }\n",
    "\n",
    "        // Next we add \\partial_j (\\alpha \\sqrt{\\gamma} A^j) to \\partial_t psi6phi:\n",
    "        psi6phi_rhsL+=dXm1*(alpha_sqrtg_Ax_interp[index] - alpha_sqrtg_Ax_interp[CCTK_GFINDEX3D(cctkGH,i+1,j,k)])\n",
    "          +           dYm1*(alpha_sqrtg_Ay_interp[index] - alpha_sqrtg_Ay_interp[CCTK_GFINDEX3D(cctkGH,i,j+1,k)])\n",
    "          +           dZm1*(alpha_sqrtg_Az_interp[index] - alpha_sqrtg_Az_interp[CCTK_GFINDEX3D(cctkGH,i,j,k+1)]);\n",
    "\n",
    "        // *GENERALIZED* LORENZ GAUGE:\n",
    "        // Finally, add damping factor to \\partial_t psi6phi\n",
    "        //subtract lambda * alpha psi^6 Phi\n",
    "        psi6phi_rhsL+=-damp_lorenz*alpha_iphjphkph[index]*psi6phiL;\n",
    "\n",
    "        psi6phi_rhs[index] = psi6phi_rhsL;\n",
    "      }\n",
    "}\n",
    "\n",
    "static inline CCTK_REAL avg(CCTK_REAL f[PLUS2+1][PLUS2+1][PLUS2+1],int imin,int imax, int jmin,int jmax, int kmin,int kmax) {\n",
    "  CCTK_REAL retval=0.0,num_in_sum=0.0;\n",
    "  for(int kk=kmin;kk<=kmax;kk++) for(int jj=jmin;jj<=jmax;jj++) for(int ii=imin;ii<=imax;ii++) {\n",
    "        retval+=f[kk][jj][ii]; num_in_sum++;\n",
    "      }\n",
    "  return retval/num_in_sum;\n",
    "}\n",
    "\n",
    "\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": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Validation test for Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C: FAILED!\n",
      "Diff:\n",
      "8c8\n",
      "<   DECLARE_CCTK_PARAMETERS;  \n",
      "---\n",
      ">   DECLARE_CCTK_PARAMETERS;\n",
      "11c11\n",
      "<    * Note that the RHS consists of a shift advection term on psi6phi and \n",
      "---\n",
      ">    * Note that the RHS consists of a shift advection term on psi6phi and\n",
      "13c13\n",
      "<    * psi6phi is defined at (i+1/2,j+1/2,k+1/2), but instead of reconstructing \n",
      "---\n",
      ">    * psi6phi is defined at (i+1/2,j+1/2,k+1/2), but instead of reconstructing\n",
      "29c29\n",
      "<     \n",
      "---\n",
      "> \n",
      "38c38\n",
      "<         //    (i,j+1/2,k+1/2)and (i+1,j+1/2,k+1/2), then taking \\partial_x (RHS1x) = \n",
      "---\n",
      ">         //    (i,j+1/2,k+1/2)and (i+1,j+1/2,k+1/2), then taking \\partial_x (RHS1x) =\n",
      "43c43\n",
      "<         num_vars_to_interp = 11; \n",
      "---\n",
      ">         num_vars_to_interp = 11;\n",
      "74c74\n",
      "<         for(int kk=PLUS0;kk<=PLUS1;kk++) for(int jj=PLUS0;jj<=PLUS1;jj++) for(int ii=PLUS0;ii<=PLUS1;ii++) { \n",
      "---\n",
      ">         for(int kk=PLUS0;kk<=PLUS1;kk++) for(int jj=PLUS0;jj<=PLUS1;jj++) for(int ii=PLUS0;ii<=PLUS1;ii++) {\n",
      "105c105\n",
      "<     \n",
      "---\n",
      "> \n",
      "123c123\n",
      "<     \n",
      "---\n",
      "> \n",
      "129c129\n",
      "<         //   We add a \"L\" suffix to shifti_iphjphkph to denote \"LOCAL\", as we set \n",
      "---\n",
      ">         //   We add a \"L\" suffix to shifti_iphjphkph to denote \"LOCAL\", as we set\n",
      "195c195\n",
      "<  \n",
      "---\n",
      "> \n",
      "212c212\n",
      "<         // *GENERALIZED* LORENZ GAUGE: \n",
      "---\n",
      ">         // *GENERALIZED* LORENZ GAUGE:\n",
      "227a228\n",
      "> \n"
     ]
    }
   ],
   "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/Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C\"\n",
    "original_IGM_file_name = \"Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs-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__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs__C  = !diff $original_IGM_file_path $outfile_path__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs__C\n",
    "\n",
    "if Validation__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs__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 Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C: PASSED!\")\n",
    "else:\n",
    "    # If the validation fails, we keep the original IGM source code file\n",
    "    print(\"Validation test for Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.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__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs__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__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.pdf](Tutorial-IllinoisGRMHD__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.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": 4,
   "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__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.ipynb\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.tex\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.tex\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.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.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}