#!/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>
# 
# # Start-to-Finish Example: Unit Testing `GiRaFFE_NRPy`: Conservative-to-Primitive and Primitive-to-Conservative Solvers
# 
# ## Author: Patrick Nelson
# 
# **Notebook Status:** <font color='green'><b>Validated</b></font>
# 
# **Validation Notes:** This module will validate the routines in [Tutorial-GiRaFFE_NRPy-C2P_P2C](Tutorial-GiRaFFE_NRPy-C2P_P2C.ipynb).
# 
# ### NRPy+ Source Code for this module: 
# * [GiRaFFE_NRPy/GiRaFFE_NRPy_C2P_P2C.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_C2P_P2C.py) [\[**tutorial**\]](Tutorial-GiRaFFE_NRPy-C2P_P2C.ipynb) Generates the conservative-to-primitive and primitive-to-conservative solvers.
# 
# ## Introduction:
# 
# This notebook validates the NRPyfied C2P and P2C solvers against the original `GiRaFFE` code. This will be done at a point with a random but realistic spacetime and a variety of Poynting fluxes and Valencia velocities to test edge cases.
# 
# When this notebook is run, the significant digits of agreement between the old `GiRaFFE` and new `GiRaFFE_NRPy` versions of the algorithm outputs of Poynting flux and Valencia three-velocity will be printed to the screen right after the code is run [here](#compile_run).

# <a id='toc'></a>
# 
# # Table of Contents
# $$\label{toc}$$
# 
# This notebook is organized as follows
# 
# 1. [Step 1](#setup): Set up core functions and parameters for unit testing the C2P and P2C algorithms
#     1. [Step 1.a](#c2p) Conservative-to-Primitive Solver
#     1. [Step 1.b](#p2c) Primitive-to-Conservative Solver
#     1. [Step 1.c](#download) Download original `GiRaFFE` files
#     1. [Step 1.d](#free_params) Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h`
#     1. [Step 1.e](#interface) Create dummy files for the CCTK version of the code
# 1. [Step 2](#mainc): `C2P_P2C_unit_test.c`: The Main C Code
#     1. [Step 2.a](#compile_run): Compile and run the code to validate the output
# 1. [Step 3](#drift_notes): Output this notebook to $\LaTeX$-formatted PDF file
# 1. [Step 4](#latex_pdf_output): Output this notebook to $\LaTeX$-formatted PDF file

# <a id='setup'></a>
# 
# # Step 1: Set up core functions and parameters for unit testing the C2P and P2C algorithms \[Back to [top](#toc)\]
# 
# $$\label{setup}$$
# 
# We'll start by appending the relevant paths to `sys.path` so that we can access sympy modules in other places. Then, we'll import NRPy+ core functionality and set up a directory in which to carry out our test. We will also declare the gridfunctions that are needed for this portion of the code.

# In[1]:


import shutil, os, sys           # Standard Python modules for multiplatform OS-level functions
# First, we'll add the parent directory to the list of directories Python will check for modules.
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)

from outputC import *            # NRPy+: Core C code output module
import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends
import finite_difference as fin  # NRPy+: Finite difference C code generation module
import NRPy_param_funcs as par   # NRPy+: Parameter interface
import grid as gri               # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface

out_dir = "Validation/"
cmd.mkdir(out_dir)

thismodule = "Start_to_Finish_UnitTest-GiRaFFE_NRPy-C2P_P2C"

# Register the gridfunctions we need for this function
StildeD = ixp.register_gridfunctions_for_single_rank1("EVOL","StildeD")
BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU")
ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU")
gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01")
alpha = gri.register_gridfunctions("AUXEVOL","alpha")
betaU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","betaU")


# <a id='c2p'></a>
# 
# ## Step 1.a: Conservative-to-Primitive Solver \[Back to [top](#toc)\]
# 
# $$\label{c2p}$$
# 
# First, we'll use NRPy+ to build the C function that will apply fixes to  $\tilde{S}_i$ and recompute the velocity to match, along with the current sheet prescription. Note that the NRPy+ version of this code also speed-limits the velocities. 

# In[2]:


import GiRaFFE_NRPy.GiRaFFE_NRPy_C2P_P2C as C2P_P2C
C2P_P2C.GiRaFFE_NRPy_C2P(StildeD,BU,gammaDD,betaU,alpha)

values_to_print = [
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD0"),rhs=C2P_P2C.outStildeD[0]),
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD1"),rhs=C2P_P2C.outStildeD[1]),
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD2"),rhs=C2P_P2C.outStildeD[2]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU0"),rhs=C2P_P2C.ValenciavU[0]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU1"),rhs=C2P_P2C.ValenciavU[1]),
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU2"),rhs=C2P_P2C.ValenciavU[2])
                  ]

