#!/usr/bin/env python
# coding: utf-8

# <script async src="https://www.googletagmanager.com/gtag/js?id=UA-59152712-8"></script>
# <script>
#   window.dataLayer = window.dataLayer || [];
#   function gtag(){dataLayer.push(arguments);}
#   gtag('js', new Date());
# 
#   gtag('config', 'UA-59152712-8');
# </script>
# 
# # Tutorial-IllinoisGRMHD: `driver_evaluate_MHD_rhs.C`
# 
# ## Authors: Leo Werneck & Zach Etienne
# 
# <font color='red'>**This module is currently under development**</font>
# 
# ## This tutorial notebook presents the driver functions that compute the right-hand side (RHS) of the MHD equations within `IllinoisGRMHD`.
# 
# ### Required and recommended citations:
# 
# * **(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)).
# * **(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)).
# * **(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)).

# <a id='toc'></a>
# 
# # Table of Contents
# $$\label{toc}$$
# 
# This module is organized as follows
# 
# 0. [Step 0](#src_dir): **Source directory creation**
# 1. [Step 1](#introduction): **Introduction**
# 1. [Step 2](#header_files): **Load up necessary ETK and `IllinoisGRMHD` header files**
# 1. [Step 3](#driver_mhd_rhs_function): **The ` IllinoisGRMHD_driver_evaluate_MHD_rhs()` function**
#     1. [Step 3.1](#eos_parameters): *Equation of state (EOS) parameters*
#     1. [Step 3.2](#set_pointers_grmhd_gfs): *Set pointers to GRMHD gridfunctions*
#     1. [Step 3.3](#admbase_to_bssnbase): *Convert ADM variables to BSSN variables*
#     1. [Step 3.4](#pointers_metric_tmunu_gfs): *Setting up pointers to the metric and stress-energy tensor gridfunctions*
#     1. [Step 3.5](#initialize_rhss): *Initialization of the RHS variables*
#     1. [Step 3.6](#tau_rhs_ext_curv_and_tupmunu): *Compute extrinsic curvature terms from the RHS of $\partial_{t}\tilde\tau$ and $T^{\mu\nu}$*
#     1. [Step 3.7](#computing_ftilde): *Computing ${\rm ftilde}$*
#     1. [Step 3.8](#rhs_mhd_and_a_i): *The RHSs of $\rho_{\star}$, $\tilde\tau$, $\tilde{S}_{i}$, and $A_{i}$*
#         1. [Step 3.8.1](#reconstructing_vx_vy_by_along_x): Reconstructing  $\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\}$ along the $x$-direction
#         1. [Step 3.8.2](#fluxes_x_dirn): Evaluating $\partial_{x}\boldsymbol{F}$
#         1. [Step 3.8.3](#reconstructing_vx_vy_by_along_y): Reconstructing $\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\}$ along the $y$-direction
#         1. [Step 3.8.4](#fluxes_y_dirn): Evaluating $\partial_{y}\boldsymbol{F}$
#         1. [Step 3.8.5](#rhs_az_no_gauge_terms): Evaluating $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$
#         1. [Step 3.8.6](#multiple_reconstructions): Multiple reconstructions
#         1. [Step 3.8.7](#fluxes_z_dirn): Evaluating $\partial_{z}\boldsymbol{F}$
#         1. [Step 3.8.8](#rhs_ax_no_gauge_terms): Evaluating $\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms}$
#         1. [Step 3.8.9](#rhs_ay_no_gauge_terms): Evaluating $\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms}$
#         1. [Step 3.8.10](#rhs_psi6phi_and_ai_gauge_terms): Evaluating $\partial_{t}\left[\psi^{6}\Phi\right]$ and $\left[\partial_{t}A_{i}\right]_{\rm gauge\ terms}$
# 1. [Step 4](#driver_evaluate_MHD_rhs__h): **The `driver_evaluate_MHD_rhs.h` header file**
# 1. [Step 5](#code_validation): **Code validation**
#     1. [Step 5.a](#code_validation_driver_evaluate_MHD_rhs__c): *`driver_evaluate_MHD_rhs.C`*
#     1. [Step 5.b](#code_validation_driver_evaluate_MHD_rhs__h): *`driver_evaluate_MHD_rhs.h`*
# 1. [Step 6](#latex_pdf_output): **Output this notebook to $\LaTeX$-formatted PDF file**

# <a id='src_dir'></a>
# 
# # Step 0: Source directory creation \[Back to [top](#toc)\]
# $$\label{src_dir}$$
# 
# 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.

# In[1]:


# Step 0: Creation of the IllinoisGRMHD source directory
# Step 0a: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os,sys
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)

# Step 0b: Load up cmdline_helper and create the directory
import cmdline_helper as cmd
IGM_src_dir_path = os.path.join("..","src")
cmd.mkdir(IGM_src_dir_path)

# Step 0c: Create the output file path
outfile_path__driver_evaluate_MHD_rhs__C = os.path.join(IGM_src_dir_path,"driver_evaluate_MHD_rhs.C")
outfile_path__driver_evaluate_MHD_rhs__h = os.path.join(IGM_src_dir_path,"driver_evaluate_MHD_rhs.h")


# <a id='introduction'></a>
# 
# # Step 1: Introduction \[Back to [top](#toc)\]
# $$\label{introduction}$$
# 
# We will start by creating the file `driver_evaluate_MHD_rhs.C` and writing down the preamble of the file, which contains useful references and information for the user.
# 
# We remind the reader of the "[generalized Lorenz gauge condition](https://arxiv.org/pdf/1207.3354.pdf)",
# 
# $$
# \nabla_{\mu}\mathcal{A^{\mu}} = \xi n_{\mu}\mathcal{A^{\mu}}\ ,
# $$
# 
# where $n_{\mu} = \left(\alpha,0,0,0\right)$ is the normal unit vector, $\mathcal{A}_{\mu}$ is the magnetic 4-vector potential, and $xi$ is a parameter with dimensions 1/Length, just like the $\eta$ parameter in the gamma-driving shift condition, so its value is chosen so that the CFL condition remains satisfied.

# In[2]:


get_ipython().run_cell_magic('writefile', '$outfile_path__driver_evaluate_MHD_rhs__C', '/*********************************************\n * Evaluate RHS of GRMHD & induction equations\n * (vector potential prescription), using the\n * generalized Lorenz gauge condition for the\n * EM gauge.\n *\n * Based originally on the Illinois GRMHD code,\n * written by Matt Duez, Yuk Tung Liu, and Branson\n * Stephens (original version), and then developed\n * primarily by Zachariah Etienne, Yuk Tung Liu,\n * and Vasileios Paschalidis.\n *\n * Rewritten for public release in 2013\n *      by Zachariah B. Etienne\n *\n * References:\n * Original unigrid GRMHD evolution prescription:\n *    http://arxiv.org/abs/astro-ph/0503420\n * Vector potential formulation in full GR:\n *    http://arxiv.org/abs/1007.2848\n * Improved EM gauge conditions for AMR grids:\n *    http://arxiv.org/abs/1110.4633\n * Generalized Lorenz gauge prescription:\n *    http://arxiv.org/abs/1207.3354\n *\n * Note that the Generalized Lorenz gauge strength\n *  parameter has units of 1/M, just like the \\eta\n *  parameter in the gamma-driving shift condition,\n *  so setting it too large will result in violation\n *  of the CFL condition.\n *\n * This version of PPM implements the standard\n * Colella & Woodward PPM, though modified as in GRHydro\n * to have 3 ghostzones instead of 4.\n *********************************************/\n')


# <a id='header_files'></a>
# 
# # Step 2: Load up necessary ETK and `IllinoisGRMHD` header files \[Back to [top](#toc)\]
# $$\label{header_files}$$
# 
# Here we load all necessary ETK and `IllinoisGRMHD` files, as well as some standard C++ libraries.

# In[3]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n#include "cctk.h"\n#include <cstdio>\n#include <cstdlib>\n#include <cmath>\n#include <sys/time.h>\n#include "cctk_Arguments.h"\n#include "cctk_Parameters.h"\n#include "IllinoisGRMHD_headers.h" /* Generic #define\'s and function prototypes */\n#include "driver_evaluate_MHD_rhs.h" /* Function prototypes for this file only */\n#include "IllinoisGRMHD_EoS_lowlevel_functs.C"\n#include "inlined_functions.C"\n')


# <a id='driver_mhd_rhs_function'></a>
# 
# # Step 3: The ` IllinoisGRMHD_driver_evaluate_MHD_rhs()` function \[Back to [top](#toc)\]
# $$\label{driver_mhd_rhs_function}$$
# 
# This is the basic function declaration. We set up basic ETK parameters and verify double precision is being used.

