#!/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: `convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C`
# 
# ## Authors: Leo Werneck & Zach Etienne
# 
# <font color='red'>**This module is currently under development**</font>
# 
# ## This tutorial module explains the algorithm used to get the BSSN variables from the ADM variables. Further, it outlines how to impose the constraint that the conformal metric has a unit determinant.
# 
# ### 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](#convert_adm_to_bssn__det_gammabar_eq_1): **`convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C`**
#     1. [Step 2.a](#physical_metric_and_its_determinant): *Loading the physical metric $\gamma_{ij}$ and computing $\gamma = \det\left(\gamma_{ij}\right)$*
#     1. [Step 2.b](#phi_and_psi): *Computing $\phi$ and $\psi$*
#     1. [Step 2.c](#enforce_det_gammabar_eq_1): *Enforcing the constraint $\bar\gamma = 1$*
#     1. [Step 2.d](#update_gfs): *Updating the gridfunctions*
# 1. [Step 3](#code_validation): **Code validation**
# 1. [Step 4](#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__ADM_to_BSSN__det_gammabar_eq_1__C  = os.path.join(IGM_src_dir_path,"convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C")


# <a id='introduction'></a>
# 
# # Step 1: Introduction \[Back to [top](#toc)\]
# $$\label{introduction}$$
# 
# In this module, we will explain the procedure used within `IllinoisGRMHD` to compute the conformal metric, $\bar\gamma_{ij}$, and its inverse, $\bar\gamma^{ij}$, from the physical metric $\gamma_{ij}$. We will also describe how to compute $\phi$, the conformal factor, and $\psi\equiv e^{\phi}$. Finally, we also explain the procedure used to enforce the constraint $\bar\gamma = \det\left(\bar\gamma_{ij}\right) = 1$.
# 
# **A note on notation**: The notation used throughout the NRPy tutorial notebooks for the conformal metric is $\bar\gamma_{ij}$. However, in the literature, the notation $\tilde\gamma_{ij}$ is also used to represent the conformal metric. It is important to note that in `IllinoisGRMHD` we refer to the conformal metric using the *latter* notation, i.e. ${\rm gtij} := \tilde\gamma_{ij}$ and ${\rm gtupij} := \tilde\gamma^{ij}$. To keep the discussion consistent with the other notebooks, however, we will still present the discussion using the notation $\tilde\gamma_{ij} \to \bar\gamma_{ij}$. Bottom line, here $\tilde\gamma_{ij}$ and $\bar\gamma_{ij}$ represent exactly the same quantity, the conformal metric.

# <a id='convert_adm_to_bssn__det_gammabar_eq_1'></a>
# 
# # Step 2: `convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C` \[Back to [top](#toc)\]
# $$\label{convert_adm_to_bssn__det_gammabar_eq_1}$$
# 
# <a id='physical_metric_and_its_determinant'></a>
# 
# ## Step 2.a: Loading the physical metric $\gamma_{ij}$ and computing $\gamma = \det\left(\gamma_{ij}\right)$ \[Back to [top](#toc)\]
# $$\label{physical_metric_and_its_determinant}$$
# 
# We start by reading in the physical metric ${\rm gij\_physL} := \gamma_{ij}$ and computing its determinant
# 
# $$
# \boxed{
# \begin{align}
# \gamma = \det\left(\gamma_{ij}\right)
# &= \gamma_{xx}\gamma_{yy}\gamma_{zz} + \gamma_{xy}\gamma_{yz}\gamma_{xz} + \gamma_{xz}\gamma_{xy}\gamma_{yz}\\
# &- \gamma_{xz}\gamma_{yy}\gamma_{xz} - \gamma_{xy}\gamma_{xy}\gamma_{zz} - \gamma_{xx}\gamma_{yz}\gamma_{yz}
# \end{align}
# }\ .
# $$
# 
# Notice that we have used the fact that $\gamma_{ij}$ is symmetric above, i.e. $\gamma_{ij}=\gamma_{ji}$.

# In[2]:


get_ipython().run_cell_magic('writefile', '$outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C', '#include <cmath>\n#include <cstdio>\n#include <cstdlib>\n#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n#include "cctk.h"\n#include "cctk_Parameters.h"\n#endif\n\nvoid IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij\n(const cGH *cctkGH,const int *cctk_lsh,\n CCTK_REAL *gxx,CCTK_REAL *gxy,CCTK_REAL *gxz,CCTK_REAL *gyy,CCTK_REAL *gyz,CCTK_REAL *gzz,CCTK_REAL *alp,\n CCTK_REAL *gtxx,CCTK_REAL *gtxy,CCTK_REAL *gtxz,CCTK_REAL *gtyy,CCTK_REAL *gtyz,CCTK_REAL *gtzz,\n CCTK_REAL *gtupxx,CCTK_REAL *gtupxy,CCTK_REAL *gtupxz,CCTK_REAL *gtupyy,CCTK_REAL *gtupyz,CCTK_REAL *gtupzz,\n CCTK_REAL *phi,CCTK_REAL *psi,CCTK_REAL *lapm1) {\n\n#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n DECLARE_CCTK_PARAMETERS;\n#endif\n\n\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        CCTK_REAL gxx_physL=gxx[index];\n        CCTK_REAL gxy_physL=gxy[index];\n        CCTK_REAL gxz_physL=gxz[index];\n        CCTK_REAL gyy_physL=gyy[index];\n        CCTK_REAL gyz_physL=gyz[index];\n        CCTK_REAL gzz_physL=gzz[index];\n\n        /**********************************************************************\n         * Compute \\tilde{\\gamma_{ij}}, phi, and psi (BSSN) from g_{ij} (ADM) *\n         **********************************************************************/\n        CCTK_REAL gijdet = gxx_physL * gyy_physL * gzz_physL + gxy_physL * gyz_physL * gxz_physL + gxz_physL * gxy_physL * gyz_physL\n          - gxz_physL * gyy_physL * gxz_physL - gxy_physL * gxy_physL * gzz_physL - gxx_physL * gyz_physL * gyz_physL;\n\n        gijdet = fabs(gijdet);\n')


