IllinoisGRMHD/doc/Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Tutorial-IllinoisGRMHD: `convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C`\n",
    "\n",
    "## Authors: Leo Werneck & Zach Etienne\n",
    "\n",
    "<font color='red'>**This module is currently under development**</font>\n",
    "\n",
    "## In this tutorial module we explain the algorithm used to get the BSSN variables from the ADM variables. We will also explain how to impose the constraint that the conformal metric has unit determinant.\n",
    "\n",
    "### Required and recommended citations:\n",
    "\n",
    "* **(Required)** Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. ([arxiv:1501.07276](http://arxiv.org/abs/1501.07276)).\n",
    "* **(Required)** Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) ([astro-ph/0512420](https://arxiv.org/abs/astro-ph/0512420)).\n",
    "* **(Recommended)** Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 ([astro-ph/0210618](https://arxiv.org/abs/astro-ph/0210618))."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This module is organized as follows\n",
    "\n",
    "0. [Step 0](#src_dir): **Source directory creation**\n",
    "1. [Step 1](#introduction): **Introduction**\n",
    "1. [Step 2](#convert_adm_to_bssn__det_gammabar_eq_1): **`convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C`**\n",
    "    1. [Step 2.a](#physical_metric_and_its_determinant): *Loading the physical metric $\\gamma_{ij}$ and computing $\\gamma = \\det\\left(\\gamma_{ij}\\right)$*\n",
    "    1. [Step 2.b](#phi_and_psi): *Computing $\\phi$ and $\\psi$*\n",
    "    1. [Step 2.c](#enforce_det_gammabar_eq_1): *Enforcing the constraint $\\bar\\gamma = 1$*\n",
    "    1. [Step 2.d](#update_gfs): *Updating the gridfunctions*\n",
    "1. [Step 3](#code_validation): **Code validation**\n",
    "1. [Step 4](#latex_pdf_output): **Output this notebook to $\\LaTeX$-formatted PDF file**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='src_dir'></a>\n",
    "\n",
    "# Step 0: Source directory creation \\[Back to [top](#toc)\\]\n",
    "$$\\label{src_dir}$$\n",
    "\n",
    "We will now use the [cmdline_helper.py NRPy+ module](Tutorial-Tutorial-cmdline_helper.ipynb) to create the source directory within the `IllinoisGRMHD` NRPy+ directory, if it does not exist yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 0: Creation of the IllinoisGRMHD source directory\n",
    "# Step 0a: Add NRPy's directory to the path\n",
    "# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n",
    "import os,sys\n",
    "nrpy_dir_path = os.path.join(\"..\",\"..\")\n",
    "if nrpy_dir_path not in sys.path:\n",
    "    sys.path.append(nrpy_dir_path)\n",
    "\n",
    "# Step 0b: Load up cmdline_helper and create the directory\n",
    "import cmdline_helper as cmd\n",
    "IGM_src_dir_path = os.path.join(\"..\",\"src\")\n",
    "cmd.mkdir(IGM_src_dir_path)\n",
    "\n",
    "# Step 0c: Create the output file path\n",
    "outfile_path__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\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='introduction'></a>\n",
    "\n",
    "# Step 1: Introduction \\[Back to [top](#toc)\\]\n",
    "$$\\label{introduction}$$\n",
    "\n",
    "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$.\n",
    "\n",
    "**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."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convert_adm_to_bssn__det_gammabar_eq_1'></a>\n",
    "\n",
    "# Step 2: `convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C` \\[Back to [top](#toc)\\]\n",
    "$$\\label{convert_adm_to_bssn__det_gammabar_eq_1}$$\n",
    "\n",
    "<a id='physical_metric_and_its_determinant'></a>\n",
    "\n",
    "## Step 2.a: Loading the physical metric $\\gamma_{ij}$ and computing $\\gamma = \\det\\left(\\gamma_{ij}\\right)$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{physical_metric_and_its_determinant}$$\n",
    "\n",
    "We start by reading in the physical metric ${\\rm gij\\_physL} := \\gamma_{ij}$ and computing its determinant\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\begin{align}\n",
    "\\gamma = \\det\\left(\\gamma_{ij}\\right)\n",
    "&= \\gamma_{xx}\\gamma_{yy}\\gamma_{zz} + \\gamma_{xy}\\gamma_{yz}\\gamma_{xz} + \\gamma_{xz}\\gamma_{xy}\\gamma_{yz}\\\\\n",
    "&- \\gamma_{xz}\\gamma_{yy}\\gamma_{xz} - \\gamma_{xy}\\gamma_{xy}\\gamma_{zz} - \\gamma_{xx}\\gamma_{yz}\\gamma_{yz}\n",
    "\\end{align}\n",
    "}\\ .\n",
    "$$\n",
    "\n",
    "Notice that we have used the fact that $\\gamma_{ij}$ is symmetric above, i.e. $\\gamma_{ij}=\\gamma_{ji}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing ../src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C\n",
    "#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",
    "\n",
    "void 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);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='phi_and_psi'></a>\n",
    "\n",
    "## Step 2.b: Computing $\\phi$ and $\\psi$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{phi_and_psi}$$\n",
    "\n",
    "The BSSN formalism relates the physical metric, $\\gamma_{ij}$, to a conformal metric, $\\bar\\gamma{ij}$, via the relation\n",
    "\n",
    "$$\n",
    "\\boxed{\\bar\\gamma_{ij} = e^{-4\\phi}\\gamma_{ij}}\\ ,\n",
    "$$\n",
    "\n",
    "with the constraint that $\\bar\\gamma \\equiv \\det\\left(\\bar\\gamma_{ij}\\right) = 1$. This immediately implies that\n",
    "\n",
    "$$\n",
    "e^{-12\\phi} \\gamma = \\bar\\gamma = 1 \\implies -12\\phi = \\log\\left(\\frac{1}{\\gamma}\\right) = - \\log\\gamma\n",
    "\\implies \\boxed{\\phi = \\frac{1}{12}\\log\\gamma}\\ .\n",
    "$$\n",
    "\n",
    "Useful quantities to compute are\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\begin{align}\n",
    "\\psi      &\\equiv e^{\\phi}\\\\\n",
    "\\psi^{-4} &\\equiv e^{-4\\phi}\n",
    "\\end{align}\n",
    "}\\ .\n",
    "$$\n",
    "\n",
    "All the boxed quantities above are evaluated below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C\n",
    "\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;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='enforce_det_gammabar_eq_1'></a>\n",
    "\n",
    "## Step 2.c: Enforcing the constraint $\\bar\\gamma = 1$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{enforce_det_gammabar_eq_1}$$\n",
    "\n",
    "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$).\n",
    "\n",
    "$$\n",
    "\\boxed{\\bar{\\gamma}_{ij} \\to \\left(\\frac{1}{\\bar{\\gamma}}\\right)^{1/3} \\bar{\\gamma}_{ij}}\\ .\n",
    "$$\n",
    "\n",
    "We then start by evaluating\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\begin{align}\n",
    "\\bar\\gamma = \\det\\left(\\bar\\gamma_{ij}\\right)\n",
    "&= \\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}\\\\\n",
    "&- \\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}\n",
    "\\end{align}\n",
    "}\\ .\n",
    "$$\n",
    "\n",
    "In the code below we also define\n",
    "\n",
    "$$\n",
    "\\boxed{{\\rm gtijdet\\_Fm1o3} = \\left(\\frac{1}{\\bar{\\gamma}}\\right)^{1/3}}\\ .\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C\n",
    "\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\",\n",
    "i,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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='update_gfs'></a>\n",
    "\n",
    "## Step 2.d: Updating the gridfunctions \\[Back to [top](#toc)\\]\n",
    "$$\\label{update_gfs}$$\n",
    "\n",
    "Finally, we update the gridfunctions for $\\phi$, $\\psi$, $\\bar\\gamma_{ij}$, $\\gamma_{ij}$, and $\\bar\\gamma^{ij}$, respectively, in the code below.\n",
    "\n",
    "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)\n",
    "\n",
    "$$\n",
    "\\mathcal{M}^{-1} = \n",
    "\\begin{pmatrix}\n",
    "a & b & c\\\\\n",
    "b & d & e\\\\\n",
    "c & e & f\n",
    "\\end{pmatrix}^{-1}\n",
    "=\n",
    "\\frac{1}{\\det(\\mathcal{M})}\n",
    "\\begin{pmatrix}\n",
    "A & B & C\\\\\n",
    "B & D & E\\\\\n",
    "C & E & F\n",
    "\\end{pmatrix}\\ ,\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\boxed{\n",
    "\\begin{align}\n",
    "A &= \\ \\ \\left(df - ee\\right)\\\\\n",
    "B &= -\\left(bf - ce\\right)\\\\\n",
    "C &= \\ \\ \\left(be - cd\\right)\\\\\n",
    "D &= \\ \\ \\left(af - cc\\right)\\\\\n",
    "E &= -\\left(ae - bc\\right)\\\\\n",
    "F &= \\ \\ \\left(ad - bb\\right)\n",
    "\\end{align}\n",
    "}\\ .\n",
    "$$\n",
    "\n",
    "Notice that if we replace $\\mathcal{M}\\to\\bar\\gamma_{ij}$, then we have also $\\det(\\mathcal{M})\\to1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to ../src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C\n",
    "\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 3: Code validation \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "First we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Validation test for convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C: FAILED!\n",
      "Diff:\n",
      "1d0\n",
      "< #include \"cctk.h\"\n",
      "4a4,5\n",
      "> #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n",
      "> #include \"cctk.h\"\n",
      "5a7\n",
      "> #endif\n",
      "13c15,20\n",
      "<   DECLARE_CCTK_PARAMETERS; \n",
      "---\n",
      "> \n",
      "> #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n",
      ">  DECLARE_CCTK_PARAMETERS;\n",
      "> #endif\n",
      "> \n",
      "> \n",
      "31a39\n",
      "> \n",
      "42c50,51\n",
      "<         \n",
      "---\n",
      "> \n",
      "> \n",
      "46c55\n",
      "<         CCTK_REAL gtijdet = gtxxL * gtyyL * gtzzL + gtxyL * gtyzL * gtxzL + gtxzL * gtxyL * gtyzL - \n",
      "---\n",
      ">         CCTK_REAL gtijdet = gtxxL * gtyyL * gtzzL + gtxyL * gtyzL * gtxzL + gtxzL * gtxyL * gtyzL -\n",
      "52,55c61,64\n",
      "<         gtxyL = gtxyL * gtijdet_Fm1o3; \n",
      "<         gtxzL = gtxzL * gtijdet_Fm1o3; \n",
      "<         gtyyL = gtyyL * gtijdet_Fm1o3; \n",
      "<         gtyzL = gtyzL * gtijdet_Fm1o3; \n",
      "---\n",
      ">         gtxyL = gtxyL * gtijdet_Fm1o3;\n",
      ">         gtxzL = gtxzL * gtijdet_Fm1o3;\n",
      ">         gtyyL = gtyyL * gtijdet_Fm1o3;\n",
      ">         gtyzL = gtyzL * gtijdet_Fm1o3;\n",
      "58,60c67,76\n",
      "<         if(gtijdet<0.0) { CCTK_VWarn(CCTK_WARN_ALERT,__LINE__, __FILE__, CCTK_THORNSTRING,\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",
      "< \t\t\t\t     i,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",
      ">         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\",\n",
      "> i,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",
      "> \n",
      "92a109\n",
      "> \n"
     ]
    }
   ],
   "source": [
    "# Verify if the code generated by this tutorial module\n",
    "# matches the original IllinoisGRMHD source code\n",
    "\n",
    "# First download the original IllinoisGRMHD source code\n",
    "import urllib\n",
    "from os import path\n",
    "\n",
    "original_IGM_file_url  = \"https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C\"\n",
    "original_IGM_file_name = \"convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij-original.C\"\n",
    "original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)\n",
    "\n",
    "# Then download the original IllinoisGRMHD source code\n",
    "# We try it here in a couple of ways in an attempt to keep\n",
    "# the code more portable\n",
    "try:\n",
    "    original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
    "    # Write down the file the original IllinoisGRMHD source code\n",
    "    with open(original_IGM_file_path,\"w\") as file:\n",
    "        file.write(original_IGM_file_code)\n",
    "except:\n",
    "    try:\n",
    "        original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
    "        # Write down the file the original IllinoisGRMHD source code\n",
    "        with open(original_IGM_file_path,\"w\") as file:\n",
    "            file.write(original_IGM_file_code)\n",
    "    except:\n",
    "        # If all else fails, hope wget does the job\n",
    "        !wget -O $original_IGM_file_path $original_IGM_file_url\n",
    "\n",
    "# Perform validation\n",
    "Validation__ADM_to_BSSN__det_gammabar_eq_1__C  = !diff $original_IGM_file_path $outfile_path__ADM_to_BSSN__det_gammabar_eq_1__C\n",
    "\n",
    "if Validation__ADM_to_BSSN__det_gammabar_eq_1__C == []:\n",
    "    # If the validation passes, we do not need to store the original IGM source code file\n",
    "    !rm $original_IGM_file_path\n",
    "    print(\"Validation test for convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C: PASSED!\")\n",
    "else:\n",
    "    # If the validation fails, we keep the original IGM source code file\n",
    "    print(\"Validation test for convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C: FAILED!\")\n",
    "    # We also print out the difference between the code generated\n",
    "    # in this tutorial module and the original IGM source code\n",
    "    print(\"Diff:\")\n",
    "    for diff_line in Validation__ADM_to_BSSN__det_gammabar_eq_1__C:\n",
    "        print(diff_line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 4: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-IllinoisGRMHD__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)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "latex_nrpy_style_path = os.path.join(nrpy_dir_path,\"latex_nrpy_style.tplx\")\n",
    "#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.ipynb\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.tex\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.tex\n",
    "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.tex\n",
    "!rm -f Tut*.out Tut*.aux Tut*.log"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}