# In[4]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\nextern "C" void IllinoisGRMHD_driver_evaluate_MHD_rhs(CCTK_ARGUMENTS) {\n  DECLARE_CCTK_ARGUMENTS;\n  DECLARE_CCTK_PARAMETERS;\n  int levelnumber = GetRefinementLevel(cctkGH);\n\n  if(CCTK_Equals(verbose, "essential+iteration output")) {\n    CCTK_VInfo(CCTK_THORNSTRING,"***** Iter. # %d, Lev: %d, Integrating to time: %e *****",cctk_iteration,levelnumber,cctk_delta_time/cctk_levfac[0]+cctk_time);\n  }\n\n  if( sizeof(CCTK_REAL) < 8 ) CCTK_VError(VERR_DEF_PARAMS,"Error: IllinoisGRMHD assumes that CCTK_REAL is a double precision number. Setting otherwise will likely cause havoc with the conserv_to_prims solver.");\n\n  if(cctk_nghostzones[0]<3 || cctk_nghostzones[1]<3 || cctk_nghostzones[2]<3) { CCTK_VError(VERR_DEF_PARAMS,"ERROR. Need at least 3 ghostzones for IllinoisGRMHD evolutions."); }\n\n  CCTK_REAL dX[3] = { CCTK_DELTA_SPACE(0), CCTK_DELTA_SPACE(1), CCTK_DELTA_SPACE(2) };\n')


# <a id='eos_parameters'></a>
# 
# ## Step 3.1: Equation of state (EOS) parameters \[Back to [top](#toc)\]
# $$\label{eos_parameters}$$
# 
# Next, we set up the EOS struct, which is defined in the `IllinoisGRMHD_headers.h` header file. We set the following parameters:
# 
# * $\rm neos$: number of EOS (currently set to 1, which is a $\Gamma$-law EOS)
# * $\rm K\_poly$: this is the constant $\kappa$ from the polytropic EOS $P = \kappa \rho_{0}^{\Gamma}$
# * $\rm rho\_poly$: $\rho_{0}$, fluid rest-mass
# * $\rm P\_poly$: $P$, pressure
# * $\rm gamma\_th$: $\Gamma_{\rm th}$, the constant parameter which determines the conversion efficiency of kinetic to thermal energy at shocks
# * $\rm eps\_poly$: $\epsilon$, specific internal energy
# * $\rm k\_poly$: $\kappa$, polytropic EOS constant
# * $\rm gamma\_poly$: $\Gamma$, the $\Gamma$-law EOS constant

# In[5]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  /**********************************\n   * Piecewise Polytropic EOS Patch *\n   *   Setting up the EOS struct    *\n   **********************************/\n  /*\n   * The short piece of code below takes care\n   * of initializing the EOS parameters.\n   * Please refer to the "inlined_functions.C"\n   * source file for the documentation on the\n   * function.\n   */\n  eos_struct eos;\n  initialize_EOS_struct_from_input(eos);\n')


# <a id='set_pointers_grmhd_gfs'></a>
# 
# ## Step 3.2: Set pointers to GRMHD gridfunctions \[Back to [top](#toc)\]
# $$\label{set_pointers_grmhd_gfs}$$

# In[6]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', "\n\n  // in_prims,out_prims_r, and out_prims_l are arrays of pointers to the actual gridfunctions.\n  gf_and_gz_struct in_prims[MAXNUMVARS],out_prims_r[MAXNUMVARS],out_prims_l[MAXNUMVARS];\n  int which_prims_to_reconstruct[MAXNUMVARS],num_prims_to_reconstruct;\n\n  /* SET POINTERS TO GRMHD GRIDFUNCTIONS */\n  // The order here MATTERS, and must be consistent with the global variable declarations in\n  //   evaluate_MHD_rhs_headers.h (look for RHOB=0, etc.)\n  //   For example, in_prims[0] _must_ be rho_b.\n  int ww=0;\n  in_prims[ww].gf=rho_b; out_prims_r[ww].gf=rho_br; out_prims_l[ww].gf=rho_bl; ww++;\n  in_prims[ww].gf=P;     out_prims_r[ww].gf=Pr;     out_prims_l[ww].gf=Pl;     ww++;\n  in_prims[ww].gf=vx;    out_prims_r[ww].gf=vxr;    out_prims_l[ww].gf=vxl;    ww++;\n  in_prims[ww].gf=vy;    out_prims_r[ww].gf=vyr;    out_prims_l[ww].gf=vyl;    ww++;\n  in_prims[ww].gf=vz;    out_prims_r[ww].gf=vzr;    out_prims_l[ww].gf=vzl;    ww++;\n  in_prims[ww].gf=Bx;    out_prims_r[ww].gf=Bxr;    out_prims_l[ww].gf=Bxl;    ww++;\n  in_prims[ww].gf=By;    out_prims_r[ww].gf=Byr;    out_prims_l[ww].gf=Byl;    ww++;\n  in_prims[ww].gf=Bz;    out_prims_r[ww].gf=Bzr;    out_prims_l[ww].gf=Bzl;    ww++;\n  in_prims[ww].gf=Bx_stagger; out_prims_r[ww].gf=Bx_staggerr; out_prims_l[ww].gf=Bx_staggerl; ww++;\n  in_prims[ww].gf=By_stagger; out_prims_r[ww].gf=By_staggerr; out_prims_l[ww].gf=By_staggerl; ww++;\n  in_prims[ww].gf=Bz_stagger; out_prims_r[ww].gf=Bz_staggerr; out_prims_l[ww].gf=Bz_staggerl; ww++;\n  in_prims[ww].gf=vxr;    out_prims_r[ww].gf=vxrr;    out_prims_l[ww].gf=vxrl;    ww++;\n  in_prims[ww].gf=vyr;    out_prims_r[ww].gf=vyrr;    out_prims_l[ww].gf=vyrl;    ww++;\n  in_prims[ww].gf=vzr;    out_prims_r[ww].gf=vzrr;    out_prims_l[ww].gf=vzrl;    ww++;\n  in_prims[ww].gf=vxl;    out_prims_r[ww].gf=vxlr;    out_prims_l[ww].gf=vxll;    ww++;\n  in_prims[ww].gf=vyl;    out_prims_r[ww].gf=vylr;    out_prims_l[ww].gf=vyll;    ww++;\n  in_prims[ww].gf=vzl;    out_prims_r[ww].gf=vzlr;    out_prims_l[ww].gf=vzll;    ww++;\n\n  // Prims are defined AT ALL GRIDPOINTS, so we set the # of ghostzones to zero:\n  for(int i=0;i<MAXNUMVARS;i++) for(int j=1;j<=3;j++) { in_prims[i].gz_lo[j]=0; in_prims[i].gz_hi[j]=0; }\n  // Left/right variables are not yet defined, yet we set the # of gz's to zero by default:\n  for(int i=0;i<MAXNUMVARS;i++) for(int j=1;j<=3;j++) { out_prims_r[i].gz_lo[j]=0; out_prims_r[i].gz_hi[j]=0; }\n  for(int i=0;i<MAXNUMVARS;i++) for(int j=1;j<=3;j++) { out_prims_l[i].gz_lo[j]=0; out_prims_l[i].gz_hi[j]=0; }\n")


# <a id='admbase_to_bssnbase'></a>
# 
# ## Step 3.3: Convert ADM variables to BSSN variables \[Back to [top](#toc)\]
# $$\label{admbase_to_bssnbase}$$
# 
# We summarize here the algorithm of the `IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij()` function, which is explained in detail in [this tutorial module]() (<font color='red'>**Link not available yet - TODO**</font>):
# 
# * First, $\gamma\equiv \det\left(\gamma_{ij}\right)$, where $\gamma_{ij}$ is the physical spatial metric, is evaluated.
# * Then, $\phi$, the conformal factor, is computed via the relation $\phi = \frac{1}{12}\log\gamma$.
# * Next, we compute $e^{-4\phi}$.
# * Then, the conformal metric, $\bar{\gamma}_{ij} = e^{-4\phi}\gamma_{ij}$, is computed.
# * Next, the condition $\bar\gamma = 1$ is enforced, by first computing $\gamma$ then performing $\bar\gamma_{ij}\to\left(\frac{1}{\bar\gamma}\right)^{1/3}\bar\gamma_{ij}$.
# * Then, $\gamma_{ij}$ is computed from $\bar\gamma_{ij}$ *after* the condition $\bar\gamma = 1$ is enforced via the inverse relation $\gamma_{ij} = e^{4\phi}\bar\gamma_{ij}$.
# * Finally, we compute the inverse conformal metric $\bar\gamma^{ij}$.
# 
# **A note on notation:** In the C code, we have the following identifications between the quantities described above and the C variables:
# 
# *Input and temporary variables:*
# * $\gamma_{ij} := {\rm gij\_physL}$
# * $\det(\gamma_{ij}) := {\rm gijdet}$
# * $\phi := {\rm phiL}$
# * $\psi \equiv e^{\phi} := {\rm psiL}$
# * $\bar\gamma_{ij} := {\rm gtijL}$ (the "t" stands for the notation where the conformal metric is written as $\tilde\gamma_{ij}$ instead of $\bar\gamma_{ij}$. In our discussion we use the latter to keep our notation consistent with other NRPy notebooks).
# * $\det(\bar\gamma_{ij}) := {\rm gtijdet}$
# * $\left(\frac{1}{\bar\gamma}\right)^{1/3} := {\rm gtijdet\_Fm1o3}$
# 
# *Output/gridfunction variables:*
# * $\gamma_{ij} := {\rm gij}$ (Physical metric)
# * $\bar\gamma_{ij} := {\rm gtij}$ (Conformal metric)
# * $\bar\gamma^{ij} := {\rm gtupij}$ (Inverse conformal metric)
# * $\phi := {\rm phi}$
# * $\psi := {\rm psi}$

