#!/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: harm_u2p_util.c
# 
# ## Authors: Leo Werneck & Zach Etienne
# 
# <font color='red'>**This module is currently under development**</font>
# 
# ## In this tutorial module we explain the utility functions needed by the conservative-to-primitive algorithm used by `HARM`. This module will likely be absorbed by another one once we finish documenting the code.
# 
# ### 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](#harm_utoprim_2d__c__eos_indep): **EOS independent routines**
#     1. [Step 2.a](#raise_g): *The `raise_g()` function*
#     1. [Step 2.b](#lower_g): *The `lower_g()` function*
#     1. [Step 2.c](#ncov_calc): *The `ncov_calc()` function*
# 1. [Step 3](#harm_utoprim_2d__c__eos_dep): **EOS dependent routines**
#     1. [Step 3.a](#pressure_rho0_u): *The `pressure_rho0_u()` function*
#     1. [Step 3.b](#pressure_rho0_w): *The `pressure_rho0_w()` function*
# 1. [Step n-1](#code_validation): **Code validation**
# 1. [Step n](#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__harm_u2p_util__c = os.path.join(IGM_src_dir_path,"harm_u2p_util.c")


# <a id='introduction'></a>
# 
# # Step 1: Introduction \[Back to [top](#toc)\]
# $$\label{introduction}$$

# <a id='harm_utoprim_2d__c__eos_indep'></a>
# 
# # Step 2: EOS independent routines` \[Back to [top](#toc)\]
# $$\label{harm_utoprim_2d__c__eos_indep}$$

# In[2]:


get_ipython().run_cell_magic('writefile', '$outfile_path__harm_u2p_util__c', '#ifndef __HARM_U2P_UTIL__C__\n#define __HARM_U2P_UTIL__C__\n/*\n  -------------------------------------------------------------------------------\n  Copyright 2005 Scott C. Noble, Charles F. Gammie,\n  Jonathan C. McKinney, and Luca Del Zanna\n\n\n  This file is part of PVS-GRMHD.\n\n  PVS-GRMHD is free software; you can redistribute it and/or modify\n  it under the terms of the GNU General Public License as published by\n  the Free Software Foundation; either version 2 of the License, or\n  (at your option) any later version.\n\n  PVS-GRMHD is distributed in the hope that it will be useful,\n  but WITHOUT ANY WARRANTY; without even the implied warranty of\n  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n  GNU General Public License for more details.\n\n  You should have received a copy of the GNU General Public License\n  along with PVS-GRMHD; if not, write to the Free Software\n  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\n\n  -------------------------------------------------------------------------------\n*/\n\n// Function prototypes for this file:\nstatic void raise_g(CCTK_REAL vcov[NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL vcon[NDIM]);\nstatic void lower_g(CCTK_REAL vcon[NDIM], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL vcov[NDIM]);\nstatic void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM]);\nstatic CCTK_REAL pressure_rho0_u(eos_struct eos, CCTK_REAL rho0, CCTK_REAL u);\nstatic CCTK_REAL pressure_rho0_w(eos_struct eos, CCTK_REAL rho0, CCTK_REAL w);\n\n// Inlined function used by this file\nstatic inline void compute_P_cold__eps_cold(eos_struct eos,CCTK_REAL rho_in, CCTK_REAL &P_cold,CCTK_REAL &eps_cold);\n')


# <a id='raise_g'></a>
# 
# ## Step 2.a: The `raise_g()` function \[Back to [top](#toc)\]
# $$\label{raise_g}$$
# 
# This is a simple function, used to *raise* the indices of a *covariant vector* using an inverse metric, $g^{\mu\nu}$. Usually the vector is 4-dimensional and $g^{\mu\nu}$ is the inverse physical ADM 4-metric, but the function can be used with arbitrary vectors and metrics. In other words, given a vector $v_{\mu}$ and an inverse metric $g^{\mu\nu}$ the function outputs
# 
# $$
# \boxed{v^{\mu} = g^{\mu\nu}v_{\nu}}\ .
# $$

# In[3]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__harm_u2p_util__c', '\n\n/**********************************************************************\n    raise_g():\n\n         -- calculates the contravariant form of a covariant tensor,\n            using the inverse of the metric;\n***********************************************************************/\nstatic void raise_g(CCTK_REAL vcov[NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL vcon[NDIM])\n{\n  int i,j;\n\n  for(i=0;i<NDIM;i++) {\n    vcon[i] = 0. ;\n    for(j=0;j<NDIM;j++)\n      vcon[i] += gcon[i][j]*vcov[j] ;\n  }\n\n  return ;\n}\n')


# <a id='lower_g'></a>
# 
# ## Step 2.b: The `lower_g()` function \[Back to [top](#toc)\]
# $$\label{lower_g}$$
# 
# This is a simple function, used to *lower* the indices of a *contravariant vector* using a metric, $g_{\mu\nu}$. Usually the vector is 4-dimensional and $g_{\mu\nu}$ is the physical ADM 4-metric, but the function can be used with arbitrary vectors and metrics. In other words, given a vector $v^{\mu}$ and a metric $g_{\mu\nu}$ the function outputs
# 
# $$
# \boxed{v_{\mu} = g_{\mu\nu}v^{\nu}}\ .
# $$