# <a id='phi_and_psi'></a>
# 
# ## Step 2.b: Computing $\phi$ and $\psi$ \[Back to [top](#toc)\]
# $$\label{phi_and_psi}$$
# 
# The BSSN formalism relates the physical metric, $\gamma_{ij}$, to a conformal metric, $\bar\gamma{ij}$, via the relation
# 
# $$
# \boxed{\bar\gamma_{ij} = e^{-4\phi}\gamma_{ij}}\ ,
# $$
# 
# with the constraint that $\bar\gamma \equiv \det\left(\bar\gamma_{ij}\right) = 1$. This immediately implies that
# 
# $$
# e^{-12\phi} \gamma = \bar\gamma = 1 \implies -12\phi = \log\left(\frac{1}{\gamma}\right) = - \log\gamma
# \implies \boxed{\phi = \frac{1}{12}\log\gamma}\ .
# $$
# 
# Useful quantities to compute are
# 
# $$
# \boxed{
# \begin{align}
# \psi      &\equiv e^{\phi}\\
# \psi^{-4} &\equiv e^{-4\phi}
# \end{align}
# }\ .
# $$
# 
# All the boxed quantities above are evaluated below.

# In[3]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C', '\n\n        CCTK_REAL phiL = (1.0/12.0) * log(gijdet);\n        CCTK_REAL psiL = exp(phiL);\n\n        CCTK_REAL Psim4 = 1.0/(psiL*psiL*psiL*psiL);\n        CCTK_REAL gtxxL = gxx_physL*Psim4;\n        CCTK_REAL gtxyL = gxy_physL*Psim4;\n        CCTK_REAL gtxzL = gxz_physL*Psim4;\n        CCTK_REAL gtyyL = gyy_physL*Psim4;\n        CCTK_REAL gtyzL = gyz_physL*Psim4;\n        CCTK_REAL gtzzL = gzz_physL*Psim4;\n')


# <a id='enforce_det_gammabar_eq_1'></a>
# 
# ## Step 2.c: Enforcing the constraint $\bar\gamma = 1$ \[Back to [top](#toc)\]
# $$\label{enforce_det_gammabar_eq_1}$$
# 
# The BSSN formalism evolves the determinant of the conformal metric through $\partial_{t}\bar\gamma = 0$. Since initially $\bar\gamma=1$, one would expect that $\bar\gamma=1$ at all times. However, numerical errors can cause the value of the determinant of the conformal metric to drift away from unity. To remedy this, we enforce the constraint $\bar\gamma=1$ by transforming the conformal metric according to the algebraic condition (cf. eq. 53 of [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/pdf/1712.07658.pdf) with $\hat\gamma = 1$).
# 
# $$
# \boxed{\bar{\gamma}_{ij} \to \left(\frac{1}{\bar{\gamma}}\right)^{1/3} \bar{\gamma}_{ij}}\ .
# $$
# 
# We then start by evaluating
# 
# $$
# \boxed{
# \begin{align}
# \bar\gamma = \det\left(\bar\gamma_{ij}\right)
# &= \bar\gamma_{xx}\bar\gamma_{yy}\bar\gamma_{zz} + \bar\gamma_{xy}\bar\gamma_{yz}\bar\gamma_{xz} + \bar\gamma_{xz}\bar\gamma_{xy}\bar\gamma_{yz}\\
# &- \bar\gamma_{xz}\bar\gamma_{yy}\bar\gamma_{xz} - \bar\gamma_{xy}\bar\gamma_{xy}\bar\gamma_{zz} - \bar\gamma_{xx}\bar\gamma_{yz}\bar\gamma_{yz}
# \end{align}
# }\ .
# $$
# 
# In the code below we also define
# 
# $$
# \boxed{{\rm gtijdet\_Fm1o3} = \left(\frac{1}{\bar{\gamma}}\right)^{1/3}}\ .
# $$