# In[7]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\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')


# <a id='pointers_metric_tmunu_gfs'></a>
# 
# ## Step 3.4: Setting up pointers to the metric and stress-energy tensor gridfunctions \[Back to [top](#toc)\]
# $$\label{pointers_metric_tmunu_gfs}$$

# In[8]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  /* SET POINTERS TO METRIC GRIDFUNCTIONS */\n  CCTK_REAL *metric[NUMVARS_FOR_METRIC_FACEVALS]; // "metric" here is array of pointers to the actual gridfunctions.\n  ww=0;\n  metric[ww]=phi_bssn;ww++;\n  metric[ww]=psi_bssn;ww++;\n  metric[ww]=gtxx;    ww++;\n  metric[ww]=gtxy;    ww++;\n  metric[ww]=gtxz;    ww++;\n  metric[ww]=gtyy;    ww++;\n  metric[ww]=gtyz;    ww++;\n  metric[ww]=gtzz;    ww++;\n  metric[ww]=lapm1;   ww++;\n  metric[ww]=betax;   ww++;\n  metric[ww]=betay;   ww++;\n  metric[ww]=betaz;   ww++;\n  metric[ww]=gtupxx;  ww++;\n  metric[ww]=gtupyy;  ww++;\n  metric[ww]=gtupzz;  ww++;\n\n  /* SET POINTERS TO STRESS-ENERGY TENSOR GRIDFUNCTIONS */\n  CCTK_REAL *TUPmunu[10];// "TUPmunu" here is array of pointers to the actual gridfunctions.\n  ww=0;\n  TUPmunu[ww]=TUPtt; ww++;\n  TUPmunu[ww]=TUPtx; ww++;\n  TUPmunu[ww]=TUPty; ww++;\n  TUPmunu[ww]=TUPtz; ww++;\n  TUPmunu[ww]=TUPxx; ww++;\n  TUPmunu[ww]=TUPxy; ww++;\n  TUPmunu[ww]=TUPxz; ww++;\n  TUPmunu[ww]=TUPyy; ww++;\n  TUPmunu[ww]=TUPyz; ww++;\n  TUPmunu[ww]=TUPzz; ww++;\n')


# <a id='initialize_rhss'></a>
# 
# ## Step 3.5: Initialization of the RHS variables \[Back to [top](#toc)\]
# $$\label{initialize_rhss}$$

# In[9]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  // 1) First initialize {rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs} to zero\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        Ax_rhs[index]=0.0;\n        Ay_rhs[index]=0.0;\n        Az_rhs[index]=0.0;\n        psi6phi_rhs[index]=0.0;\n\n        tau_rhs[index]=0.0;\n        rho_star_rhs[index]=0.0;\n        st_x_rhs[index]=0.0;\n        st_y_rhs[index]=0.0;\n        st_z_rhs[index]=0.0;\n\n        //if(i==17 && j==19 && k==26) CCTK_VInfo(CCTK_THORNSTRING,"CONSSS: %.15e %.15e %.15e %.15e %.15e | %.15e",rho_star[index],mhd_st_x[index],mhd_st_y[index],mhd_st_z[index],tau[index],P[index]);\n      }\n')


# <a id='tau_rhs_ext_curv_and_tupmunu'></a>
# 
# ## Step 3.6: Compute extrinsic curvature terms from the RHS of $\partial_{t}\tilde\tau$ and $T^{\mu\nu}$ \[Back to [top](#toc)\]
# $$\label{tau_rhs_ext_curv_and_tupmunu}$$

# In[10]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n  // Here, we:\n  // 1) Compute tau_rhs extrinsic curvature terms, and\n  // 2) Compute TUPmunu.\n  // This function is housed in the file: "compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C"\n  compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu(cctkGH,cctk_lsh,cctk_nghostzones,dX,metric,in_prims,TUPmunu,\n                                                        eos, Gamma_th,\n                                                        gtupxy,gtupxz,gtupyz,\n                                                        kxx,kxy,kxz,kyy,kyz,kzz,\n                                                        tau_rhs);\n')


# <a id='computing_ftilde'></a>
# 
# ## Step 3.7: Computing ${\rm ftilde}$ \[Back to [top](#toc)\]
# $$\label{computing_ftilde}$$
# 
# This is part of the flattening scheme of the PPM algorithm. The main reference to look at is [Colella & Woodward (1983)](https://crd.lbl.gov/assets/pubs_presos/AMCS/ANAG/A141984.pdf). The equations implemented can be found in Appendix A (particularly Eqs. (A.1) and (A.2)), while the flattening method is introduced and discussed in Section 4. More will follow when we talk about the `reconstruct_set_of_prims_PPM.C` file of `IllinoisGRMHD`.

# In[11]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n  int flux_dirn;\n  flux_dirn=1;\n  // First compute ftilde, which is used for flattening left and right face values\n  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"\n  ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);\n')


# <a id='rhs_mhd_and_a_i'></a>
# 
# ## Step 3.8: The RHSs of $\rho_{\star}$, $\tilde\tau$, $\tilde{S}_{i}$, and $A_{i}$ \[Back to [top](#toc)\]
# $$\label{rhs_mhd_and_a_i}$$
# 
# This part of the code evaluates the RHSs of $\rho_{\star}$, $\tilde\tau$, and $\tilde{S}_{i}$, i.e.
# 
# $$
# \partial_{t}
# \begin{bmatrix}
# \rho_{\star}\\
# \tilde\tau\\
# \tilde{S}_{i}
# \end{bmatrix}
# =
# -\partial_{j}
# \underbrace{\begin{bmatrix}
# \rho_{\star}v^{j}\\
# \alpha^{2}\sqrt{\gamma}T^{0j} - \rho_{\star}v^{j}\\
# \alpha\sqrt{\gamma}T^{j}_{\ i}
# \end{bmatrix}}_{\rm Flux\ terms}
# +
# \underbrace{\begin{bmatrix}
# 0\\
# s\\
# \frac{1}{2}\alpha\sqrt{\gamma}T^{\alpha\beta}\partial_{i}g_{\alpha\beta}
# \end{bmatrix}}_{\rm Source\ terms}\ .
# $$
# 
# At the same time, we are also interested in evaluating the RHS of the evolution equation for $A_{i}$, namely
# 
# $$
# \partial_{t}A_{i} = \epsilon_{ijk}v^{j}\tilde{B}^{k} - \partial_{i}\left(\alpha\Phi - \beta^{j}A_{j}\right) = \psi^{6}\epsilon_{ijk}v^{j}B^{k} - \underbrace{\partial_{i}\left(\alpha\Phi - \beta^{j}A_{j}\right)}_{\rm Gauge\ terms}
# $$
# 
# The following summary greatly oversimplifies what the code below does, but it is enough for the user to understand the purpose of the algorithm:
# 
# 1. Compute $\partial_{x}\boldsymbol{F}$, then $\partial_{y}\boldsymbol{F}$, and finally $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$
# 2. Compute $\partial_{y}\boldsymbol{F}$, then $\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms}$
# 3. Compute $\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms}$
# 4. Add gauge terms to $\partial_{t}A_{i}$
# 
# Now, in between every step of the summary above, care must be taken to evaluate the gridfunctions at the appropriate gridpoints (see table below for the location of each variable in the computational grid).
# 
# <a id='table_staggerings'></a>
# 
# |                           Variable(s)                            | Gridpoint location in the computational grid  |
# |------------------------------------------------------------------|-----------------------------------------------|
# | Metric terms, $\vec{P}$, $\rho_*$, $\tilde{S}_i$, $\tilde{\tau}$ |                   $(i,j,k)$                   |
# |                       $B^x$, $\tilde{B}^x$                       |             $(i+\frac{1}{2},j,k)$             |
# |                       $B^y$, $\tilde{B}^y$                       |             $(i,j+\frac{1}{2},k)$             |
# |                       $B^z$, $\tilde{B}^z$                       |             $(i,j,k+\frac{1}{2})$             |
# |                              $A_x$                               |       $(i,j+\frac{1}{2},k+\frac{1}{2})$       |
# |                              $A_y$                               |       $(i+\frac{1}{2},j,k+\frac{1}{2})$       |
# |                              $A_z$                               |       $(i+\frac{1}{2},j+\frac{1}{2},k)$       |
# |                       $\sqrt{\gamma}\Phi$                        | $(i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2})$ |
# 
# $$\label{table_staggerings}$$
# 
# For example, we know that
# 
# $$
# \left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms} \equiv \psi^{6}\left(v^{x}B^{y} - v^{y}B^{x}\right)\ .
# $$
# 
# But once we evaluate $v^{x}$ and $v^{y}$, we know them at the point $(i,j,k)$. Similarly, the gridfunction $B^{x}$ is known at $(i+\frac{1}{2},j,k)$, while $B^{y}$ is known at $(i,j+\frac{1}{2},k)$. This means that we are not able to immediately evaluate the equation above, since determining $A_{z}$ at $(i+\frac{1}{2},j+\frac{1}{2},k)$ requires knowing $\left\{v^{x},v^{y},B^{x},B^{y}\right\}$ at $(i+\frac{1}{2},j+\frac{1}{2},k)$ as well. To this end, we reconstruct the variables $\left\{v^{x},v^{y},B^{x},B^{y}\right\}$ using the PPM method at the desired staggered point. An analogous procedure is required in order to determine the RHS of $\partial_{t}A_{x}$ and $\partial_{t}A_{y}$.