# In[4]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__harm_u2p_util__c', '\n\n/**********************************************************************\n     lower_g():\n\n          -- calculates the ocvariant form of a contravariant tensor\n             using the metric;\n***********************************************************************/\nstatic void lower_g(CCTK_REAL vcon[NDIM], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL vcov[NDIM])\n{\n  int i,j;\n\n  for(i=0;i<NDIM;i++) {\n    vcov[i] = 0. ;\n    for(j=0;j<NDIM;j++)\n      vcov[i] += gcov[i][j]*vcon[j] ;\n  }\n\n  return ;\n}\n')


# <a id='ncov_calc'></a>
# 
# ## Step 2.c: The `ncov_calc()` function \[Back to [top](#toc)\]
# $$\label{ncov_calc}$$
# 
# This simple function sets the covariant normal vector $n_{\mu} = \left(-\alpha,0,0,0\right)$.

# In[5]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__harm_u2p_util__c', '\n\n/**********************************************************************\n     ncov_calc():\n\n         -- calculates the covariant form of the normal vector to our\n            spacelike hypersurfaces ala the ADM formalism.\n\n         -- requires the inverse metric;\n***********************************************************************/\nstatic void ncov_calc(CCTK_REAL gcon[NDIM][NDIM],CCTK_REAL ncov[NDIM])\n{\n  CCTK_REAL lapse ;\n  int i;\n\n  lapse = sqrt(-1./gcon[0][0]) ;\n\n  ncov[0] = -lapse ;\n  for( i = 1; i < NDIM; i++) {\n    ncov[i] = 0. ;\n  }\n\n  return ;\n}\n')


# <a id='harm_utoprim_2d__c__eos_dep'></a>
# 
# # Step 3: EOS dependent routines
# $$\label{harm_utoprim_2d__c__eos_dep}$$

# <a id='pressure_rho0_u'></a>
# 
# ## Step 3.a: The `pressure_rho0_u()` function \[Back to [top](#toc)\]
# $$\label{pressure_rho0_u}$$
# 
# The $\Gamma$-law EOS implemented in `HARM` is
# 
# $$
# p_{\Gamma}\left(\rho_{b},u\right) = \left(\Gamma-1\right)u\ ,
# $$
# 
# where
# 
# $$
# u = \rho_{b}\epsilon\ .
# $$
# 
# Thus, the pre-PPEOS Patch version of this function was
# 
# ```c
# /**************************************************
#   The following functions assume a Gamma-law EOS:
# ***************************************************/
# 
# /* 
#    pressure as a function of rho0 and u 
#    this is used by primtoU and Utoprim_?D
# */
# static CCTK_REAL pressure_rho0_u(CCTK_REAL rho0, CCTK_REAL u)
# {
#   return((GAMMA  - 1.)*u) ;
# }
# ```
# 
# In the case of a hybrid EOS, however, we have $p_{\rm hybrid}=p_{\rm hybrid}\left(\rho_{b},\epsilon\right)$. To obtain $p_{\rm hybrid}\left(\rho_{b},u\right)$ we use:
# 
# $$
# p_{\rm hybrid} = P_{\rm cold}\left(\rho_{b}\right) + \left(\Gamma_{\rm th}-1\right)\rho_{b}\left[\epsilon - \epsilon_{\rm cold}\left(\rho_{b}\right)\right]
# \implies
# \boxed{
# p_{\rm hybrid}\left(\rho_{b},u\right) = P_{\rm cold}\left(\rho_{b}\right) + \left(\Gamma_{\rm th}-1\right)\left[u - \rho_{b}\epsilon_{\rm cold}\left(\rho_{b}\right)\right]\ ,
# }
# $$
# 
# where
# 
# $$
# \left\{
# \begin{align}
# P_{\rm cold}\left(\rho_{b}\right) = K_{\rm poly}\rho_{b}^{\Gamma_{\rm poly}}\ ,\\
# \epsilon_{\rm cold}\left(\rho_{b}\right) = C + \frac{P_{\rm cold}}{\rho_{b}\left(\Gamma_{\rm poly}-1\right)}\ ,
# \end{align}
# \right.
# $$
# 
# where $C$ is a precomputed integration constant which guarantees the continuity of $\epsilon_{\rm cold}\left(\rho_{b}\right)$ and the subscript "poly" indicates the local polytropic EOS quantity.

