{
 "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: add_fluxes_and_source_terms_to_hydro_rhss.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 demonstrates the addition of the flux and source terms to the right-hand side of the hydrodynamic variables. \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](#add_fluxes_and_source_terms_to_hydro_rhss__c): **`add_fluxes_and_source_terms_to_hydro_rhss.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__add_fluxes_and_source_terms_to_hydro_rhss__C = os.path.join(IGM_src_dir_path,\"add_fluxes_and_source_terms_to_hydro_rhss.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='add_fluxes_and_source_terms_to_hydro_rhss__c'></a>\n",
    "\n",
    "# Step 2: `add_fluxes_and_source_terms_to_hydro_rhss.C` \\[Back to [top](#toc)\\]\n",
    "$$\\label{add_fluxes_and_source_terms_to_hydro_rhss__c}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing ../src/add_fluxes_and_source_terms_to_hydro_rhss.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outfile_path__add_fluxes_and_source_terms_to_hydro_rhss__C\n",
    "// Side note: the following values could be used for cell averaged gfs:\n",
    "//     am2=-1.0/12.0, am1=7.0/12.0, a0=7.0/12.0, a1=-1.0/12.0\n",
    "// However, since the metric gfs store the grid point values instead of the cell average,\n",
    "//     the following coefficients should be used:\n",
    "//     am2 = -1/16, am1 = 9/16, a0 = 9/16, a1 = -1/16\n",
    "// This will yield the third-order-accurate face values at m-1/2,\n",
    "//      using values specified at {m-2,m-1,m,m+1}\n",
    "#define AM2 -0.0625\n",
    "#define AM1  0.5625\n",
    "#define A0   0.5625\n",
    "#define A1  -0.0625\n",
    "#define COMPUTE_FCVAL(METRICm2,METRICm1,METRIC,METRICp1) (AM2*(METRICm2) + AM1*(METRICm1) + A0*(METRIC) + A1*(METRICp1))\n",
    "\n",
    "static inline void mhdflux(int i,int j,int k,const int flux_dirn,CCTK_REAL *Ul,CCTK_REAL *Ur,  CCTK_REAL *FACEVAL,CCTK_REAL *FACEVAL_LAPSE_PSI4,eos_struct &eos, CCTK_REAL Gamma_th,\n",
    "                           CCTK_REAL &cmax,CCTK_REAL &cmin,\n",
    "                           CCTK_REAL &rho_star_flux,CCTK_REAL &tau_flux,CCTK_REAL &st_x_flux,CCTK_REAL &st_y_flux,CCTK_REAL &st_z_flux);\n",
    "\n",
    "#define COMPUTE_FOURMETRIC(g4tt,g4tx,g4ty,g4tz,g4xx,g4xy,g4xz,g4yy,g4yz,g4zz,METRIC,METRIC_AUX)  ( { \\\n",
    "      /* g_{0i} = beta_i */                                             \\\n",
    "      g4tx = METRIC_AUX[PSI4]*(METRIC[GXX]*METRIC[SHIFTX] + METRIC[GXY]*METRIC[SHIFTY] + METRIC[GXZ]*METRIC[SHIFTZ]); \\\n",
    "      g4ty = METRIC_AUX[PSI4]*(METRIC[GXY]*METRIC[SHIFTX] + METRIC[GYY]*METRIC[SHIFTY] + METRIC[GYZ]*METRIC[SHIFTZ]); \\\n",
    "      g4tz = METRIC_AUX[PSI4]*(METRIC[GXZ]*METRIC[SHIFTX] + METRIC[GYZ]*METRIC[SHIFTY] + METRIC[GZZ]*METRIC[SHIFTZ]); \\\n",
    "      /* g_{00} = -alpha^2 + beta^i beta^j gamma_{ij} = -alpha^2 + beta^i beta_i = -alpha^2 + beta^i g_{0i} */ \\\n",
    "      g4tt = -SQR(METRIC_AUX[LAPSE]) + g4tx*METRIC[SHIFTX] + g4ty*METRIC[SHIFTY] + g4tz*METRIC[SHIFTZ]; \\\n",
    "      g4xx = METRIC_AUX[PSI4]*METRIC[GXX];                              \\\n",
    "      g4xy = METRIC_AUX[PSI4]*METRIC[GXY];                              \\\n",
    "      g4xz = METRIC_AUX[PSI4]*METRIC[GXZ];                              \\\n",
    "      g4yy = METRIC_AUX[PSI4]*METRIC[GYY];                              \\\n",
    "      g4yz = METRIC_AUX[PSI4]*METRIC[GYZ];                              \\\n",
    "      g4zz = METRIC_AUX[PSI4]*METRIC[GZZ];                              \\\n",
    "    } )\n",
    "\n",
    "\n",
    "static void add_fluxes_and_source_terms_to_hydro_rhss(const int flux_dirn,const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,CCTK_REAL *dX,\n",
    "                                                      CCTK_REAL **metric,gf_and_gz_struct *in_prims,CCTK_REAL **TUPmunu,\n",
    "                                                      int numvars_reconstructed,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,eos_struct &eos,\n",
    "                                                      CCTK_REAL *cmax,CCTK_REAL *cmin,\n",
    "                                                      CCTK_REAL *rho_star_flux,CCTK_REAL *tau_flux,CCTK_REAL *st_x_flux,CCTK_REAL *st_y_flux,CCTK_REAL *st_z_flux,\n",
    "                                                      CCTK_REAL *rho_star_rhs,CCTK_REAL *tau_rhs,CCTK_REAL *st_x_rhs,CCTK_REAL *st_y_rhs,CCTK_REAL *st_z_rhs) {\n",
    "\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  CCTK_REAL dxi[4] = { 1e100,1.0/dX[0],1.0/dX[1],1.0/dX[2] };\n",
    "\n",
    "  // Notice in the loop below that we go from 3 to cctk_lsh-2 for i, j, AND k, even though\n",
    "  //   we are only computing the flux in one direction at a time. This is because in the end,\n",
    "  //   we only need the rhs's from 3 to cctk_lsh-3 for i, j, and k.\n",
    "#pragma omp parallel for\n",
    "  for(int k=cctk_nghostzones[2];k<cctk_lsh[2]-(cctk_nghostzones[2]-1);k++) for(int j=cctk_nghostzones[1];j<cctk_lsh[1]-(cctk_nghostzones[1]-1);j++) for(int i=cctk_nghostzones[0];i<cctk_lsh[0]-(cctk_nghostzones[0]-1);i++) {\n",
    "\tint index = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "\n",
    "\t// Set metric and associated variables\n",
    "\tCCTK_REAL METRIC[NUMVARS_FOR_METRIC_FACEVALS]; for(int ii=0;ii<NUMVARS_FOR_METRIC_FACEVALS;ii++) METRIC[ii] = metric[ii][index];\n",
    "\tCCTK_REAL METRIC_LAP_PSI4[NUMVARS_METRIC_AUX]; SET_LAPSE_PSI4(METRIC_LAP_PSI4,METRIC);\n",
    "\n",
    "\tCCTK_REAL Ur[MAXNUMVARS]; for(int ii=0;ii<numvars_reconstructed;ii++) Ur[ii] = out_prims_r[ii].gf[index];\n",
    "\tCCTK_REAL Ul[MAXNUMVARS]; for(int ii=0;ii<numvars_reconstructed;ii++) Ul[ii] = out_prims_l[ii].gf[index];\n",
    "\n",
    "\t// Read the T^{\\mu \\nu} gridfunction from memory, since computing T^{\\mu \\nu} is expensive\n",
    "\tCCTK_REAL TUP[4][4]; int counter=0;\n",
    "\tfor(int ii=0;ii<4;ii++) for(int jj=ii;jj<4;jj++) { TUP[ii][jj] = TUPmunu[counter][index]; counter++; }\n",
    "\n",
    "\t// Next set metric on the faces, applying a 3rd-order lopsided stencil.\n",
    "\tint indexm2 = CCTK_GFINDEX3D(cctkGH,i-2*kronecker_delta[flux_dirn][0],j-2*kronecker_delta[flux_dirn][1],k-2*kronecker_delta[flux_dirn][2]);\n",
    "\tint indexm1 = CCTK_GFINDEX3D(cctkGH,i-  kronecker_delta[flux_dirn][0],j-  kronecker_delta[flux_dirn][1],k-  kronecker_delta[flux_dirn][2]);\n",
    "\tint indexp1 = CCTK_GFINDEX3D(cctkGH,i+  kronecker_delta[flux_dirn][0],j+  kronecker_delta[flux_dirn][1],k+  kronecker_delta[flux_dirn][2]);\n",
    "\tint indexp2 = CCTK_GFINDEX3D(cctkGH,i+2*kronecker_delta[flux_dirn][0],j+2*kronecker_delta[flux_dirn][1],k+2*kronecker_delta[flux_dirn][2]);\n",
    "\t// The \"vector\" METRIC stores needed metric-related quantities.\n",
    "\tCCTK_REAL METRICm2[NUMVARS_FOR_METRIC_FACEVALS]; for(int ii=0;ii<NUMVARS_FOR_METRIC_FACEVALS;ii++) METRICm2[ii] = metric[ii][indexm2];\n",
    "\tCCTK_REAL METRICm1[NUMVARS_FOR_METRIC_FACEVALS]; for(int ii=0;ii<NUMVARS_FOR_METRIC_FACEVALS;ii++) METRICm1[ii] = metric[ii][indexm1];\n",
    "\tCCTK_REAL METRICp1[NUMVARS_FOR_METRIC_FACEVALS]; for(int ii=0;ii<NUMVARS_FOR_METRIC_FACEVALS;ii++) METRICp1[ii] = metric[ii][indexp1];\n",
    "\tCCTK_REAL METRICp2[NUMVARS_FOR_METRIC_FACEVALS]; for(int ii=0;ii<NUMVARS_FOR_METRIC_FACEVALS;ii++) METRICp2[ii] = metric[ii][indexp2];\n",
    "\n",
    "\t// Next compute the metric values at the {i,j,k} +/- 1/2 faces (i.e., the \"face values\" of the metric)\n",
    "\tCCTK_REAL FACEVAL[NUMVARS_FOR_METRIC_FACEVALS],FACEVALp1[NUMVARS_FOR_METRIC_FACEVALS];\n",
    "\tfor(int w=0;w<NUMVARS_FOR_METRIC_FACEVALS;w++) FACEVAL[w]   = COMPUTE_FCVAL(METRICm2[w],METRICm1[w],METRIC[w],METRICp1[w]);\n",
    "\tfor(int w=0;w<NUMVARS_FOR_METRIC_FACEVALS;w++) FACEVALp1[w] = COMPUTE_FCVAL(METRICm1[w],METRIC[w],METRICp1[w],METRICp2[w]);\n",
    "\t// Then compute the lapse and Psi4 = exp(4*phi)\n",
    "\tCCTK_REAL FACEVAL_LAPSE_PSI4[NUMVARS_METRIC_AUX],FACEVAL_LAPSE_PSI4p1[NUMVARS_METRIC_AUX];\n",
    "\tSET_LAPSE_PSI4(FACEVAL_LAPSE_PSI4,  FACEVAL);\n",
    "\tSET_LAPSE_PSI4(FACEVAL_LAPSE_PSI4p1,FACEVALp1);\n",
    "\n",
    "\t//-----------------------------------------------------------------------------\n",
    "\t// Next compute fluxes for \\tilde{S}_i, tau, and rho_*\n",
    "\tmhdflux(i,j,k,flux_dirn,Ul  ,Ur  ,FACEVAL  ,FACEVAL_LAPSE_PSI4  ,eos, Gamma_th, cmax[index],cmin[index],\n",
    "\t\trho_star_flux[index],tau_flux[index],st_x_flux[index],st_y_flux[index],st_z_flux[index]);\n",
    "\n",
    "\n",
    "\t//-----------------------------------------------------------------------------\n",
    "\t// If we are not in the ghostzones, then add third-order accurate curvature terms to \\tilde{S}_i RHS's\n",
    "\t//    Without this if() statement, _rhs variables are in general set to nonzero values in ghostzones, which messes up frozen BC's.\n",
    "\t//    Also, this if() statement should speed up the computation slightly.\n",
    "\tif(k<cctk_lsh[2]-cctk_nghostzones[2] && j<cctk_lsh[1]-cctk_nghostzones[1] && i<cctk_lsh[0]-cctk_nghostzones[0]) {\n",
    "\n",
    "\t  CCTK_REAL Psi6 = METRIC_LAP_PSI4[PSI2]*METRIC_LAP_PSI4[PSI4];\n",
    "\t  CCTK_REAL half_alpha_sqrtgamma = 0.5*METRIC_LAP_PSI4[LAPSE]*Psi6;\n",
    "\n",
    "\t  // First compute four metric.\n",
    "\t  CCTK_REAL g4tt_f,g4tx_f,g4ty_f,g4tz_f,g4xx_f,g4xy_f,g4xz_f,g4yy_f,g4yz_f,g4zz_f;\n",
    "\t  COMPUTE_FOURMETRIC(g4tt_f,g4tx_f,g4ty_f,g4tz_f,g4xx_f,g4xy_f,g4xz_f,g4yy_f,g4yz_f,g4zz_f,FACEVAL,FACEVAL_LAPSE_PSI4);\n",
    "\n",
    "\t  CCTK_REAL g4tt_fp1,g4tx_fp1,g4ty_fp1,g4tz_fp1,g4xx_fp1,g4xy_fp1,g4xz_fp1,g4yy_fp1,g4yz_fp1,g4zz_fp1;\n",
    "\t  COMPUTE_FOURMETRIC(g4tt_fp1,g4tx_fp1,g4ty_fp1,g4tz_fp1,g4xx_fp1,g4xy_fp1,g4xz_fp1,g4yy_fp1,g4yz_fp1,g4zz_fp1,FACEVALp1,FACEVAL_LAPSE_PSI4p1);\n",
    "\n",
    "\t  // Compute \\partial_i g_{\\mu \\nu} at m+1/2\n",
    "\t  CCTK_REAL partial_i_gmunu[4][4];\n",
    "\t  partial_i_gmunu[0][0] = (g4tt_fp1 - g4tt_f)*dxi[flux_dirn];\n",
    "\t  partial_i_gmunu[0][1] = (g4tx_fp1 - g4tx_f)*dxi[flux_dirn];\n",
    "\t  partial_i_gmunu[0][2] = (g4ty_fp1 - g4ty_f)*dxi[flux_dirn];\n",
    "\t  partial_i_gmunu[0][3] = (g4tz_fp1 - g4tz_f)*dxi[flux_dirn];\n",
    "\t  partial_i_gmunu[1][1] = (g4xx_fp1 - g4xx_f)*dxi[flux_dirn];\n",
    "\t  partial_i_gmunu[1][2] = (g4xy_fp1 - g4xy_f)*dxi[flux_dirn];\n",
    "\t  partial_i_gmunu[1][3] = (g4xz_fp1 - g4xz_f)*dxi[flux_dirn];\n",
    "\t  partial_i_gmunu[2][2] = (g4yy_fp1 - g4yy_f)*dxi[flux_dirn];\n",
    "\t  partial_i_gmunu[2][3] = (g4yz_fp1 - g4yz_f)*dxi[flux_dirn];\n",
    "\t  partial_i_gmunu[3][3] = (g4zz_fp1 - g4zz_f)*dxi[flux_dirn];\n",
    "\n",
    "\t  // Needed for tau_rhs computation:\n",
    "\t  CCTK_REAL lapse_deriv[4] = { 0,0,0,0 };\n",
    "\t  lapse_deriv[flux_dirn] = (FACEVALp1[LAPM1] - FACEVAL[LAPM1])*dxi[flux_dirn];\n",
    "\n",
    "\t  // Needed for st_i_rhs computation:\n",
    "\t  CCTK_REAL st_i_curvature_terms[4] = { 0,0,0,0 };\n",
    "\t  // add \\frac{1}{2} \\alpha \\sqrt{\\gamma} T^{\\mu \\nu} \\partial_i g_{\\mu \\nu} . Note that i is given by the flux direction.\n",
    "\t  //   (Source term of Eq 43 in http://arxiv.org/pdf/astro-ph/0503420.pdf)\n",
    "\t  st_i_curvature_terms[flux_dirn] = half_alpha_sqrtgamma * ( TUP[0][0]*partial_i_gmunu[0][0] +\n",
    "\t\t\t\t\t\t\t\t     TUP[1][1]*partial_i_gmunu[1][1] +\n",
    "\t\t\t\t\t\t\t\t     TUP[2][2]*partial_i_gmunu[2][2] +\n",
    "\t\t\t\t\t\t\t\t     TUP[3][3]*partial_i_gmunu[3][3] +\n",
    "\t\t\t\t\t\t\t\t     2.0*(TUP[0][1]*partial_i_gmunu[0][1] +\n",
    "\t\t\t\t\t\t\t\t\t  TUP[0][2]*partial_i_gmunu[0][2] +\n",
    "\t\t\t\t\t\t\t\t\t  TUP[0][3]*partial_i_gmunu[0][3] +\n",
    "\t\t\t\t\t\t\t\t\t  TUP[1][2]*partial_i_gmunu[1][2] +\n",
    "\t\t\t\t\t\t\t\t\t  TUP[1][3]*partial_i_gmunu[1][3] +\n",
    "\t\t\t\t\t\t\t\t\t  TUP[2][3]*partial_i_gmunu[2][3]) );\n",
    "\n",
    "\t  // add - ( T^{00} \\beta^i + T^{0i} ) \\partial_i \\alpha.\n",
    "\t  //   (Last part of Eq. 39 source term in http://arxiv.org/pdf/astro-ph/0503420.pdf)\n",
    "\t  CCTK_REAL alpha_sqrtgamma = 2.0*half_alpha_sqrtgamma;\n",
    "\t  tau_rhs[index]  += alpha_sqrtgamma*(-(TUP[0][0]*METRIC[SHIFTX+(flux_dirn-1)] + TUP[0][flux_dirn])*lapse_deriv[flux_dirn]);\n",
    "\n",
    "\t  // Eq 43 in http://arxiv.org/pdf/astro-ph/0503420.pdf:\n",
    "\t  // \\partial_t \\tilde{S}_i = - \\partial_i (\\alpha \\sqrt{\\gamma} T^j_i) + \\frac{1}{2}\\alpha \\sqrt{\\gamma} T^{\\mu \\nu}g_{\\mu \\nu,i}\n",
    "\t  // Notice that st_i_curvature_terms[N]=0 for N!=flux_dirn.\n",
    "\t  st_x_rhs[index] += st_i_curvature_terms[1];\n",
    "\t  st_y_rhs[index] += st_i_curvature_terms[2];\n",
    "\t  st_z_rhs[index] += st_i_curvature_terms[3];\n",
    "\t}\n",
    "\n",
    "      }\n",
    "\n",
    "  // Notice in the loop below that we go from 3 to cctk_lsh-3 for i, j, AND k, even though\n",
    "  //   we are only computing the flux in one direction. This is because in the end,\n",
    "  //   we only need the rhs's from 3 to cctk_lsh-3 for i, j, and k.\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",
    "\tint index = CCTK_GFINDEX3D(cctkGH,i,j,k);\n",
    "\tint indexp1 = CCTK_GFINDEX3D(cctkGH,i+kronecker_delta[flux_dirn][0],j+kronecker_delta[flux_dirn][1],k+kronecker_delta[flux_dirn][2]);\n",
    "\n",
    "\trho_star_rhs[index] += (rho_star_flux[index] - rho_star_flux[indexp1]) * dxi[flux_dirn];\n",
    "\ttau_rhs[index]      += (tau_flux[index]      - tau_flux[indexp1]     ) * dxi[flux_dirn];\n",
    "\tst_x_rhs[index]     += (st_x_flux[index]     - st_x_flux[indexp1]    ) * dxi[flux_dirn];\n",
    "\tst_y_rhs[index]     += (st_y_flux[index]     - st_y_flux[indexp1]    ) * dxi[flux_dirn];\n",
    "\tst_z_rhs[index]     += (st_z_flux[index]     - st_z_flux[indexp1]    ) * dxi[flux_dirn];\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 add_fluxes_and_source_terms_to_hydro_rhss.C: FAILED!\n",
      "Diff:\n",
      "1c1\n",
      "< // Side note: the following values could be used for cell averaged gfs: \n",
      "---\n",
      "> // Side note: the following values could be used for cell averaged gfs:\n",
      "3,4c3,4\n",
      "< // However, since the metric gfs store the grid point values instead of the cell average, \n",
      "< //     the following coefficients should be used: \n",
      "---\n",
      "> // However, since the metric gfs store the grid point values instead of the cell average,\n",
      "> //     the following coefficients should be used:\n",
      "6c6\n",
      "< // This will yield the third-order-accurate face values at m-1/2, \n",
      "---\n",
      "> // This will yield the third-order-accurate face values at m-1/2,\n",
      "14c14\n",
      "< static inline void mhdflux(int i,int j,int k,const int flux_dirn,CCTK_REAL *Ul,CCTK_REAL *Ur,  CCTK_REAL *FACEVAL,CCTK_REAL *FACEVAL_LAPSE_PSI4,eos_struct &eos,\n",
      "---\n",
      "> static inline void mhdflux(int i,int j,int k,const int flux_dirn,CCTK_REAL *Ul,CCTK_REAL *Ur,  CCTK_REAL *FACEVAL,CCTK_REAL *FACEVAL_LAPSE_PSI4,eos_struct &eos, CCTK_REAL Gamma_th,\n",
      "85c85\n",
      "< \tmhdflux(i,j,k,flux_dirn,Ul  ,Ur  ,FACEVAL  ,FACEVAL_LAPSE_PSI4  ,eos, cmax[index],cmin[index],\n",
      "---\n",
      "> \tmhdflux(i,j,k,flux_dirn,Ul  ,Ur  ,FACEVAL  ,FACEVAL_LAPSE_PSI4  ,eos, Gamma_th, cmax[index],cmin[index],\n",
      "88c88\n",
      "<     \n",
      "---\n",
      "> \n",
      "96c96\n",
      "< \t  CCTK_REAL half_alpha_sqrtgamma = 0.5*METRIC_LAP_PSI4[LAPSE]*Psi6;    \n",
      "---\n",
      "> \t  CCTK_REAL half_alpha_sqrtgamma = 0.5*METRIC_LAP_PSI4[LAPSE]*Psi6;\n",
      "101c101\n",
      "<       \n",
      "---\n",
      "> \n",
      "126,127c126,127\n",
      "< \t  st_i_curvature_terms[flux_dirn] = half_alpha_sqrtgamma * ( TUP[0][0]*partial_i_gmunu[0][0] + \n",
      "< \t\t\t\t\t\t\t\t     TUP[1][1]*partial_i_gmunu[1][1] + \n",
      "---\n",
      "> \t  st_i_curvature_terms[flux_dirn] = half_alpha_sqrtgamma * ( TUP[0][0]*partial_i_gmunu[0][0] +\n",
      "> \t\t\t\t\t\t\t\t     TUP[1][1]*partial_i_gmunu[1][1] +\n",
      "130,131c130,131\n",
      "< \t\t\t\t\t\t\t\t     2.0*(TUP[0][1]*partial_i_gmunu[0][1] + \n",
      "< \t\t\t\t\t\t\t\t\t  TUP[0][2]*partial_i_gmunu[0][2] + \n",
      "---\n",
      "> \t\t\t\t\t\t\t\t     2.0*(TUP[0][1]*partial_i_gmunu[0][1] +\n",
      "> \t\t\t\t\t\t\t\t\t  TUP[0][2]*partial_i_gmunu[0][2] +\n",
      "133,134c133,134\n",
      "< \t\t\t\t\t\t\t\t\t  TUP[1][2]*partial_i_gmunu[1][2] + \n",
      "< \t\t\t\t\t\t\t\t\t  TUP[1][3]*partial_i_gmunu[1][3] + \n",
      "---\n",
      "> \t\t\t\t\t\t\t\t\t  TUP[1][2]*partial_i_gmunu[1][2] +\n",
      "> \t\t\t\t\t\t\t\t\t  TUP[1][3]*partial_i_gmunu[1][3] +\n",
      "140c140\n",
      "< \t  tau_rhs[index]  += alpha_sqrtgamma*(-(TUP[0][0]*METRIC[SHIFTX+(flux_dirn-1)] + TUP[0][flux_dirn])*lapse_deriv[flux_dirn]); \n",
      "---\n",
      "> \t  tau_rhs[index]  += alpha_sqrtgamma*(-(TUP[0][0]*METRIC[SHIFTX+(flux_dirn-1)] + TUP[0][flux_dirn])*lapse_deriv[flux_dirn]);\n",
      "166a167\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/add_fluxes_and_source_terms_to_hydro_rhss.C\"\n",
    "original_IGM_file_name = \"add_fluxes_and_source_terms_to_hydro_rhss-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__add_fluxes_and_source_terms_to_hydro_rhss__C  = !diff $original_IGM_file_path $outfile_path__add_fluxes_and_source_terms_to_hydro_rhss__C\n",
    "\n",
    "if Validation__add_fluxes_and_source_terms_to_hydro_rhss__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 add_fluxes_and_source_terms_to_hydro_rhss.C: PASSED!\")\n",
    "else:\n",
    "    # If the validation fails, we keep the original IGM source code file\n",
    "    print(\"Validation test for add_fluxes_and_source_terms_to_hydro_rhss.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__add_fluxes_and_source_terms_to_hydro_rhss__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__add_fluxes_and_source_terms_to_hydro_rhss.pdf](Tutorial-IllinoisGRMHD__add_fluxes_and_source_terms_to_hydro_rhss.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__add_fluxes_and_source_terms_to_hydro_rhss.ipynb\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__add_fluxes_and_source_terms_to_hydro_rhss.tex\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__add_fluxes_and_source_terms_to_hydro_rhss.tex\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__add_fluxes_and_source_terms_to_hydro_rhss.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
}