# <a id='reconstructing_vx_vy_by_along_x'></a>
# 
# ### Step 3.8.1: Reconstructing  $\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\}$ along the $x$-direction \[Back to [top](#toc)\]
# $$\label{reconstructing_vx_vy_by_along_x}$$
# 
# We want to evaluate $\partial_{x}\boldsymbol{F}$. It is important that we keep in the back of our minds our intention of evaluating the RHS of $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$ as well, since then we can reconstruct $\left\{v^{x},v^{y},B^{x},B^{y}\right\}$ cleverly, as we need them at the same gridpoint as $A_{z}$ (see the table at the end of [step 3.8](#table_staggerings)).
# 
# We start by reconstructing $\left\{\rho_{0},P,v^{i},B^{i}, B^{y}_{\rm stagger}\right\}$ in the $x$-direction, keeping in mind that after the reconstruction we will know:
# 
# 1. The flux variables at $\left(i-\frac{1}{2},j,k\right)$, so that we can evaluate $\partial_{x}\boldsymbol{F}_{i,j,k}=dx^{-1}\left(\boldsymbol{F}_{i+1/2,j,k}-\boldsymbol{F}_{i-1/2,j,k}\right)$
# 1. The velocities $\left\{v^{x},v^{y}\right\}$ at $\left(i-\frac{1}{2},j,k\right)$
# 1. The staggered value $B^{y}_{\rm stagger}$ at $\left(i-\frac{1}{2},j+\frac{1}{2},k\right)$

# In[12]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n\n  /* There are two stories going on here:\n   * 1) Computation of \\partial_x on RHS of \\partial_t {rho_star,tau,mhd_st_{x,y,z}},\n   *    via PPM reconstruction onto (i-1/2,j,k), so that\n   *    \\partial_x F = [ F(i+1/2,j,k) - F(i-1/2,j,k) ] / dx\n   * 2) Computation of \\partial_t A_i, where A_i are *staggered* gridfunctions,\n   *    where A_x is defined at (i,j+1/2,k+1/2), A_y at (i+1/2,j,k+1/2), etc.\n   *    Ai_rhs = \\partial_t A_i = \\epsilon_{ijk} \\psi^{6} v^j B^k,\n   *    where \\epsilon_{ijk} is the flat-space antisymmetric operator.\n   * 2A) Az_rhs is defined at (i+1/2,j+1/2,k), and it depends on {Bx,By,vx,vy},\n   *     so the trick is to reconstruct {Bx,By,vx,vy} cleverly to get to these\n   *     staggered points. For example:\n   * 2Aa) vx and vy are at (i,j,k), and we reconstruct them to (i-1/2,j,k) below. After\n   *      this, we\'ll reconstruct again in the y-dir\'n to get {vx,vy} at (i-1/2,j-1/2,k)\n   * 2Ab) By_stagger is at (i,j+1/2,k), and we reconstruct below to (i-1/2,j+1/2,k). */\n  ww=0;\n  which_prims_to_reconstruct[ww]=RHOB;      ww++;\n  which_prims_to_reconstruct[ww]=PRESSURE;  ww++;\n  which_prims_to_reconstruct[ww]=VX;        ww++;\n  which_prims_to_reconstruct[ww]=VY;        ww++;\n  which_prims_to_reconstruct[ww]=VZ;        ww++;\n  //which_prims_to_reconstruct[ww]=BX_CENTER; ww++;\n  which_prims_to_reconstruct[ww]=BY_CENTER; ww++;\n  which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;\n  which_prims_to_reconstruct[ww]=BY_STAGGER;ww++;\n  num_prims_to_reconstruct=ww;\n  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"\n  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);\n')


# <a id='fluxes_x_dirn'></a>
# 
# ### Step 3.8.2: Evaluating $\partial_{x}\boldsymbol{F}$ \[Back to [top](#toc)\]
# $$\label{fluxes_x_dirn}$$
# 
# Next, we set the face values of $B^{x}$ (which are needed for the computation of MHD flux terms) by making them consistent with $B^{x}_{\rm stagger}$.
# 
# After that, we evaluate $\partial_{x}\boldsymbol{F}$ and add it to the RHS of $\partial_{t}\left(\rho_{\star},\tilde\tau,\tilde{S}_{i}\right)$. It is important to notice that, as we mentioned, $A_{z}$ is defined at $\left(i+\frac{1}{2},j+\frac{1}{2},k\right)$, but other functions, like $v^{x}$ and $v^{y}$, are now known only at $\left(i-\frac{1}{2},j-\frac{1}{2},k\right)$. The function `add_fluxes_and_source_terms_to_hydro_rhss()` below takes care of this, and we will study the process in more detail when we look at the `add_fluxes_and_source_terms_to_hydro_rhss.C` file of `IllinoisGRMHD`.

# In[13]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n  //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a).\n  //   Instead of reconstructing, we simply set B^x face values to be consistent with BX_STAGGER.\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), indexim1=CCTK_GFINDEX3D(cctkGH,i-1+(i==0),j,k); /* indexim1=0 when i=0 */\n        out_prims_r[BX_CENTER].gf[index]=out_prims_l[BX_CENTER].gf[index]=in_prims[BX_STAGGER].gf[indexim1]; }\n  // Then add fluxes to RHS for hydro variables {rho_b,P,vx,vy,vz}:\n  // This function is housed in the file: "add_fluxes_and_source_terms_to_hydro_rhss.C"\n  add_fluxes_and_source_terms_to_hydro_rhss(flux_dirn,cctkGH,cctk_lsh,cctk_nghostzones,dX,   metric,in_prims,TUPmunu,\n                                            num_prims_to_reconstruct,out_prims_r,out_prims_l,eos,\n                                            cmax_x,cmin_x,\n                                            rho_star_flux,tau_flux,st_x_flux,st_y_flux,st_z_flux,\n                                            rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs);\n')


# <a id='reconstructing_vx_vy_by_along_y'></a>
# 
# ### Step 3.8.3: Reconstructing $\left\{v^{x}, v^{y}, B^{y}_{\rm stagger}\right\}$ along the $y$-direction \[Back to [top](#toc)\]
# $$\label{reconstructing_vx_vy_by_along_y}$$
# 
# We want to evaluate $\partial_{y}\boldsymbol{F}$. At this point we must remember that $v^{x}$ and $v^{y}$ have already been reconstruct along the $x$-direction and are now known at $\left(i-\frac{1}{2},j,k\right)$. Our goal is to reconstruct these quantities at $\left(i+\frac{1}{2},j+\frac{1}{2},k\right)$.
# 
# We then reconstruct $\left\{\rho_{0},P,v^{i},B^{i}, B^{i}_{\rm stagger}\right\}$ in the $y$-direction, keeping in mind that after the reconstruction we will know:
# 
# 1. The flux variables at $\left(i-\frac{1}{2},j-\frac{1}{2},k\right)$, so that we can evaluate $\partial_{y}\boldsymbol{F}_{i,j,k}=dy^{-1}\left(\boldsymbol{F}_{i,j+1/2,k}-\boldsymbol{F}_{i,j-1/2,k}\right)$
# 1. The velocities $\left\{v^{x},v^{y}\right\}$ at $\left(i-\frac{1}{2},j-\frac{1}{2},k\right)$
# 1. The staggered value $B^{x}_{\rm stagger}$ at $\left(i+\frac{1}{2},j-\frac{1}{2},k\right)$
# 1. The staggered value $B^{y}_{\rm stagger}$ at $\left(i-\frac{1}{2},j+\frac{1}{2},k\right)$
# 1. The staggered value $B^{z}_{\rm stagger}$ at $\left(i,j-\frac{1}{2},k+\frac{1}{2}\right)$