# In[4]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C', '\n\n        /*********************************\n         * Apply det gtij = 1 constraint *\n         *********************************/\n        CCTK_REAL gtijdet = gtxxL * gtyyL * gtzzL + gtxyL * gtyzL * gtxzL + gtxzL * gtxyL * gtyzL -\n          gtxzL * gtyyL * gtxzL - gtxyL * gtxyL * gtzzL - gtxxL * gtyzL * gtyzL;\n\n        CCTK_REAL gtijdet_Fm1o3 = fabs(1.0/cbrt(gtijdet));\n\n        gtxxL = gtxxL * gtijdet_Fm1o3;\n        gtxyL = gtxyL * gtijdet_Fm1o3;\n        gtxzL = gtxzL * gtijdet_Fm1o3;\n        gtyyL = gtyyL * gtijdet_Fm1o3;\n        gtyzL = gtyzL * gtijdet_Fm1o3;\n        gtzzL = gtzzL * gtijdet_Fm1o3;\n\n        if(gtijdet<0.0) {\n#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n          CCTK_VWarn(CCTK_WARN_ALERT,__LINE__, __FILE__, CCTK_THORNSTRING,\n#else\n          fprintf(stderr,\n#endif\n"WARNING: det[3-metric]<0.0 at point  %d %d %d | cctk_lsh: %d %d %d. Hopefully this is occurring in gz\'s! gtij_phys = %.2e %.2e %.2e %.2e %.2e %.2e gtij_new = %.2e %.2e %.2e %.2e %.2e %.2e | gijdet = %.2e | gtijdet = %.2e\\n",\ni,j,k,cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],gxx_physL,gxy_physL,gxz_physL,gyy_physL,gyz_physL,gzz_physL,gtxxL,gtxyL,gtxzL,gtyyL,gtyzL,gtzzL,-gijdet,gtijdet);\n}\n')


# <a id='update_gfs'></a>
# 
# ## Step 2.d: Updating the gridfunctions \[Back to [top](#toc)\]
# $$\label{update_gfs}$$
# 
# Finally, we update the gridfunctions for $\phi$, $\psi$, $\bar\gamma_{ij}$, $\gamma_{ij}$, and $\bar\gamma^{ij}$, respectively, in the code below.
# 
# If $\mathcal{M}$ is a $3\times3$ symmetric matrix and $\det(\mathcal{M})\neq0$, [recall that we have](https://en.wikipedia.org/wiki/Invertible_matrix#Inversion_of_3_%C3%97_3_matrices)
# 
# $$
# \mathcal{M}^{-1} = 
# \begin{pmatrix}
# a & b & c\\
# b & d & e\\
# c & e & f
# \end{pmatrix}^{-1}
# =
# \frac{1}{\det(\mathcal{M})}
# \begin{pmatrix}
# A & B & C\\
# B & D & E\\
# C & E & F
# \end{pmatrix}\ ,
# $$
# 
# where
# 
# $$
# \boxed{
# \begin{align}
# A &= \ \ \left(df - ee\right)\\
# B &= -\left(bf - ce\right)\\
# C &= \ \ \left(be - cd\right)\\
# D &= \ \ \left(af - cc\right)\\
# E &= -\left(ae - bc\right)\\
# F &= \ \ \left(ad - bb\right)
# \end{align}
# }\ .
# $$
# 
# Notice that if we replace $\mathcal{M}\to\bar\gamma_{ij}$, then we have also $\det(\mathcal{M})\to1$.

# In[5]:


get_ipython().run_cell_magic('writefile', '-a $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C', '\n\n        CCTK_REAL Psi4 = psiL*psiL*psiL*psiL;\n        /*****************************************\n         * Set all the needed BSSN gridfunctions *\n         *****************************************/\n        phi[index] = phiL;\n        psi[index] = psiL;\n\n        lapm1[index] = alp[index] - 1.0;\n\n        gtxx[index] = gtxxL;\n        gtxy[index] = gtxyL;\n        gtxz[index] = gtxzL;\n        gtyy[index] = gtyyL;\n        gtyz[index] = gtyzL;\n        gtzz[index] = gtzzL;\n\n        gxx[index] = gtxxL*Psi4;\n        gxy[index] = gtxyL*Psi4;\n        gxz[index] = gtxzL*Psi4;\n        gyy[index] = gtyyL*Psi4;\n        gyz[index] = gtyzL*Psi4;\n        gzz[index] = gtzzL*Psi4;\n\n        gtupxx[index] =   ( gtyyL * gtzzL - gtyzL * gtyzL );\n        gtupxy[index] = - ( gtxyL * gtzzL - gtyzL * gtxzL );\n        gtupxz[index] =   ( gtxyL * gtyzL - gtyyL * gtxzL );\n        gtupyy[index] =   ( gtxxL * gtzzL - gtxzL * gtxzL );\n        gtupyz[index] = - ( gtxxL * gtyzL - gtxyL * gtxzL );\n        gtupzz[index] =   ( gtxxL * gtyyL - gtxyL * gtxyL );\n      }\n}\n\n')


# <a id='code_validation'></a>
# 
# # Step 3: 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[6]:


# 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/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C"
original_IGM_file_name = "convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij-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__ADM_to_BSSN__det_gammabar_eq_1__C  = get_ipython().getoutput('diff $original_IGM_file_path $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C')

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


# <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-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.pdf](Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means).

# In[7]:


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