desc = "Apply fixes to \tilde{S}_i and recompute the velocity to match with current sheet prescription."
name = "GiRaFFE_NRPy_cons_to_prims"
outCfunction(
    outfile  = os.path.join(out_dir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs,REAL *in_gfs",
    body     = fin.FD_outputC("returnstring",values_to_print,params="outCverbose=False").replace("IDX4","IDX4S"),
    loopopts ="AllPoints,Read_xxs")


# <a id='p2c'></a>
# 
# ## Step 1.b: Primitive-to-Conservative Solver \[Back to [top](#toc)\]
# 
# $$\label{p2c}$$
# 
# Now, we'll output the function to solve for $\tilde{S}_i$ from the Valencia 3-velocity after the current sheet prescription is applied. 

# In[3]:


# Declare this symbol:
sqrt4pi = par.Cparameters("REAL",thismodule,"sqrt4pi","sqrt(4.0*M_PI)")
TINYDOUBLE = par.Cparameters("REAL",thismodule,"TINYDOUBLE",1e-100)

C2P_P2C.GiRaFFE_NRPy_P2C(gammaDD,betaU,alpha,  ValenciavU,BU, sqrt4pi)

values_to_print = [
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD0"),rhs=sp.simplify(C2P_P2C.StildeD[0])),
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD1"),rhs=sp.simplify(C2P_P2C.StildeD[1])),
                   lhrh(lhs=gri.gfaccess("in_gfs","StildeD2"),rhs=sp.simplify(C2P_P2C.StildeD[2])),
                  ]

desc = "Recompute StildeD after current sheet fix to Valencia 3-velocity to ensure consistency between conservative & primitive variables."
name = "GiRaFFE_NRPy_prims_to_cons"
outCfunction(
    outfile  = os.path.join(out_dir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs,REAL *in_gfs",
    body     = fin.FD_outputC("returnstring",values_to_print,params="outCverbose=False").replace("IDX4","IDX4S"),
    loopopts ="AllPoints")


# <a id='download'></a>
# 
# ## Step 1.c: Download original `GiRaFFE` files \[Back to [top](#toc)\]
# 
# $$\label{download}$$
# 
# Here, we download the relevant portion of the original `GiRaFFE` code from Bitbucket. 

# In[4]:


# First download the original IllinoisGRMHD source code
import urllib
original_file_url  = ["https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/GiRaFFE/src/GiRaFFE_headers.h",
                      "https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/GiRaFFE/src/inlined_functions.C",
                      "https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFE/src/driver_conserv_to_prims_FFE.C",
                      "https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFE/src/compute_conservatives_FFE.C",
                      "https://bitbucket.org/zach_etienne/wvuthorns/raw/231af720ccf3f1af50f7cce4a86b410fc8ea2e51/GiRaFFE/src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C"
                     ]
original_file_name = ["GiRaFFE_headers.h",
                      "inlined_functions.C",
                      "driver_conserv_to_prims_FFE.C",
                      "compute_conservatives_FFE.C",
                      "convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C"
                     ]

for i in range(len(original_file_url)):
    original_file_path = os.path.join(out_dir,original_file_name[i])

    # 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_file_code = urllib.request.urlopen(original_file_url[i]).read().decode('utf-8')
    except:
        original_file_code = urllib.urlopen(original_file_url[i]).read().decode('utf-8')

    # Write down the file the original IllinoisGRMHD source code
    with open(original_file_path,"w") as file:
        file.write(original_file_code)


# <a id='free_params'></a>
# 
# ## Step 1.d: Output C codes needed for declaring and setting Cparameters; also set `free_parameters.h` \[Back to [top](#toc)\]
# 
# $$\label{free_params}$$
# 
# Based on declared NRPy+ Cparameters, first we generate `declare_Cparameters_struct.h`, `set_Cparameters_default.h`, and `set_Cparameters[-SIMD].h`.
# 
# Then we output `free_parameters.h`, which sets some basic grid parameters as well as the speed limit parameter we need for this function.

# In[5]:


# Step 3.d.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
# par.generate_Cparameters_Ccodes(os.path.join(out_dir))

# Step 3.d.ii: Set free_parameters.h
with open(os.path.join(out_dir,"free_parameters.h"),"w") as file:
    file.write("""
// Set free-parameter values.

const int NGHOSTS = 0;

// Set free-parameter values for the initial data.
// Override parameter defaults with values based on command line arguments and NGHOSTS.
const int Nx0x1x2 = 1;
params.Nxx0 = Nx0x1x2;
params.Nxx1 = Nx0x1x2;
params.Nxx2 = Nx0x1x2;
params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*NGHOSTS;
params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*NGHOSTS;
params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*NGHOSTS;
// Step 0d: Set up space and time coordinates
// Step 0d.i: Declare \Delta x^i=dxx{0,1,2} and invdxx{0,1,2}, as well as xxmin[3] and xxmax[3]:
const REAL xxmin[3] = {0.0,0.0,0.0};
const REAL xxmax[3] = {0.1,0.1,0.1};

params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx0);
params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx1);
params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx2);
params.invdx0 = 1.0 / params.dxx0;
params.invdx1 = 1.0 / params.dxx1;
params.invdx2 = 1.0 / params.dxx2;

params.GAMMA_SPEED_LIMIT = 2000.0;
\n""")

# Generates declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
par.generate_Cparameters_Ccodes(os.path.join(out_dir))


# <a id='interface'></a>
# 
# ## Step 1.e: Create dummy files for the CCTK version of the code \[Back to [top](#toc)\]
# 
# $$\label{interface}$$
# 
# The original `GiRaFFE` code depends on some functionalities of the CCTK. Since we only care about this one small function, we can get around this by creating some nearly-empty, non-functional files that can be included to satisfy the pre-processor without changing functionality. We will later replace what little functionality we need with some basic global variables and macros.

# In[6]:


#incldue "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "Symmetry.h"
with open(os.path.join(out_dir,"cctk.h"),"w") as file:
    file.write("""//""")

with open(os.path.join(out_dir,"cctk_Arguments.h"),"w") as file:
    file.write("""#define DECLARE_CCTK_ARGUMENTS //
#define CCTK_ARGUMENTS void
""")

with open(os.path.join(out_dir,"cctk_Parameters.h"),"w") as file:
    file.write("""#define DECLARE_CCTK_PARAMETERS //
""")

with open(os.path.join(out_dir,"Symmetry.h"),"w") as file:
    file.write("""//""")


# <a id='mainc'></a>
# 
# # Step 2: `C2P_P2C_unit_test.c`: The Main C Code \[Back to [top](#toc)\]
# 
# $$\label{mainc}$$
# 
# Now that we have our vector potential and analytic magnetic field to compare against, we will start writing our unit test. We'll also import common C functionality, define `REAL`, the number of ghost zones, and the faces, and set the standard macros for NRPy+ style memory access.

# In[7]:


get_ipython().run_cell_magic('writefile', '$out_dir/C2P_P2C_unit_test.C', '// These are common packages that we are likely to need.\n#include "stdio.h"\n#include "stdlib.h"\n#include "math.h"\n#include "string.h" // Needed for strncmp, etc.\n#include "stdint.h" // Needed for Windows GCC 6.x compatibility\n#include <time.h>   // Needed to set a random seed.\n\n#define REAL double\n#include "declare_Cparameters_struct.h"\n\n// Standard NRPy+ memory access:\n#define IDX4S(g,i,j,k) \\\n( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )\n\n// Memory access definitions for NRPy+\n#define GAMMADD00GF 0\n#define GAMMADD01GF 1\n#define GAMMADD02GF 2\n#define GAMMADD11GF 3\n#define GAMMADD12GF 4\n#define GAMMADD22GF 5\n#define BETAU0GF 6\n#define BETAU1GF 7\n#define BETAU2GF 8\n#define ALPHAGF 9\n#define BU0GF 10\n#define BU1GF 11\n#define BU2GF 12\n#define VALENCIAVU0GF 13\n#define VALENCIAVU1GF 14\n#define VALENCIAVU2GF 15\n#define NUM_AUXEVOL_GFS 16\n\n#define STILDED0GF 0\n#define STILDED1GF 1\n#define STILDED2GF 2\n#define NUM_EVOL_GFS 3\n\n// Include the functions that we want to test:\n#include "GiRaFFE_NRPy_cons_to_prims.h"\n#include "GiRaFFE_NRPy_prims_to_cons.h"\n\n// Define CCTK macros\n#define CCTK_REAL double\n#define CCTK_INT int\nstruct cGH{};\nconst cGH* cctkGH;\n\n\n// GiRaFFE parameters in ETK\nconst CCTK_REAL min_radius_inside_of_which_conserv_to_prims_FFE_and_FFE_evolution_is_DISABLED = -1;\nconst int current_sheet_null_v = 1;\n\n// More definitions to interface with ETK code:\nconst int cctk_lsh[3] = {1,1,1};\nCCTK_REAL gxx[1];\nCCTK_REAL gxy[1];\nCCTK_REAL gxz[1];\nCCTK_REAL gyy[1];\nCCTK_REAL gyz[1];\nCCTK_REAL gzz[1];\nCCTK_REAL alp[1];\nCCTK_REAL gtxx[1];\nCCTK_REAL gtxy[1];\nCCTK_REAL gtxz[1];\nCCTK_REAL gtyy[1];\nCCTK_REAL gtyz[1];\nCCTK_REAL gtzz[1];\nCCTK_REAL gtupxx[1];\nCCTK_REAL gtupxy[1];\nCCTK_REAL gtupxz[1];\nCCTK_REAL gtupyy[1];\nCCTK_REAL gtupyz[1];\nCCTK_REAL gtupzz[1];\nCCTK_REAL phi_bssn[1];\nCCTK_REAL psi_bssn[1];\nCCTK_REAL lapm1[1];\nCCTK_REAL betax[1];\nCCTK_REAL betay[1];\nCCTK_REAL betaz[1];\nCCTK_REAL mhd_st_x[1];\nCCTK_REAL mhd_st_y[1];\nCCTK_REAL mhd_st_z[1];\nCCTK_REAL vx[1];\nCCTK_REAL vy[1];\nCCTK_REAL vz[1];\nCCTK_REAL Bx[1];\nCCTK_REAL By[1];\nCCTK_REAL Bz[1];\nCCTK_REAL x[1];\nCCTK_REAL y[1];\nCCTK_REAL z[1];\nCCTK_REAL r[1];\n\n// Define dz in CCTK\nCCTK_REAL cactus_dz;\n#define CCTK_DELTA_SPACE(i) cactus_dz\n\nCCTK_REAL GAMMA_SPEED_LIMIT = 2000.0;\n\n// Dummy ETK function:\n#define CCTK_GFINDEX3D(cctkGH,i,j,k) 0\n#define CCTK_VInfo(a01,a02,a03,a04,a05,a06,a07,a08,a09,a10,a11,a12) //\n#define CCTK_VWarn(b01,b02,b03,b04,b05,b06,b07,b08,b09,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19,b20,b21,b22,b23,b24,b25) //\n\n#include "driver_conserv_to_prims_FFE.C"\n#include "compute_conservatives_FFE.C"\n#include "convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C"\n\nint main() {\n    paramstruct params;\n#include "set_Cparameters_default.h"\n\n    // Step 0c: Set free parameters, overwriting Cparameters defaults\n    //          by hand or with command-line input, as desired.\n#include "free_parameters.h"\n#include "set_Cparameters-nopointer.h"\n\n    // Now that we\'ve calculated dxx2,  we can define a cactus equivalent\n    cactus_dz = dxx2;\n    // We\'ll define all indices to be 0. No need to complicate memory access\n    const int i0 = 0;\n    const int i1 = 0;\n    const int i2 = 0;\n\n    // This is the array to which we\'ll write the NRPy+ variables.\n    REAL *auxevol_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS);\n    REAL *evol_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS);\n\n    // Step 0e: Set up cell-centered Cartesian coordinate grids\n    REAL *xx[3];\n    xx[0] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS0);\n    xx[1] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS1);\n    xx[2] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS2);\n    for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] = xxmin[0] + (j-NGHOSTS)*dxx0;\n    for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + (j-NGHOSTS)*dxx1;\n    for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + (j-NGHOSTS)*dxx2;\n    x[0] = xx[0][0];\n    y[0] = xx[1][0];\n    z[0] = xx[2][0];\n    r[0] = sqrt(xx[0][0]*xx[0][0] + xx[1][0]*xx[1][0] + xx[2][0]*xx[2][0]);\n\n    // Now, it\'s time to make the random numbers.\n    //const long int seed = time(NULL); // seed = 1570632212; is an example of a seed that produces\n                                      // bad agreement for high speeds\n    //srand(seed); // Set the seed\n    //printf("seed for random number generator = %ld; RECORD IF AGREEMENT IS BAD\\\\n\\\\n",seed);\n    // We take care to make sure the corresponding quantities have the SAME value.\n    auxevol_gfs[IDX4S(ALPHAGF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;\n    const double alphaL = auxevol_gfs[IDX4S(ALPHAGF, i0,i1,i2)];\n    alp[0] = alphaL;\n\n    auxevol_gfs[IDX4S(GAMMADD00GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;\n    auxevol_gfs[IDX4S(GAMMADD01GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;\n    auxevol_gfs[IDX4S(GAMMADD02GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;\n    auxevol_gfs[IDX4S(GAMMADD11GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;\n    auxevol_gfs[IDX4S(GAMMADD12GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;\n    auxevol_gfs[IDX4S(GAMMADD22GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*0.2-0.1;\n\n    // Generated by NRPy+:\n       const double gammaDD00 = auxevol_gfs[IDX4S(GAMMADD00GF, i0,i1,i2)];\n       const double gammaDD01 = auxevol_gfs[IDX4S(GAMMADD01GF, i0,i1,i2)];\n       const double gammaDD02 = auxevol_gfs[IDX4S(GAMMADD02GF, i0,i1,i2)];\n       const double gammaDD11 = auxevol_gfs[IDX4S(GAMMADD11GF, i0,i1,i2)];\n       const double gammaDD12 = auxevol_gfs[IDX4S(GAMMADD12GF, i0,i1,i2)];\n       const double gammaDD22 = auxevol_gfs[IDX4S(GAMMADD22GF, i0,i1,i2)];\n       /*\n        * NRPy+ Finite Difference Code Generation, Step 2 of 1: Evaluate SymPy expressions and write to main memory:\n        */\n       const double tmp0 = gammaDD11*gammaDD22;\n       const double tmp1 = pow(gammaDD12, 2);\n       const double tmp2 = gammaDD02*gammaDD12;\n       const double tmp3 = pow(gammaDD01, 2);\n       const double tmp4 = pow(gammaDD02, 2);\n       const double tmp5 = gammaDD00*tmp0 - gammaDD00*tmp1 + 2*gammaDD01*tmp2 - gammaDD11*tmp4 - gammaDD22*tmp3;\n       const double tmp6 = 1.0/tmp5;\n\n    // Set the ETK metric:\n    gxx[0] = gammaDD00;\n    gxy[0] = gammaDD01;\n    gxz[0] = gammaDD02;\n    gyy[0] = gammaDD11;\n    gyz[0] = gammaDD12;\n    gzz[0] = gammaDD22;\n\n    auxevol_gfs[IDX4S(BETAU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;\n    const double betaxL = auxevol_gfs[IDX4S(BETAU0GF, i0,i1,i2)];\n    betax[0] = betaxL;\n    auxevol_gfs[IDX4S(BETAU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.2-0.1;\n    const double betayL = auxevol_gfs[IDX4S(BETAU1GF, i0,i1,i2)];\n    betay[0] = betayL;\n    auxevol_gfs[IDX4S(BETAU2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*0.02-0.01;\n    const double betazL = auxevol_gfs[IDX4S(BETAU2GF, i0,i1,i2)];\n    betaz[0] = betazL;\n\n    /* Generate physically meaningful speeds */\n    auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;\n    vx[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)]-betaxL;\n    auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;\n    vy[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)]-betayL;\n    auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;\n    vz[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]-betazL;\n    /* Superluminal speeds for testing */\n    /*auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;\n    vx[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)]-betaxL;\n    auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;\n    vy[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)]-betayL;\n    auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)] = 1.0+(double)rand()/RAND_MAX*9.0;\n    vz[0] = alphaL*auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]-betazL;*/\n\n    auxevol_gfs[IDX4S(BU0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;\n    Bx[0] = auxevol_gfs[IDX4S(BU0GF, i0,i1,i2)];\n    auxevol_gfs[IDX4S(BU1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;\n    By[0] = auxevol_gfs[IDX4S(BU1GF, i0,i1,i2)];\n    auxevol_gfs[IDX4S(BU2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;\n    Bz[0] = auxevol_gfs[IDX4S(BU2GF, i0,i1,i2)];\n\n    evol_gfs[IDX4S(STILDED0GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;\n    mhd_st_x[0] = evol_gfs[IDX4S(STILDED0GF, i0,i1,i2)];\n    evol_gfs[IDX4S(STILDED1GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;\n    mhd_st_y[0] = evol_gfs[IDX4S(STILDED1GF, i0,i1,i2)];\n    evol_gfs[IDX4S(STILDED2GF, i0,i1,i2)] = (double)rand()/RAND_MAX*2.0-1.0;\n    mhd_st_z[0] = evol_gfs[IDX4S(STILDED2GF, i0,i1,i2)];\n\n    GiRaFFE_conserv_to_prims_FFE();\n    GiRaFFE_NRPy_cons_to_prims(&params,xx,auxevol_gfs,evol_gfs);\n    GiRaFFE_NRPy_prims_to_cons(&params,xx,auxevol_gfs,evol_gfs);\n\n    //printf("Checking the Poynting Fluxes:\\n");\n    printf("%.1f, %.1f, %.1f\\n",1.0-log10(2.0*fabs(evol_gfs[IDX4S(STILDED0GF, i0,i1,i2)]-mhd_st_x[0])/(fabs(evol_gfs[IDX4S(STILDED0GF, i0,i1,i2)])+fabs(mhd_st_x[0]))),\n                                     1.0-log10(2.0*fabs(evol_gfs[IDX4S(STILDED1GF, i0,i1,i2)]-mhd_st_y[0])/(fabs(evol_gfs[IDX4S(STILDED1GF, i0,i1,i2)])+fabs(mhd_st_y[0]))),\n                                     1.0-log10(2.0*fabs(evol_gfs[IDX4S(STILDED2GF, i0,i1,i2)]-mhd_st_z[0])/(fabs(evol_gfs[IDX4S(STILDED2GF, i0,i1,i2)])+fabs(mhd_st_z[0]))));\n    //printf("NRPy: %.15e,%.15e,%.15e\\n",evol_gfs[IDX4S(STILDED0GF, i0,i1,i2)],evol_gfs[IDX4S(STILDED1GF, i0,i1,i2)],evol_gfs[IDX4S(STILDED2GF, i0,i1,i2)]);\n    //printf("CCTK: %.15e,%.15e,%.15e\\n",mhd_st_x[0],mhd_st_y[0],mhd_st_z[0]);\n    //printf("\\n\\n");\n    //printf("Checking the Valencia Velocities:\\n");\n    printf("%.1f, %.1f, %.1f\\n",1.0-log10(2.0*fabs(auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)]-(vx[0]+betaxL)/alphaL)/(fabs(auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)])+fabs((vx[0]+betaxL)/alphaL))),\n                                     1.0-log10(2.0*fabs(auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)]-(vy[0]+betayL)/alphaL)/(fabs(auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)])+fabs((vy[0]+betayL)/alphaL))),\n                                     1.0-log10(2.0*fabs(auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]-(vz[0]+betazL)/alphaL)/(fabs(auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)])+fabs((vz[0]+betazL)/alphaL))));\n    //printf("NRPy: %.15e,%.15e,%.15e\\n",auxevol_gfs[IDX4S(VALENCIAVU0GF, i0,i1,i2)],auxevol_gfs[IDX4S(VALENCIAVU1GF, i0,i1,i2)],auxevol_gfs[IDX4S(VALENCIAVU2GF, i0,i1,i2)]);\n    //printf("CCTK: %.15e,%.15e,%.15e\\n",(vx[0]+betaxL)/alphaL,(vy[0]+betayL)/alphaL,(vz[0]+betazL)/alphaL);\n}\n')


# <a id='compile_run'></a>
# 
# ## Step 2.a: Compile and run the code to validate the output \[Back to [top](#toc)\]
# 
# $$\label{compile_run}$$
# 
# Finally, we can compile and run the code we have written. Once run, this code will output the level of agreement between the two codes and some information to help interpret those numbers.

# In[8]:


import time
results_file = "out.txt"

print("Now compiling, should take ~2 seconds...\n")
start = time.time()
# Can't use C here; must use C++
# This command should work on all systems:
get_ipython().system('g++ -Ofast -fopenmp -march=native -funroll-loops Validation/C2P_P2C_unit_test.C -o Validation/C2P_P2C_unit_test -lm')
end = time.time()
print("Finished in "+str(end-start)+" seconds.\n\n")

os.chdir(os.path.join("Validation"))
cmd.Execute(os.path.join("C2P_P2C_unit_test"),file_to_redirect_stdout=os.path.join(results_file))
os.chdir(os.path.join(".."))


# Below are the numbers we care about. These are the Significant Digits of Agreement between the corrected StildeD and ValenciavU as computed by NRPy+ and ETK. Each row represents a flux direction; each entry therein corresponds to a component of StildeD. Each pair of outputs should show at least 10 significant digits of agreement. 

# In[9]:


import numpy as np
with open(os.path.join(out_dir,results_file),"r") as file:
    for i in range(2):
        output = file.readline()
        substrings = output.split(", ")
        sda = [float(i) for i in substrings]
        enough_digits = np.array(sda)>14
        if not enough_digits.all():
            sys.exit(1)


# <a id='drift_notes'></a>
# 
# ## Step 3: Some notes on the drift velocity \[Back to [top](#toc)\]
# $$\label{drift_notes}$$
# 
# While we have chosen to use the Valencia three-velocity in this version of `GiRaFFE`, we have also frequently made use of a relationship expressing this in terms of the drift velocity, which was used in the original `GiRaFFE`. The usefulness of this relationship to drift velocity extends beyond merely translating the original code. As discussed in [Paschalidis, et al.](https://arxiv.org/pdf/1310.3274.pdf), Sec. III.A (just above Eq. 45, with a proof in Appendix A), there is a one-parameter family of velocity definitions that fulfill the GRFFE conditions. The drift velocity sets this parameter to 0, which minimizes the Lorentz factor and *guarantees* that the four-velocity and magnetic fields are orthogonal to each other. This simplifies the form of $b^\mu$ and quantities that depend on it. 
# 
# This must be taken into account in developing this unit test, because NRPy+'s GRFFE module defaults to using a definition of $b^\mu$ that does not assume that this criterion is met, while the original `GiRaFFE` code assumes this in its C2P and P2C solvers. So, if we do not guarantee that our test data fulfills this criterion, these two different routines will produce different results. We will now go through the derivation of the equation used by `GiRaFFE` from first principles to show where this extra term appears.
# 
# This is the equation used by `GiRaFFE`, pulled from Eqs. 47 and 85 of [this paper](https://arxiv.org/abs/1310.3274):
# $$\tilde{S}_i = \gamma_{ij} \frac{\bar{v}^j \sqrt{\gamma}B^2}{4 \pi},$$
# or $$\tilde{S}_i = \gamma_{ij} \frac{(v^j + \beta^j) \sqrt{\gamma}B^2}{4 \pi \alpha},$$
# where $\bar{v}^j$ is the Valencia 3-velocity and $v^j$ is the drift velocity.
# 
# In IllinoisGRMHD (IGM), the expression used is $\tilde{S}_i = \alpha \sqrt{\gamma} T^0_{{\rm EM}i},$
# where 
# \begin{align}
# T^{\mu\nu}_{\rm EM} &= b^2 u^{\mu} u^{\nu} + \frac{1}{2} b^2 g^{\mu\nu} - b^\mu b^\nu \\
# b^0 &= \frac{u_j B^j}{\sqrt{4\pi} \alpha} \\
# b^i &= \frac{B^i + (u_j B^j) u^i}{\sqrt{4\pi} \alpha u^0} \\
# u^i &= u^0 v^i.
# \end{align}
# 

# Now, we'll work to bridge the gap between the expression in terms of the stress-energy tensor and that used in `GiRaFFE`. Consider that $T^0_{{\rm EM}i} = g_{i \mu} T^{0 \mu}_{{\rm EM}}.$ We'll thus start with the index lowering operation.
# \begin{align}
# g_{\mu \xi} T^{\xi \nu}_{\rm EM} &= g_{\mu \xi} \left( b^2 u^{\xi} u^{\nu} + \frac{1}{2} b^2 g^{\xi\nu} - b^\xi b^\nu \right) \\
# &= b^2 u_{\mu} u^{\nu} + \frac{1}{2} b^2 \delta^\nu_\mu - b_\mu b^\nu.
# \end{align}
# 
# After doing so, it is still apparent that we will need $b^2 = g_{\mu\nu} b^\mu b^\nu$, where 
# $$
# g_{\mu\nu} = \begin{pmatrix} 
# -\alpha^2 + \beta^k \beta_k & \beta_i \\
# \beta_j & \gamma_{ij}
# \end{pmatrix}.
# $$
# Expanding out the implied sum, we see that
# \begin{align}
# b^2 &= g_{00} b^0 b^0 + g_{i0} b^i b^0 + g_{0j} b^0 b^j + g_{ij} b^i b^j \\
# &= \left(-\alpha^2 + \beta^k \beta_k\right) b^0 b^0 + \beta_i b^i b^0 + \beta_j b^0 b^j + \gamma_{ij} b^i b^j \\
# &= \left(-\alpha^2 + \beta^k \beta_k\right) b^0 b^0 + 2 \beta_i b^i b^0 + \gamma_{ij} b^i b^j \\
# \end{align}
# 
# Now, it will be useful to plug in the definition of $b^\mu$. Recall that
# \begin{align}
# b^0 &= \frac{u_j B^j}{\sqrt{4\pi} \alpha} \\
# b^i &= \frac{B^i + (u_j B^j) u^i}{\sqrt{4\pi} \alpha u^0} \\
# \end{align}
# 
# Then 
# \begin{align}
# b^2 &= \frac{1}{4 \pi \alpha^2} \left( -\alpha^2 (u_k B^k)^2 + \beta^k \beta_k (u_k B^k)^2 + 2 \beta_i (u_k B^k) (B^i + (u_k B^k) u^i)/u^0 + \gamma_{ij}(B^i + (u_k B^k) u^i)(B^j + (u_k B^k) u^j)/(u^0)^2 \right) \\
# \end{align}
# 

# This ultimately reduces to (from Paschalidis, et al.):
# \begin{align}
# b^2 &= \frac{B^2 + (u_\mu B^\mu)^2}{4 \pi \alpha^2 (u^0)^2} \\
# \end{align}
# 
# So, we insert this into the stress energy tensor:
# $$
# T^0_{{\rm EM}i} = \frac{B^2 + (u_\mu B^\mu)^2}{4 \pi \alpha^2 (u^0)^2} u_i u^0 - \gamma_{ij} b^j b^0
# $$
# Note that since $B^\mu$ is purely spatial, $B^0=0$, thus $u_\mu B^\mu = u_k B^k$
# 
# \begin{align}
# T^0_{{\rm EM}i} &= \frac{B^2 + (u_k B^k)^2}{4 \pi \alpha^2 u^0} u_i - \gamma_{ij} \frac{(B^j + (u_k B^k) u^j)(u_k B^k)}{4 \pi \alpha^2 u^0} \\
# &= \frac{1}{4 \pi \alpha^2 u^0} \left( B^2 u_i + (u_k B^k)^2 u_i - \gamma_{ij}(B^j + (u_k B^k) u^j)(u_k B^k)\right) \\
# &= \frac{1}{4 \pi \alpha^2 u^0} \left( B^2 u_i + (u_k B^k)^2 u_i - (B_i + (u_k B^k) u_i)(u_k B^k)\right) \\
# \end{align}
# 
# Now, we will substitute in $u_j/u_0 = \beta_j + \gamma_{lj} v^l$:
# \begin{align}
# &= \frac{1}{4 \pi \alpha^2} \left( B^2 (\beta_i + \gamma_{li} v^l) - B_i (u_k B^k)\right) \\
# &= \frac{1}{4 \pi \alpha^2} \left( B^2 (\gamma_{ij} \beta^j + \gamma_{ij} v^j) - B_i (u_k B^k)\right) \\
# \end{align}
# 
# Now, we multiply by $\alpha \sqrt{\gamma}$: 
# \begin{align}
# \tilde{S}_i = \alpha \sqrt{\gamma} T^0_{{\rm EM}i} &= \gamma_{ij} \frac{(v^j + \beta^j) \sqrt{\gamma}B^2}{4 \pi \alpha} - \frac{\sqrt{\gamma}B_i (u_k B^k)}{4 \pi \alpha} 
# \end{align}
# 
# We note that there is an extra term here beyond the one that `GiRaFFE` uses for its P2C solver. This term exists because we did not make the assumption that $(u_k B^k) = 0$, which follows from the definition of drift velocity ([see here](https://arxiv.org/pdf/1310.3274.pdf), Sec. III.A, just above Eq. 45). Making that assumption, we get the expression used in `GiRaFFE`:
# \begin{align}
# \tilde{S}_i = \alpha \sqrt{\gamma} T^0_{{\rm EM}i} &= \gamma_{ij} \frac{(v^j + \beta^j) \sqrt{\gamma}B^2}{4 \pi \alpha}
# \end{align}
# 
# Normally, this relation does not come into play because $u_i \propto v_i \propto \tilde{S}_i$ and $\tilde{S}_i$ is calculated initially as the cross product of the electric and magnetic fields, so is guaranteed to be orthogonal to $B^i$. Here, however, we did nothing to guarantee this beforehand, so the difference shows up. Note also that the entire derivation becomes far easier if you assume this sooner since $b^0 = 0$ and $b^i = B^i / \left(\sqrt{4 \pi} \alpha u^0\right)$
# 
# To handle this more smoothly moving forward, we have added an extra function to the GRFFE module that computes $b^\mu$ under the assumption that $u_j B^j = 0$ for use only in GRFFE codes (i.e., without GRHD). This is desirable for more than just this small case since the extra term might cause small errors to propagate.

# <a id='latex_pdf_output'></a>
# 
# # Step 4: 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-Start_to_Finish_UnitTest-GiRaFFE_NRPy-C2P_P2C.pdf](Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-C2P_P2C.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

# In[10]:


import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
# cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-C2P_P2C",location_of_template_file=os.path.join(".."))