# In[14]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  // Note that we have already reconstructed vx and vy along the x-direction,\n  //   at (i-1/2,j,k). That result is stored in v{x,y}{r,l}.  Bx_stagger data\n  //   are defined at (i+1/2,j,k).\n  // Next goal: reconstruct Bx, vx and vy at (i+1/2,j+1/2,k).\n  flux_dirn=2;\n  // First compute ftilde, which is used for flattening left and right face values\n  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"\n  ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);\n\n  // in_prims[{VXR,VXL,VYR,VYL}].gz_{lo,hi} ghostzones are set to all zeros, which\n  //    is incorrect. We fix this below.\n  // [Note that this is a cheap operation, copying only 8 integers and a pointer.]\n  in_prims[VXR]=out_prims_r[VX];\n  in_prims[VXL]=out_prims_l[VX];\n  in_prims[VYR]=out_prims_r[VY];\n  in_prims[VYL]=out_prims_l[VY];\n\n  /* There are two stories going on here:\n   * 1) Computation of \\partial_y on RHS of \\partial_t {rho_star,tau,mhd_st_{x,y,z}},\n   *    via PPM reconstruction onto (i,j-1/2,k), so that\n   *    \\partial_y F = [ F(i,j+1/2,k) - F(i,j-1/2,k) ] / dy\n   * 2) Computation of \\partial_t A_i, where A_i are *staggered* gridfunctions,\n   *    where A_x is defined at (i,j+1/2,k+1/2), A_y at (i+1/2,j,k+1/2), etc.\n   *    Ai_rhs = \\partial_t A_i = \\epsilon_{ijk} \\psi^{6} v^j B^k,\n   *    where \\epsilon_{ijk} is the flat-space antisymmetric operator.\n   * 2A) Az_rhs is defined at (i+1/2,j+1/2,k), and it depends on {Bx,By,vx,vy},\n   *     so the trick is to reconstruct {Bx,By,vx,vy} cleverly to get to these\n   *     staggered points. For example:\n   * 2Aa) VXR = [right-face of vx reconstructed along x-direction above] is at (i-1/2,j,k),\n   *      and we reconstruct it to (i-1/2,j-1/2,k) below. Similarly for {VXL,VYR,VYL}\n   * 2Ab) Bx_stagger is at (i+1/2,j,k), and we reconstruct to (i+1/2,j-1/2,k) below\n   * 2Ac) By_stagger is at (i-1/2,j+1/2,k) already for Az_rhs, from the previous step.\n   * 2B) Ax_rhs is defined at (i,j+1/2,k+1/2), and it depends on {By,Bz,vy,vz}.\n   *     Again the trick is to reconstruct these onto these staggered points.\n   * 2Ba) Bz_stagger is at (i,j,k+1/2), and we reconstruct to (i,j-1/2,k+1/2) below */\n  ww=0;\n  // NOTE! The order of variable reconstruction is important here,\n  //   as we don\'t want to overwrite {vxr,vxl,vyr,vyl}!\n  which_prims_to_reconstruct[ww]=VXR;       ww++;\n  which_prims_to_reconstruct[ww]=VYR;       ww++;\n  which_prims_to_reconstruct[ww]=VXL;       ww++;\n  which_prims_to_reconstruct[ww]=VYL;       ww++;\n  num_prims_to_reconstruct=ww;\n  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"\n  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);\n  ww=0;\n  // Reconstruct other primitives last!\n  which_prims_to_reconstruct[ww]=RHOB;      ww++;\n  which_prims_to_reconstruct[ww]=PRESSURE;  ww++;\n  which_prims_to_reconstruct[ww]=VX;        ww++;\n  which_prims_to_reconstruct[ww]=VY;        ww++;\n  which_prims_to_reconstruct[ww]=VZ;        ww++;\n  which_prims_to_reconstruct[ww]=BX_CENTER; ww++;\n  //which_prims_to_reconstruct[ww]=BY_CENTER; ww++;\n  which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;\n  which_prims_to_reconstruct[ww]=BX_STAGGER;ww++;\n  which_prims_to_reconstruct[ww]=BZ_STAGGER;ww++;\n  num_prims_to_reconstruct=ww;\n  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"\n  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);\n')


# <a id='fluxes_y_dirn'></a>
# 
# ### Step 3.8.4: Evaluating $\partial_{y}\boldsymbol{F}$ \[Back to [top](#toc)\]
# $$\label{fluxes_y_dirn}$$

# In[15]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n  //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a).\n  //   Instead of reconstructing, we simply set B^y face values to be consistent with BY_STAGGER.\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), indexjm1=CCTK_GFINDEX3D(cctkGH,i,j-1+(j==0),k); /* indexjm1=0 when j=0 */\n        out_prims_r[BY_CENTER].gf[index]=out_prims_l[BY_CENTER].gf[index]=in_prims[BY_STAGGER].gf[indexjm1]; }\n  // Then add fluxes to RHS for hydro variables {rho_b,P,vx,vy,vz}:\n  // This function is housed in the file: "add_fluxes_and_source_terms_to_hydro_rhss.C"\n  add_fluxes_and_source_terms_to_hydro_rhss(flux_dirn,cctkGH,cctk_lsh,cctk_nghostzones,dX,   metric,in_prims,TUPmunu,\n                                            num_prims_to_reconstruct,out_prims_r,out_prims_l,eos,\n                                            cmax_y,cmin_y,\n                                            rho_star_flux,tau_flux,st_x_flux,st_y_flux,st_z_flux,\n                                            rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs);\n')


# <a id='rhs_az_no_gauge_terms'></a>
# 
# ### Step 3.8.5: Evaluating $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$ \[Back to [top](#toc)\]
# $$\label{rhs_az_no_gauge_terms}$$
# 
# As a friendly reminder, we summarize the known gridpoint location of the needed gridfunctions here:
# 
# | Staggered and unstaggered variables                 | Gridpoint location at which the variable is known  |
# |-----------------------------------------------------|----------------------------------------------------|
# | $\left(v^{x}\right)_{r,l},\left(v^{y}\right)_{r,l}$ |          $(i-\frac{1}{2},j-\frac{1}{2},k)$         |
# | $\left(B^{z}_{\rm stagger}\right)_{r,l}$            |          $(i+\frac{1}{2},j-\frac{1}{2},k)$         |
# | $\left(B^{y}_{\rm stagger}\right)_{r,l}$            |          $(i-\frac{1}{2},j+\frac{1}{2},k)$         |
# |                  $\phi$                             |                      $(i,j,k)$                     |
# 
# We start by interpolating $\phi$ to $\left(i+\frac{1}{2},j,k\right)$, followed by a second interpolation so that $\phi$ is known at $\left(i+\frac{1}{2},j+\frac{1}{2},k\right)$.

# In[16]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  /*****************************************\n   * COMPUTING RHS OF A_z, BOOKKEEPING NOTE:\n   * We want to compute\n   * \\partial_t A_z - [gauge terms] = \\psi^{6} (v^x B^y - v^y B^x).\n   * A_z is defined at (i+1/2,j+1/2,k).\n   * ==========================\n   * Where defined  | Variables\n   * (i-1/2,j-1/2,k)| {vxrr,vxrl,vxlr,vxll,vyrr,vyrl,vylr,vyll}\n   * (i+1/2,j-1/2,k)| {Bx_stagger_r,Bx_stagger_l} (see Table 1 in arXiv:1007.2848)\n   * (i-1/2,j+1/2,k)| {By_stagger_r,By_stagger_l} (see Table 1 in arXiv:1007.2848)\n   * (i,j,k)        | {phi}\n   * ==========================\n   ******************************************/\n  // Interpolates to i+1/2\n#define IPH(METRICm1,METRICp0,METRICp1,METRICp2) (-0.0625*((METRICm1) + (METRICp2)) + 0.5625*((METRICp0) + (METRICp1)))\n  // Next compute phi at (i+1/2,j+1/2,k):\n#pragma omp parallel for\n  for(int k=0;k<cctk_lsh[2];k++) for(int j=1;j<cctk_lsh[1]-2;j++) for(int i=1;i<cctk_lsh[0]-2;i++) {\n        temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=\n          IPH(IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j-1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j-1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j-1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j-1,k)]),\n              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j  ,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j  ,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j  ,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j  ,k)]),\n              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j+1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j+1,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j+1,k)]),\n              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j+2,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+2,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j+2,k)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j+2,k)]));\n      }\n')


# Then we update the RHS of $\left[\partial_{t}A_{z}\right]_{\rm no\ gauge\ terms}$. Keep in mind that the function `A_i_rhs_no_gauge_terms()` takes care of determining $\left\{v^{i},B^{i},B^{i}_{\rm stagger}\right\}$ at $\left(i+\frac{1}{2},j+\frac{1}{2},k\right)$. We will look at it in more detail when we see the `A_i_rhs_no_gauge_terms.C` file from `IllinoisGRMHD`.

# In[17]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  int A_directionz=3;\n  A_i_rhs_no_gauge_terms(A_directionz,cctkGH,cctk_lsh,cctk_nghostzones,out_prims_r,out_prims_l,temporary,cmax_x,cmin_x,cmax_y,cmin_y, Az_rhs);\n')