# In[6]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__harm_u2p_util__c', '\n/**************************************************\n  The following functions assume a Gamma-law EOS:\n***************************************************/\n\n/*\npressure as a function of rho0 and u\nthis is used by primtoU and Utoprim_?D\n*/\nstatic CCTK_REAL pressure_rho0_u(eos_struct eos, CCTK_REAL rho0, CCTK_REAL u)\n{\n\n  // Set up Gamma_th:\n#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n  DECLARE_CCTK_PARAMETERS;\n#endif\n\n  // Compute P_cold, eps_cold\n  CCTK_REAL P_cold, eps_cold;\n  compute_P_cold__eps_cold(eos,rho0, P_cold,eps_cold);\n\n  /* Compute the pressure as a function of rho_b (rho0) and\n   * u = rho_b * eps, using our hybrid EOS:\n   * .-------------------------------------------------------------.\n   * | p(rho_b,u) = P_cold + (Gamma_th - 1)*(u - rho_b * eps_cold) |\n   * .-------------------------------------------------------------.\n   */\n  return( P_cold + (Gamma_th - 1.0)*(u - rho0*eps_cold) );\n\n}\n')


# <a id='pressure_rho0_w'></a>
# 
# ## Step 3.2: The `pressure_rho0_w()` function \[Back to [top](#toc)\]
# $$\label{pressure_rho0_w}$$
# 
# The $\Gamma$-law EOS implemented in `HARM` is
# 
# $$
# p_{\Gamma} = \left(\Gamma-1\right)u\ .
# $$
# 
# We want now to obtain $p_{\Gamma}\left(\rho_{b},w\right)$, where
# 
# $$
# w = u + \rho_{b} + p\ .
# $$
# 
# Then
# 
# $$
# p_{\Gamma} = \left(\Gamma-1\right)\left(w - \rho_{b} - p_{\Gamma}\right)
# \implies
# \boxed{
# p_{\Gamma}\left(\rho_{b},w\right) = \frac{\left(\Gamma-1\right)}{\Gamma}\left(w-\rho_{b}\right)
# }\ .
# $$
# 
# Thus, the pre-PPEOS Patch version of this function was
# 
# ```c
# /* 
#    pressure as a function of rho0 and w = rho0 + u + p 
#    this is used by primtoU and Utoprim_1D
# */
# static CCTK_REAL pressure_rho0_w(CCTK_REAL rho0, CCTK_REAL w)
# {
#   return((GAMMA-1.)*(w - rho0)/GAMMA) ;
# }
# ```
# 
# For our hybrid EOS, we have
# 
# $$
# \begin{align}
# p_{\rm hybrid} &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\rho_{b}\left[\epsilon - \epsilon_{\rm cold}\right]\\
# &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\left[u - \rho_{b}\epsilon_{\rm cold}\right]\\
# &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\left[w-\rho_{b}-p_{\rm hybrid} - \rho_{b}\epsilon_{\rm cold}\right]\\
# &= P_{\rm cold} + \left(\Gamma_{\rm th}-1\right)\left[w - \rho_{b}\left(1+\epsilon_{\rm cold}\right)\right]- \left(\Gamma_{\rm th}-1\right)p_{\rm hybrid}\\
# \implies
# &\boxed{
# p_{\rm hybrid}\left(\rho_{b},w\right) = \frac{P_{\rm cold}}{\Gamma_{\rm th}} + \frac{\left(\Gamma_{\rm th}-1\right)}{\Gamma_{\rm th}}\left[w - \rho_{b}\left(1+\epsilon_{\rm cold}\right)\right]
# }
# \end{align}
# $$

# In[7]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__harm_u2p_util__c', '\n\n/*\n   pressure as a function of rho0 and w = rho0 + u + p\n   this is used by primtoU and Utoprim_1D\n*/\nstatic CCTK_REAL pressure_rho0_w(eos_struct eos, CCTK_REAL rho0, CCTK_REAL w)\n{\n\n  // Set up Gamma_th:\n#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n  DECLARE_CCTK_PARAMETERS;\n#endif\n\n  // Compute P_cold, eps_cold\n  CCTK_REAL P_cold, eps_cold;\n  compute_P_cold__eps_cold(eos,rho0, P_cold,eps_cold);\n\n  /* Compute the pressure as a function of rho_b (rho0) and\n   * w = u + rho_b + p, using our hybrid EOS:\n   *  ----------------------------------------------------------------------------\n   * | p(rho_b,w) = ( P_cold + (Gamma_th-1)*( w - rho_b*(1+eps_cold) ) )/Gamma_th |\n   *  ----------------------------------------------------------------------------\n   */\n  return( (P_cold + (Gamma_th-1.0)*( w - rho0*(1.0+eps_cold) ) )/Gamma_th );\n}\n#endif\n')


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

# In[8]:


# 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/harm_u2p_util.c"
original_IGM_file_name = "harm_u2p_util-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__harm_u2p_util__c  = get_ipython().getoutput('diff $original_IGM_file_path $outfile_path__harm_u2p_util__c')

if Validation__harm_u2p_util__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 harm_u2p_util.c: PASSED!")
else:
    # If the validation fails, we keep the original IGM source code file
    print("Validation test for harm_u2p_util.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__harm_u2p_util__c:
        print(diff_line)


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

# In[9]:


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__harm_u2p_util.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_u2p_util.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_u2p_util.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_u2p_util.tex
get_ipython().system('rm -f Tut*.out Tut*.aux Tut*.log')