# <a id='multiple_reconstructions'></a>
# 
# ### Step 3.8.6: Multiple reconstructions \[Back to [top](#toc)\]
# $$\label{multiple_reconstructions}$$
# 
# 1. $\left\{\rho_{\star},P,v^{i}, B^{i}\right\}$ from $\left(i,j,k\right)$ to $\left(i,j,k-\frac{1}{2}\right)$
# 1. $\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms}$ is defined at $\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)$
#     1. $\left(v^{y}\right)_{r,l}$ and $\left(v^{z}\right)_{r,l}$ are at $\left(i,j-\frac{1}{2},k\right)$, so we reconstruct them to $\left(i,j-\frac{1}{2},k-\frac{1}{2}\right)$
#     1. $\left(B^{z}_{\rm stagger}\right)_{r,l}$ is already known at $\left(i,j-\frac{1}{2},k+\frac{1}{2}\right)$
#     1. $B^{y}_{\rm stagger}$ is at $\left(i,j+\frac{1}{2},k\right)$, so we reconstruct it to $\left(i,j+\frac{1}{2},k-\frac{1}{2}\right)$
# 1. $\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms}$ is defined at $\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)$
#     1. $v^{x}$ and $v^{z}$ are reconstructed to $\left(i,j,k-\frac{1}{2}\right)$. We'll reconstruct them to $\left(i,j-\frac{1}{2},k-\frac{1}{2}\right)$ later
#     1. $B^{z}_{\rm stagger}$ is already known at $\left(i,j,k+\frac{1}{2}\right)$. We'll reconstruct them to $\left(i-\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}\right)$ later
#     1. $B^{x}_{\rm stagger}$ is at $\left(i,j+\frac{1}{2},k\right)$, so we reconstruct it to $\left(i,j+\frac{1}{2},k-\frac{1}{2}\right)$

# In[18]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  // in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not correct, so we fix\n  //    this below.\n  // [Note that this is a cheap operation, copying only 8 integers and a pointer.]\n  in_prims[VYR]=out_prims_r[VY];\n  in_prims[VYL]=out_prims_l[VY];\n  in_prims[VZR]=out_prims_r[VZ];\n  in_prims[VZL]=out_prims_l[VZ];\n\n  flux_dirn=3;\n  // First compute ftilde, which is used for flattening left and right face values\n  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"\n  ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);\n\n  /* There are two stories going on here:\n   * 1) Single reconstruction to (i,j,k-1/2) for {rho,P,vx,vy,vz,Bx,By,Bz} to compute\n   *    z-dir\'n advection terms in \\partial_t {rho_star,tau,mhd_st_{x,y,z}} at (i,j,k)\n   * 2) Multiple reconstructions for *staggered* gridfunctions A_i:\n   *    Ai_rhs = \\partial_t A_i = \\epsilon_{ijk} \\psi^{6} v^j B^k,\n   *    where \\epsilon_{ijk} is the flat-space antisymmetric operator.\n   * 2A) Ax_rhs is defined at (i,j+1/2,k+1/2), depends on v{y,z} and B{y,z}\n   * 2Aa) v{y,z}{r,l} are at (i,j-1/2,k), so we reconstruct here to (i,j-1/2,k-1/2)\n   * 2Ab) Bz_stagger{r,l} are at (i,j-1/2,k+1/2) already.\n   * 2Ac) By_stagger is at (i,j+1/2,k), and below we reconstruct its value at (i,j+1/2,k-1/2)\n   * 2B) Ay_rhs is defined at (i+1/2,j,k+1/2), depends on v{z,x} and B{z,x}.\n   * 2Ba) v{x,z} are reconstructed to (i,j,k-1/2). Later we\'ll reconstruct again to (i-1/2,j,k-1/2).\n   * 2Bb) Bz_stagger is at (i,j,k+1/2). Later we will reconstruct to (i-1/2,j,k+1/2).\n   * 2Bc) Bx_stagger is at (i+1/2,j,k), and below we reconstruct its value at (i+1/2,j,k-1/2)\n   */\n  ww=0;\n  // NOTE! The order of variable reconstruction is important here,\n  //   as we don\'t want to overwrite {vxr,vxl,vyr,vyl}!\n  which_prims_to_reconstruct[ww]=VYR;       ww++;\n  which_prims_to_reconstruct[ww]=VZR;       ww++;\n  which_prims_to_reconstruct[ww]=VYL;       ww++;\n  which_prims_to_reconstruct[ww]=VZL;       ww++;\n  num_prims_to_reconstruct=ww;\n  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"\n  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);\n  // Reconstruct other primitives last!\n  ww=0;\n  which_prims_to_reconstruct[ww]=RHOB;      ww++;\n  which_prims_to_reconstruct[ww]=PRESSURE;  ww++;\n  which_prims_to_reconstruct[ww]=VX;        ww++;\n  which_prims_to_reconstruct[ww]=VY;        ww++;\n  which_prims_to_reconstruct[ww]=VZ;        ww++;\n  which_prims_to_reconstruct[ww]=BX_CENTER; ww++;\n  which_prims_to_reconstruct[ww]=BY_CENTER; ww++;\n  //which_prims_to_reconstruct[ww]=BZ_CENTER; ww++;\n  which_prims_to_reconstruct[ww]=BX_STAGGER; ww++;\n  which_prims_to_reconstruct[ww]=BY_STAGGER; ww++;\n  num_prims_to_reconstruct=ww;\n  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"\n  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);\n  //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a).\n  //   Instead of reconstructing, we simply set B^z face values to be consistent with BZ_STAGGER.\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), indexkm1=CCTK_GFINDEX3D(cctkGH,i,j,k-1+(k==0)); /* indexkm1=0 when k=0 */\n        out_prims_r[BZ_CENTER].gf[index]=out_prims_l[BZ_CENTER].gf[index]=in_prims[BZ_STAGGER].gf[indexkm1]; }\n')


# <a id='fluxes_z_dirn'></a>
# 
# ### Step 3.8.7: Evaluating $\partial_{z}\boldsymbol{F}$ \[Back to [top](#toc)\]
# $$\label{fluxes_z_dirn}$$

# In[19]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n  // Then add fluxes to RHS for hydro variables {rho_b,P,vx,vy,vz}:\n  // This function is housed in the file: "add_fluxes_and_source_terms_to_hydro_rhss.C"\n  add_fluxes_and_source_terms_to_hydro_rhss(flux_dirn,cctkGH,cctk_lsh,cctk_nghostzones,dX,   metric,in_prims,TUPmunu,\n                                            num_prims_to_reconstruct,out_prims_r,out_prims_l,eos,\n                                            cmax_z,cmin_z,\n                                            rho_star_flux,tau_flux,st_x_flux,st_y_flux,st_z_flux,\n                                            rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs);\n\n  // in_prims[{VYR,VYL,VZR,VZL}].gz_{lo,hi} ghostzones are not set correcty.\n  //    We fix this below.\n  // [Note that this is a cheap operation, copying only 8 integers and a pointer.]\n  in_prims[VXR]=out_prims_r[VX];\n  in_prims[VZR]=out_prims_r[VZ];\n  in_prims[VXL]=out_prims_l[VX];\n  in_prims[VZL]=out_prims_l[VZ];\n  // FIXME: lines above seem to be inconsistent with lines below.... Possible bug, not major enough to affect evolutions though.\n  in_prims[VZR].gz_lo[1]=in_prims[VZR].gz_hi[1]=0;\n  in_prims[VXR].gz_lo[1]=in_prims[VXR].gz_hi[1]=0;\n  in_prims[VZL].gz_lo[1]=in_prims[VZL].gz_hi[1]=0;\n  in_prims[VXL].gz_lo[1]=in_prims[VXL].gz_hi[1]=0;\n')


# <a id='rhs_ax_no_gauge_terms'></a>
# 
# ### Step 3.8.8: Evaluating $\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms}$ \[Back to [top](#toc)\]
# $$\label{rhs_ax_no_gauge_terms}$$
# 
# As a friendly reminder, we summarize the known gridpoint location of the needed gridfunctions here:
# 
# | Staggered and unstaggered variables                 | Gridpoint location at which the variable is known  |
# |-----------------------------------------------------|----------------------------------------------------|
# | $\left(v^{y}\right)_{r,l},\left(v^{z}\right)_{r,l}$ |          $(i,j-\frac{1}{2},k-\frac{1}{2})$         |
# | $\left(B^{y}_{\rm stagger}\right)_{r,l}$            |          $(i,j+\frac{1}{2},k-\frac{1}{2})$         |
# | $\left(B^{z}_{\rm stagger}\right)_{r,l}$            |          $(i,j-\frac{1}{2},k+\frac{1}{2})$         |
# |                  $\phi$                             |                      $(i,j,k)$                     |
# 
# We start by interpolating $\phi$ to  $\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)$.

# In[20]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n  /*****************************************\n   * COMPUTING RHS OF A_x, BOOKKEEPING NOTE:\n   * We want to compute\n   * \\partial_t A_x - [gauge terms] = \\psi^{6} (v^y B^z - v^z B^y).\n   * A_x is defined at (i,j+1/2,k+1/2).\n   * ==========================\n   * Where defined  | Variables\n   * (i,j-1/2,k-1/2)| {vyrr,vyrl,vylr,vyll,vzrr,vzrl,vzlr,vzll}\n   * (i,j+1/2,k-1/2)| {By_stagger_r,By_stagger_l} (see Table 1 in arXiv:1007.2848)\n   * (i,j-1/2,k+1/2)| {Bz_stagger_r,Bz_stagger_l} (see Table 1 in arXiv:1007.2848)\n   * (i,j,k)        | {phi}\n   * ==========================\n   ******************************************/\n  // Next compute phi at (i,j+1/2,k+1/2):\n#pragma omp parallel for\n  for(int k=1;k<cctk_lsh[2]-2;k++) for(int j=1;j<cctk_lsh[1]-2;j++) for(int i=0;i<cctk_lsh[0];i++) {\n        temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=\n          IPH(IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j-1,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+1,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+2,k-1)]),\n              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j-1,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+1,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+2,k  )]),\n              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j-1,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+1,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+2,k+1)]),\n              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j-1,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+1,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j+2,k+2)]));\n      }\n')


# Then we update the RHS of $\left[\partial_{t}A_{x}\right]_{\rm no\ gauge\ terms}$. Keep in mind that the function `A_i_rhs_no_gauge_terms()` takes care of determining $\left\{v^{i},B^{i},B^{i}_{\rm stagger}\right\}$ at $\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)$. We will look at it in more detail when we see the `A_i_rhs_no_gauge_terms.C` file from `IllinoisGRMHD`.

# In[21]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  int A_directionx=1;\n  A_i_rhs_no_gauge_terms(A_directionx,cctkGH,cctk_lsh,cctk_nghostzones,out_prims_r,out_prims_l,temporary,cmax_y,cmin_y,cmax_z,cmin_z, Ax_rhs);\n')


# <a id='rhs_ay_no_gauge_terms'></a>
# 
# ### Step 3.8.9: Evaluating $\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms}$ \[Back to [top](#toc)\]
# $$\label{rhs_ay_no_gauge_terms}$$
# 
# As a friendly reminder, we summarize the known gridpoint location of the needed gridfunctions here:
# 
# | Staggered and unstaggered variables                 | Gridpoint location at which the variable is known  |
# |-----------------------------------------------------|----------------------------------------------------|
# | $\left(v^{x}\right)_{r,l},\left(v^{z}\right)_{r,l}$ |          $(i-\frac{1}{2},j,k-\frac{1}{2})$         |
# | $\left(B^{x}_{\rm stagger}\right)_{r,l}$            |          $(i+\frac{1}{2},j,k-\frac{1}{2})$         |
# | $\left(B^{z}_{\rm stagger}\right)_{r,l}$            |          $(i-\frac{1}{2},j,k+\frac{1}{2})$         |
# |                  $\phi$                             |                      $(i,j,k)$                     |
# 
# We start by interpolating $\phi$ to  $\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)$.

# In[22]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  // We reprise flux_dirn=1 to finish up computations of Ai_rhs\'s!\n  flux_dirn=1;\n  // First compute ftilde, which is used for flattening left and right face values\n  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"\n  ftilde_gf_compute(cctkGH,cctk_lsh,flux_dirn, in_prims, ftilde_gf);\n\n  ww=0;\n  // NOTE! The order of variable reconstruction is important here,\n  //   as we don\'t want to overwrite {vxr,vxl,vyr,vyl}!\n  which_prims_to_reconstruct[ww]=VXR;       ww++;\n  which_prims_to_reconstruct[ww]=VZR;       ww++;\n  which_prims_to_reconstruct[ww]=VXL;       ww++;\n  which_prims_to_reconstruct[ww]=VZL;       ww++;\n  which_prims_to_reconstruct[ww]=BZ_STAGGER;ww++;\n  num_prims_to_reconstruct=ww;\n  // This function is housed in the file: "reconstruct_set_of_prims_PPM.C"\n  reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct,\n                               eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary);\n\n  /*****************************************\n   * COMPUTING RHS OF A_y, BOOKKEEPING NOTE:\n   * We want to compute\n   * \\partial_t A_y - [gauge terms] = \\psi^{6} (v^z B^x - v^x B^z).\n   * A_y is defined at (i+1/2,j,k+1/2).\n   * ==========================\n   * Where defined  | Variables\n   * (i-1/2,j,k-1/2)| {vxrr,vxrl,vxlr,vxll,vzrr,vzrl,vzlr,vzll}\n   * (i+1/2,j,k-1/2)| {Bx_stagger_r,Bx_stagger_l} (see Table 1 in arXiv:1007.2848)\n   * (i-1/2,j,k+1/2)| {Bz_stagger_r,Bz_stagger_l} (see Table 1 in arXiv:1007.2848)\n   * (i,j,k)        | {phi}\n   * ==========================\n   ******************************************/\n  // Next compute phi at (i+1/2,j,k+1/2):\n#pragma omp parallel for\n  for(int k=1;k<cctk_lsh[2]-2;k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=1;i<cctk_lsh[0]-2;i++) {\n        temporary[CCTK_GFINDEX3D(cctkGH,i,j,k)]=\n          IPH(IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j,k-1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j,k-1)]),\n              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j,k  )],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j,k  )]),\n              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j,k+1)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j,k+1)]),\n              IPH(phi_bssn[CCTK_GFINDEX3D(cctkGH,i-1,j,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i,j,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+1,j,k+2)],phi_bssn[CCTK_GFINDEX3D(cctkGH,i+2,j,k+2)]));\n      }\n')


# Then we update the RHS of $\left[\partial_{t}A_{y}\right]_{\rm no\ gauge\ terms}$. Keep in mind that the function `A_i_rhs_no_gauge_terms()` takes care of determining $\left\{v^{i},B^{i},B^{i}_{\rm stagger}\right\}$ at $\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)$. We will look at it in more detail when we see the `A_i_rhs_no_gauge_terms.C` file from `IllinoisGRMHD`.

# In[23]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  int A_directiony=2;\n  A_i_rhs_no_gauge_terms(A_directiony,cctkGH,cctk_lsh,cctk_nghostzones,out_prims_r,out_prims_l,temporary,cmax_z,cmin_z,cmax_x,cmin_x, Ay_rhs);\n')


# <a id='rhs_psi6phi_and_ai_gauge_terms'></a>
# 
# ### Step 3.8.10: Evaluating $\partial_{t}\left[\psi^{6}\Phi\right]$ and $\left[\partial_{t}A_{i}\right]_{\rm gauge\ terms}$ \[Back to [top](#toc)\]
# $$\label{rhs_psi6phi_and_ai_gauge_terms}$$
# 
# Finally, we compute
# 
# $$
# \partial_{t}\left[\sqrt{\gamma}\Phi\right] =
# \partial_{t}\left[\psi^{6}\Phi\right] = 
# -\partial_{j}\left(\alpha\sqrt{\gamma}A^{j} - \beta^{j}\left[\sqrt{\gamma}\Phi\right]\right)
# -\xi\alpha\left[\sqrt{\gamma}\Phi\right]\ ,
# $$
# 
# and
# 
# $$
# \left[\partial_{t}A_{i}\right]_{\rm gauge\ terms} = -\partial_{i}\left(\alpha\Phi - \beta^{j}A_{j}\right)\ .
# $$
# 
# Notice that we will need $A^{i}$ to compute $\partial_{t}\left[\psi^{6}\Phi\right]$, but we only have $A_{i}$, so we need to determine $\bar\gamma^{ij}$ (${\rm gtupij}$).

# In[24]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__driver_evaluate_MHD_rhs__C', '\n\n  // Next compute psi6phi_rhs, and add gauge terms to A_i_rhs terms!\n  //   Note that in the following function, we don\'t bother with reconstruction, instead interpolating.\n  // We need A^i, but only have A_i. So we add gtupij to the list of input variables.\n  CCTK_REAL *interp_vars[MAXNUMINTERP];\n  ww=0;\n  interp_vars[ww]=betax;   ww++;\n  interp_vars[ww]=betay;   ww++;\n  interp_vars[ww]=betaz;   ww++;\n  interp_vars[ww]=gtupxx;  ww++;\n  interp_vars[ww]=gtupxy;  ww++;\n  interp_vars[ww]=gtupxz;  ww++;\n  interp_vars[ww]=gtupyy;  ww++;\n  interp_vars[ww]=gtupyz;  ww++;\n  interp_vars[ww]=gtupzz;  ww++;\n  interp_vars[ww]=psi_bssn;ww++;\n  interp_vars[ww]=lapm1;   ww++;\n  interp_vars[ww]=Ax;      ww++;\n  interp_vars[ww]=Ay;      ww++;\n  interp_vars[ww]=Az;      ww++;\n  int max_num_interp_variables=ww;\n  if(max_num_interp_variables>MAXNUMINTERP) {CCTK_VError(VERR_DEF_PARAMS,"Error: Didn\'t allocate enough space for interp_vars[]."); }\n  // We are FINISHED with v{x,y,z}{r,l} and P{r,l} so we use these 8 gridfunctions\' worth of space as temp storage.\n  Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs(cctkGH,cctk_lsh,cctk_nghostzones,dX,interp_vars,psi6phi,\n                                                 vxr,vyr,vzr,vxl,vyl,vzl,Pr,Pl,\n                                                 psi6phi_rhs,Ax_rhs,Ay_rhs,Az_rhs);\n\n\n  return;\n  /*\n  // FUN DEBUGGING TOOL (trust me!):\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  //st_x_rhs[index]=0.0;\n  //st_y_rhs[index]=0.0;\n  //st_z_rhs[index]=0.0;\n  //rho_star_rhs[index]=0.0;\n  //tau_rhs[index]=0.0;\n\n  psi6phi_rhs[index] = 0.0;\n  Ax_rhs[index] = 0.0;\n  Ay_rhs[index] = 0.0;\n  Az_rhs[index] = 0.0;\n  }\n  */\n}\n\n// We add #include\'s here instead of compiling these separately to help ensure that functions are properly inlined.\n//    These files only include about 800 lines of code in total (~1200 lines in total), but it\'s arguably more\n//    convenient to edit a 600 line file than an 1800 line file, so I\'d prefer to leave this unconventional structure\n//    alone.\n#include "reconstruct_set_of_prims_PPM.C"\n#include "compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C"\n#include "add_fluxes_and_source_terms_to_hydro_rhss.C"\n#include "mhdflux.C"\n#include "A_i_rhs_no_gauge_terms.C"\n#include "Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C"\n\n')


# <a id='driver_evaluate_MHD_rhs__h'></a>
# 
# # Step 4: The `driver_evaluate_MHD_rhs.h` header file \[Back to [top](#toc)\]
# $$\label{driver_evaluate_MHD_rhs__h}$$
# 
# Now we generate the header file for the `driver_evaluate_MHD_rhs.C` file.

# In[25]:


get_ipython().run_cell_magic('writefile', '$outfile_path__driver_evaluate_MHD_rhs__h', '#ifndef DRIVER_EVALUATE_MHD_RHS_H_\n#define DRIVER_EVALUATE_MHD_RHS_H_\n\n/* PRIVATE FUNCTIONS, Called within driver_evaluate_MHD_rhs.C ONLY */\nstatic void ftilde_gf_compute(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,gf_and_gz_struct *input,CCTK_REAL *ftilde_gf);\nstatic void reconstruct_set_of_prims_PPM(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,const int num_prims_to_reconstruct,const int *which_prims_to_reconstruct,\n                                         eos_struct &eosi,gf_and_gz_struct *in_prims,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,\n                                         CCTK_REAL *ftilde_gf,CCTK_REAL *temporary);\n\nstatic void compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu\n(const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,CCTK_REAL *dX,CCTK_REAL **metric,gf_and_gz_struct *prims,\n CCTK_REAL **TUPmunu,eos_struct &eos, CCTK_REAL Gamma_th,\n CCTK_REAL *gupxy,CCTK_REAL *gupxz,CCTK_REAL *gupyz,\n CCTK_REAL *kxx,CCTK_REAL *kxy,CCTK_REAL *kxz,CCTK_REAL *kyy,CCTK_REAL *kyz,CCTK_REAL *kzz,\n CCTK_REAL *tau_rhs);\nstatic void A_i_rhs_no_gauge_terms(const int A_dirn, const cGH *cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,\n                                   CCTK_REAL *phi_interped,CCTK_REAL *cmax_1,CCTK_REAL *cmin_1,CCTK_REAL *cmax_2,CCTK_REAL *cmin_2, CCTK_REAL *A3_rhs);\n\nstatic 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 **interp_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\nstatic 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#include "harm_primitives_headers.h"\n\n#endif /* DRIVER_EVALUATE_MHD_RHS_H_ */\n\n')


# <a id='code_validation'></a>
# 
# # Step 5: Code validation \[Back to [top](#toc)\]
# $$\label{code_validation}$$
# 
# <a id='code_validation_driver_evaluate_MHD_rhs__c'></a>
# 
# ## Step 5.a: `driver_evaluate_MHD_rhs.C` \[Back to [top](#toc)\]
# $$\label{code_validation_driver_evaluate_MHD_rhs__c}$$
# 
# First, we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook.

# In[26]:


# Verify if the code generated by this tutorial module
# matches the original IllinoisGRMHD source code

# First download the original IllinoisGRMHD source code
import urllib
from os import path

original_IGM_file_url  = "https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/driver_evaluate_MHD_rhs.C"
original_IGM_file_name = "driver_evaluate_MHD_rhs-original.C"
original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)

# Then download the original IllinoisGRMHD source code
# We try it here in a couple of ways in an attempt to keep
# the code more portable
try:
    original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode("utf-8")
    # Write down the file the original IllinoisGRMHD source code
    with open(original_IGM_file_path,"w") as file:
        file.write(original_IGM_file_code)
except:
    try:
        original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode("utf-8")
        # Write down the file the original IllinoisGRMHD source code
        with open(original_IGM_file_path,"w") as file:
            file.write(original_IGM_file_code)
    except:
        # If all else fails, hope wget does the job
        get_ipython().system('wget -O $original_IGM_file_path $original_IGM_file_url')

# Perform validation
Validation__driver_evaluate_MHD_rhs__C  = get_ipython().getoutput('diff $original_IGM_file_path $outfile_path__driver_evaluate_MHD_rhs__C')

if Validation__driver_evaluate_MHD_rhs__C == []:
    # If the validation passes, we do not need to store the original IGM source code file
    get_ipython().system('rm $original_IGM_file_path')
    print("Validation test for driver_evaluate_MHD_rhs.C: PASSED!")
else:
    # If the validation fails, we keep the original IGM source code file
    print("Validation test for driver_evaluate_MHD_rhs.C: FAILED!")
    # We also print out the difference between the code generated
    # in this tutorial module and the original IGM source code
    print("Diff:")
    for diff_line in Validation__driver_evaluate_MHD_rhs__C:
        print(diff_line)


# <a id='code_validation_driver_evaluate_MHD_rhs__h'></a>
# 
# ## Step 5.b: `driver_evaluate_MHD_rhs.h` \[Back to [top](#toc)\]
# $$\label{code_validation_driver_evaluate_MHD_rhs__h}$$
# 
# First we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook.

# In[27]:


# Verify if the code generated by this tutorial module
# matches the original IllinoisGRMHD source code

# First download the original IllinoisGRMHD source code
import urllib
from os import path

original_IGM_file_url  = "https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/driver_evaluate_MHD_rhs.h"
original_IGM_file_name = "driver_evaluate_MHD_rhs-original.h"
original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)

# Then download the original IllinoisGRMHD source code
# We try it here in a couple of ways in an attempt to keep
# the code more portable
try:
    original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode("utf-8")
    # Write down the file the original IllinoisGRMHD source code
    with open(original_IGM_file_path,"w") as file:
        file.write(original_IGM_file_code)
except:
    try:
        original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode("utf-8")
        # Write down the file the original IllinoisGRMHD source code
        with open(original_IGM_file_path,"w") as file:
            file.write(original_IGM_file_code)
    except:
        # If all else fails, hope wget does the job
        get_ipython().system('wget -O $original_IGM_file_path $original_IGM_file_url')

# Perform validation
Validation__driver_evaluate_MHD_rhs__h  = get_ipython().getoutput('diff $original_IGM_file_path $outfile_path__driver_evaluate_MHD_rhs__h')

if Validation__driver_evaluate_MHD_rhs__h == []:
    # If the validation passes, we do not need to store the original IGM source code file
    get_ipython().system('rm $original_IGM_file_path')
    print("Validation test for driver_evaluate_MHD_rhs.h: PASSED!")
else:
    # If the validation fails, we keep the original IGM source code file
    print("Validation test for driver_evaluate_MHD_rhs.h: FAILED!")
    # We also print out the difference between the code generated
    # in this tutorial module and the original IGM source code
    print("Diff:")
    for diff_line in Validation__driver_evaluate_MHD_rhs__h:
        print(diff_line)


# <a id='latex_pdf_output'></a>
# 
# # Step 6: Output this notebook to $\LaTeX$-formatted PDF file \[Back to [top](#toc)\]
# $$\label{latex_pdf_output}$$
# 
# 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
# [Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.pdf](Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means).

# In[28]:


latex_nrpy_style_path = os.path.join(nrpy_dir_path,"latex_nrpy_style.tplx")
#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__driver_evaluate_MHD_rhs.tex
get_ipython().system('rm -f Tut*.out Tut*.aux Tut*.log')