Tutorial-BSSN_in_terms_of_ADM.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",
    "# BSSN Quantities in terms of ADM Quantities\n",
    "## Author: Zach Etienne\n",
    "\n",
    "## This notebook showcases the conversion of ADM to BSSN variables using NRPy+, preparing them for solving Einstein's equations with the BSSN formulation. It covers the rescaling of 3-metric, extrinsic curvature, and gauge quantities, with an outline of the rescaling processes from covariant BSSN formulation and the underlying references.\n",
    "\n",
    "**Notebook Status:** <font color='orange'><b> Self-Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). **Additional validation tests may have been performed, but are as yet, undocumented. (TODO)**\n",
    "\n",
    "### NRPy+ Source Code for this module: [BSSN_in_terms_of_ADM.py](../edit/BSSN/BSSN_in_terms_of_ADM.py)\n",
    "\n",
    "## Introduction:\n",
    "This module documents the conversion of ADM variables:\n",
    "\n",
    "$$\\left\\{\\gamma_{ij}, K_{ij}, \\alpha, \\beta^i\\right\\}$$\n",
    "\n",
    "into BSSN variables\n",
    "\n",
    "$$\\left\\{\\bar{\\gamma}_{i j},\\bar{A}_{i j},\\phi, K, \\bar{\\Lambda}^{i}, \\alpha, \\beta^i, B^i\\right\\},$$ \n",
    "\n",
    "in the desired curvilinear basis (given by `reference_metric::CoordSystem`). Then it rescales the resulting BSSNCurvilinear variables (as defined in [the covariant BSSN formulation tutorial](Tutorial-BSSN_formulation.ipynb)) into the form needed for solving Einstein's equations with the BSSN formulation:\n",
    "\n",
    "$$\\left\\{h_{i j},a_{i j},\\phi, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table of Contents\n",
    "$$\\label{toc}$$ \n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#initializenrpy): Initialize core Python/NRPy+ modules; set desired output BSSN Curvilinear coordinate system set to Spherical\n",
    "1. [Step 2](#adm2bssn): Perform the ADM-to-BSSN conversion for 3-metric, extrinsic curvature, and gauge quantities\n",
    "    1. [Step 2.a](#adm2bssn_gamma): Convert ADM $\\gamma_{ij}$ to BSSN $\\bar{\\gamma}_{ij}$; rescale to get $h_{ij}$\n",
    "    1. [Step 2.b](#admexcurv_convert): Convert the ADM extrinsic curvature $K_{ij}$ to BSSN $\\bar{A}_{ij}$ and $K$; rescale to get $a_{ij}$, $K$.\n",
    "    1. [Step 2.c](#lambda): Define $\\bar{\\Lambda}^i$\n",
    "    1. [Step 2.d](#conformal): Define the conformal factor variable `cf`\n",
    "1. [Step 3](#code_validation): Code Validation against `BSSN.BSSN_in_terms_of_ADM` NRPy+ module\n",
    "1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Initialize core Python/NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:27.103925Z",
     "iopub.status.busy": "2021-03-07T17:21:27.103154Z",
     "iopub.status.idle": "2021-03-07T17:21:27.427720Z",
     "shell.execute_reply": "2021-03-07T17:21:27.428217Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Import needed core NRPy+ modules\n",
    "import sympy as sp                 # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import NRPy_param_funcs as par     # NRPy+: Parameter interface\n",
    "import indexedexp as ixp           # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm     # NRPy+: Reference metric support\n",
    "import sys                         # Standard Python modules for multiplatform OS-level functions\n",
    "import BSSN.BSSN_quantities as Bq  # NRPy+: This module depends on the parameter EvolvedConformalFactor_cf,\n",
    "                                   #        which is defined in BSSN.BSSN_quantities\n",
    "\n",
    "# Step 1.a: Set DIM=3, as we're using a 3+1 decomposition of Einstein's equations\n",
    "DIM=3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='adm2bssn'></a>\n",
    "\n",
    "# Step 2: Perform the ADM-to-BSSN conversion for 3-metric, extrinsic curvature, and gauge quantities \\[Back to [top](#toc)\\]\n",
    "$$\\label{adm2bssn}$$\n",
    "\n",
    "Here we convert ADM quantities to their BSSN Curvilinear counterparts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='adm2bssn_gamma'></a>\n",
    "\n",
    "## Step 2.a: Convert ADM $\\gamma_{ij}$ to BSSN $\\bar{\\gamma}_{ij}$; rescale to get $h_{ij}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{adm2bssn_gamma}$$\n",
    "\n",
    "We have (Eqs. 2 and 3 of [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):\n",
    "$$\n",
    "\\bar{\\gamma}_{i j} = \\left(\\frac{\\bar{\\gamma}}{\\gamma}\\right)^{1/3} \\gamma_{ij},\n",
    "$$\n",
    "where we always make the choice $\\bar{\\gamma} = \\hat{\\gamma}$.\n",
    "\n",
    "After constructing $\\bar{\\gamma}_{ij}$, we rescale to get $h_{ij}$ according to the prescription described in the [the covariant BSSN formulation tutorial](Tutorial-BSSN_formulation.ipynb) (also [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):\n",
    "\n",
    "$$\n",
    "h_{ij} = (\\bar{\\gamma}_{ij} - \\hat{\\gamma}_{ij})/\\text{ReDD[i][j]}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:27.436089Z",
     "iopub.status.busy": "2021-03-07T17:21:27.435407Z",
     "iopub.status.idle": "2021-03-07T17:21:27.437595Z",
     "shell.execute_reply": "2021-03-07T17:21:27.438163Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2: All ADM quantities were input into this function in the Spherical or Cartesian\n",
    "#         basis, as functions of r,th,ph or x,y,z, respectively. In Steps 1 and 2 above,\n",
    "#         we converted them to the xx0,xx1,xx2 basis, and as functions of xx0,xx1,xx2.\n",
    "#         Here we convert ADM quantities to their BSSN Curvilinear counterparts:\n",
    "\n",
    "# Step 2.a: Convert ADM $\\gamma_{ij}$ to BSSN $\\bar{gamma}_{ij}$:\n",
    "#           We have (Eqs. 2 and 3 of [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):\n",
    "def gammabarDD_hDD(gammaDD):\n",
    "    global gammabarDD,hDD\n",
    "    if rfm.have_already_called_reference_metric_function == False:\n",
    "        print(\"BSSN.BSSN_in_terms_of_ADM.hDD_given_ADM(): Must call reference_metric() first!\")\n",
    "        sys.exit(1)\n",
    "    # \\bar{gamma}_{ij} = (\\frac{\\bar{gamma}}{gamma})^{1/3}*gamma_{ij}.\n",
    "    _gammaUU, gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)\n",
    "    gammabarDD = ixp.zerorank2()\n",
    "    hDD        = ixp.zerorank2()\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            gammabarDD[i][j] = (rfm.detgammahat/gammaDET)**(sp.Rational(1,3))*gammaDD[i][j]\n",
    "            hDD[i][j] = (gammabarDD[i][j] - rfm.ghatDD[i][j]) / rfm.ReDD[i][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='admexcurv_convert'></a>\n",
    "\n",
    "## Step 2.b: Convert the ADM extrinsic curvature $K_{ij}$ to BSSN quantities $\\bar{A}_{ij}$ and $K={\\rm tr}(K_{ij})$; rescale $\\bar{A}_{ij}$ to get $a_{ij}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{admexcurv_convert}$$\n",
    "\n",
    "Convert the ADM extrinsic curvature $K_{ij}$ to the trace-free extrinsic curvature $\\bar{A}_{ij}$, plus the trace of the extrinsic curvature $K$, where (Eq. 3 of [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf)):\n",
    "\\begin{align}\n",
    "K &= \\gamma^{ij} K_{ij} \\\\\n",
    "\\bar{A}_{ij} &= \\left(\\frac{\\bar{\\gamma}}{\\gamma}\\right)^{1/3} \\left(K_{ij} - \\frac{1}{3} \\gamma_{ij} K \\right)\n",
    "\\end{align}\n",
    "\n",
    "After constructing $\\bar{A}_{ij}$, we rescale to get $a_{ij}$ according to the prescription described in the [the covariant BSSN formulation tutorial](Tutorial-BSSN_formulation.ipynb) (also [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):\n",
    "\n",
    "$$\n",
    "a_{ij} = \\bar{A}_{ij}/\\text{ReDD[i][j]}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:27.447400Z",
     "iopub.status.busy": "2021-03-07T17:21:27.446610Z",
     "iopub.status.idle": "2021-03-07T17:21:27.448700Z",
     "shell.execute_reply": "2021-03-07T17:21:27.449265Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.b: Convert the extrinsic curvature K_{ij} to the trace-free extrinsic\n",
    "#           curvature \\bar{A}_{ij}, plus the trace of the extrinsic curvature K,\n",
    "#           where (Eq. 3 of [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf)):\n",
    "def trK_AbarDD_aDD(gammaDD,KDD):\n",
    "    global trK,AbarDD,aDD\n",
    "    if rfm.have_already_called_reference_metric_function == False:\n",
    "        print(\"BSSN.BSSN_in_terms_of_ADM.trK_AbarDD(): Must call reference_metric() first!\")\n",
    "        sys.exit(1)\n",
    "    # \\bar{gamma}_{ij} = (\\frac{\\bar{gamma}}{gamma})^{1/3}*gamma_{ij}.\n",
    "    gammaUU, gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)\n",
    "    # K = gamma^{ij} K_{ij}, and\n",
    "    # \\bar{A}_{ij} &= (\\frac{\\bar{gamma}}{gamma})^{1/3}*(K_{ij} - \\frac{1}{3}*gamma_{ij}*K)\n",
    "    trK = sp.sympify(0)\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            trK += gammaUU[i][j]*KDD[i][j]\n",
    "\n",
    "    AbarDD = ixp.zerorank2()\n",
    "    aDD    = ixp.zerorank2()\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            AbarDD[i][j] = (rfm.detgammahat/gammaDET)**(sp.Rational(1,3))*(KDD[i][j] - sp.Rational(1,3)*gammaDD[i][j]*trK)\n",
    "            aDD[i][j]    = AbarDD[i][j] / rfm.ReDD[i][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='lambda'></a>\n",
    "\n",
    "## Step 2.c: Assuming the ADM 3-metric $\\gamma_{ij}$ is given as an explicit function of `(xx0,xx1,xx2)`, convert to BSSN $\\bar{\\Lambda}^i$; rescale to compute $\\lambda^i$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{lambda}$$\n",
    "\n",
    "To define $\\bar{\\Lambda}^i$ we implement Eqs. 4 and 5 of [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf):\n",
    "$$\n",
    "\\bar{\\Lambda}^i = \\bar{\\gamma}^{jk}\\left(\\bar{\\Gamma}^i_{jk} - \\hat{\\Gamma}^i_{jk}\\right).\n",
    "$$\n",
    "\n",
    "The [reference_metric.py](../edit/reference_metric.py) module provides us with exact, closed-form expressions for $\\hat{\\Gamma}^i_{jk}$, so here we need only compute exact expressions for $\\bar{\\Gamma}^i_{jk}$, based on $\\gamma_{ij}$ given as an explicit function of `(xx0,xx1,xx2)`. This is particularly useful when setting up initial data.\n",
    "\n",
    "After constructing $\\bar{\\Lambda}^i$, we rescale to get $\\lambda^i$ according to the prescription described in the [the covariant BSSN formulation tutorial](Tutorial-BSSN_formulation.ipynb) (also [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):\n",
    "\n",
    "$$\n",
    "\\lambda^i = \\bar{\\Lambda}^i/\\text{ReU[i]}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:27.459613Z",
     "iopub.status.busy": "2021-03-07T17:21:27.458929Z",
     "iopub.status.idle": "2021-03-07T17:21:27.461197Z",
     "shell.execute_reply": "2021-03-07T17:21:27.461662Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.c: Define \\bar{Lambda}^i (Eqs. 4 and 5 of [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf)):\n",
    "def LambdabarU_lambdaU__exact_gammaDD(gammaDD):\n",
    "    global LambdabarU,lambdaU\n",
    "\n",
    "    # \\bar{Lambda}^i = \\bar{gamma}^{jk}(\\bar{Gamma}^i_{jk} - \\hat{Gamma}^i_{jk}).\n",
    "    gammabarDD_hDD(gammaDD)\n",
    "    gammabarUU, _gammabarDET = ixp.symm_matrix_inverter3x3(gammabarDD)\n",
    "\n",
    "    # First compute Christoffel symbols \\bar{Gamma}^i_{jk}, with respect to barred metric:\n",
    "    GammabarUDD = ixp.zerorank3()\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            for k in range(DIM):\n",
    "                for l in range(DIM):\n",
    "                    GammabarUDD[i][j][k] += sp.Rational(1,2)*gammabarUU[i][l]*( sp.diff(gammabarDD[l][j],rfm.xx[k]) +\n",
    "                                                                                sp.diff(gammabarDD[l][k],rfm.xx[j]) -\n",
    "                                                                                sp.diff(gammabarDD[j][k],rfm.xx[l]) )\n",
    "    # Next evaluate \\bar{Lambda}^i, based on GammabarUDD above and GammahatUDD\n",
    "    #       (from the reference metric):\n",
    "    LambdabarU = ixp.zerorank1()\n",
    "    for i in range(DIM):\n",
    "        for j in range(DIM):\n",
    "            for k in range(DIM):\n",
    "                LambdabarU[i] += gammabarUU[j][k] * (GammabarUDD[i][j][k] - rfm.GammahatUDD[i][j][k])\n",
    "    for i in range(DIM):\n",
    "        # We evaluate LambdabarU[i] here to ensure proper cancellations. If these cancellations\n",
    "        #   are not applied, certain expressions (e.g., lambdaU[0] in StaticTrumpet) will\n",
    "        #   cause SymPy's (v1.5+) CSE algorithm to hang\n",
    "        LambdabarU[i] = LambdabarU[i].doit()\n",
    "    lambdaU    = ixp.zerorank1()\n",
    "    for i in range(DIM):\n",
    "        lambdaU[i] = LambdabarU[i] / rfm.ReU[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='conformal'></a>\n",
    "\n",
    "## Step 2.d: Define the conformal factor variable `cf` \\[Back to [top](#toc)\\]\n",
    "$$\\label{conformal}$$\n",
    "\n",
    "We define the conformal factor variable `cf` based on the setting of the `\"BSSN_quantities::EvolvedConformalFactor_cf\"` parameter.\n",
    "\n",
    "For example if `\"BSSN_quantities::EvolvedConformalFactor_cf\"` is set to `\"phi\"`, we can use Eq. 3 of [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf), which in arbitrary coordinates is written:\n",
    "\n",
    "$$\n",
    "\\phi = \\frac{1}{12} \\log\\left(\\frac{\\gamma}{\\bar{\\gamma}}\\right).\n",
    "$$\n",
    "\n",
    "Alternatively if `\"BSSN_quantities::EvolvedConformalFactor_cf\"` is set to `\"chi\"`, then\n",
    "$$\n",
    "\\chi = e^{-4 \\phi} = \\exp\\left(-4 \\frac{1}{12} \\left(\\frac{\\gamma}{\\bar{\\gamma}}\\right)\\right) \n",
    "= \\exp\\left(-\\frac{1}{3} \\log\\left(\\frac{\\gamma}{\\bar{\\gamma}}\\right)\\right) = \\left(\\frac{\\gamma}{\\bar{\\gamma}}\\right)^{-1/3}.\n",
    "$$\n",
    "\n",
    "Finally if `\"BSSN_quantities::EvolvedConformalFactor_cf\"` is set to `\"W\"`, then\n",
    "$$\n",
    "W = e^{-2 \\phi} = \\exp\\left(-2 \\frac{1}{12} \\log\\left(\\frac{\\gamma}{\\bar{\\gamma}}\\right)\\right) = \n",
    "\\exp\\left(-\\frac{1}{6} \\log\\left(\\frac{\\gamma}{\\bar{\\gamma}}\\right)\\right) = \n",
    "\\left(\\frac{\\gamma}{\\bar{\\gamma}}\\right)^{-1/6}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:27.469865Z",
     "iopub.status.busy": "2021-03-07T17:21:27.468755Z",
     "iopub.status.idle": "2021-03-07T17:21:27.471869Z",
     "shell.execute_reply": "2021-03-07T17:21:27.472581Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.d: Set the conformal factor variable cf, which is set\n",
    "#           by the \"BSSN_quantities::EvolvedConformalFactor_cf\" parameter. For example if\n",
    "#           \"EvolvedConformalFactor_cf\" is set to \"phi\", we can use Eq. 3 of\n",
    "#           [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf),\n",
    "#           which in arbitrary coordinates is written:\n",
    "def cf_from_gammaDD(gammaDD):\n",
    "    global cf\n",
    "\n",
    "    # \\bar{Lambda}^i = \\bar{gamma}^{jk}(\\bar{Gamma}^i_{jk} - \\hat{Gamma}^i_{jk}).\n",
    "    gammabarDD_hDD(gammaDD)\n",
    "    _gammabarUU, gammabarDET = ixp.symm_matrix_inverter3x3(gammabarDD)\n",
    "    _gammaUU, gammaDET       = ixp.symm_matrix_inverter3x3(gammaDD)\n",
    "\n",
    "    cf = sp.sympify(0)\n",
    "\n",
    "    if par.parval_from_str(\"EvolvedConformalFactor_cf\") == \"phi\":\n",
    "        # phi = \\frac{1}{12} log(\\frac{gamma}{\\bar{gamma}}).\n",
    "        cf = sp.Rational(1,12)*sp.log(gammaDET/gammabarDET)\n",
    "    elif par.parval_from_str(\"EvolvedConformalFactor_cf\") == \"chi\":\n",
    "        # chi = exp(-4*phi) = exp(-4*\\frac{1}{12}*(\\frac{gamma}{\\bar{gamma}}))\n",
    "        #      = exp(-\\frac{1}{3}*log(\\frac{gamma}{\\bar{gamma}})) = (\\frac{gamma}{\\bar{gamma}})^{-1/3}.\n",
    "        #\n",
    "        cf = (gammaDET/gammabarDET)**(-sp.Rational(1,3))\n",
    "    elif par.parval_from_str(\"EvolvedConformalFactor_cf\") == \"W\":\n",
    "        # W = exp(-2*phi) = exp(-2*\\frac{1}{12}*log(\\frac{gamma}{\\bar{gamma}}))\n",
    "        #   = exp(-\\frac{1}{6}*log(\\frac{gamma}{\\bar{gamma}})) = (\\frac{gamma}{bar{gamma}})^{-1/6}.\n",
    "        cf = (gammaDET/gammabarDET)**(-sp.Rational(1,6))\n",
    "    else:\n",
    "        print(\"Error EvolvedConformalFactor_cf type = \\\"\"+par.parval_from_str(\"EvolvedConformalFactor_cf\")+\"\\\" unknown.\")\n",
    "        sys.exit(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='betvet'></a>\n",
    "\n",
    "## Step 2.e: Rescale $\\beta^i$ and $B^i$ to compute $\\mathcal{V}^i={\\rm vet}^i$ and $\\mathcal{B}^i={\\rm bet}^i$, respectively \\[Back to [top](#toc)\\]\n",
    "$$\\label{betvet}$$\n",
    "\n",
    "We rescale $\\beta^i$ and $B^i$ according to the prescription described in the [the covariant BSSN formulation tutorial](Tutorial-BSSN_formulation.ipynb) (also [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):\n",
    "\\begin{align}\n",
    "\\mathcal{V}^i &= \\beta^i/\\text{ReU[i]}\\\\\n",
    "\\mathcal{B}^i &= B^i/\\text{ReU[i]}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:27.480259Z",
     "iopub.status.busy": "2021-03-07T17:21:27.479612Z",
     "iopub.status.idle": "2021-03-07T17:21:27.481888Z",
     "shell.execute_reply": "2021-03-07T17:21:27.482389Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2.e: Rescale beta^i and B^i according to the prescription described in\n",
    "#         the [BSSN in curvilinear coordinates tutorial notebook](Tutorial-BSSNCurvilinear.ipynb)\n",
    "#         (also [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):\n",
    "#\n",
    "# \\mathcal{V}^i &= beta^i/(ReU[i])\n",
    "# \\mathcal{B}^i &= B^i/(ReU[i])\n",
    "def betU_vetU(betaU,BU):\n",
    "    global vetU,betU\n",
    "\n",
    "    if rfm.have_already_called_reference_metric_function == False:\n",
    "        print(\"BSSN.BSSN_in_terms_of_ADM.bet_vet(): Must call reference_metric() first!\")\n",
    "        sys.exit(1)\n",
    "    vetU    = ixp.zerorank1()\n",
    "    betU    = ixp.zerorank1()\n",
    "    for i in range(DIM):\n",
    "        vetU[i]    =      betaU[i] / rfm.ReU[i]\n",
    "        betU[i]    =         BU[i] / rfm.ReU[i]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 3: Code Validation against `BSSN.BSSN_in_terms_of_ADM` module \\[Back to [top](#toc)\\] \n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for [UIUC initial data](Tutorial-ADM_Initial_Data-UIUC_BlackHole.ipynb) between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [BSSN.BSSN_in_terms_of_ADM](../edit/BSSN/BSSN_in_terms_of_ADM.py) module.\n",
    "\n",
    "As no basis transformation is performed, we analyze these expressions in their native, Spherical coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:27.497813Z",
     "iopub.status.busy": "2021-03-07T17:21:27.497089Z",
     "iopub.status.idle": "2021-03-07T17:21:28.240964Z",
     "shell.execute_reply": "2021-03-07T17:21:28.241439Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.a: Set the desired *output* coordinate system to Spherical:\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Spherical\")\n",
    "rfm.reference_metric()\n",
    "\n",
    "# Step 3.b: Set up initial data; assume UIUC spinning black hole initial data\n",
    "import BSSN.UIUCBlackHole as uibh\n",
    "uibh.UIUCBlackHole()\n",
    "\n",
    "# Step 3.c: Call above functions to convert ADM to BSSN curvilinear\n",
    "gammabarDD_hDD(                   uibh.gammaDD)\n",
    "trK_AbarDD_aDD(                   uibh.gammaDD,uibh.KDD)\n",
    "LambdabarU_lambdaU__exact_gammaDD(uibh.gammaDD)\n",
    "cf_from_gammaDD(                  uibh.gammaDD)\n",
    "betU_vetU(                        uibh.betaU,uibh.BU)\n",
    "\n",
    "# Step 3.d: Now load the BSSN_in_terms_of_ADM module and perform the same conversion\n",
    "import BSSN.BSSN_in_terms_of_ADM as BitoA\n",
    "BitoA.gammabarDD_hDD(                   uibh.gammaDD)\n",
    "BitoA.trK_AbarDD_aDD(                   uibh.gammaDD,uibh.KDD)\n",
    "BitoA.LambdabarU_lambdaU__exact_gammaDD(uibh.gammaDD)\n",
    "BitoA.cf_from_gammaDD(                  uibh.gammaDD)\n",
    "BitoA.betU_vetU(                        uibh.betaU,uibh.BU)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:27.497813Z",
     "iopub.status.busy": "2021-03-07T17:21:27.497089Z",
     "iopub.status.idle": "2021-03-07T17:21:28.240964Z",
     "shell.execute_reply": "2021-03-07T17:21:28.241439Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between this tutorial notebook and BSSN.BSSN_in_terms_of_ADM NRPy+ module:\n",
      "ALL TESTS PASS\n"
     ]
    }
   ],
   "source": [
    "# Step 3.e: Perform the consistency check\n",
    "def compare(q1, q2, q1name, q2name):\n",
    "    if sp.simplify(q1 - q2) != 0:\n",
    "        print(\"Error: \"+q1name+\" - \"+q2name+\" = \"+str(sp.simplify(q1 - q2)))\n",
    "        sys.exit(1) \n",
    "\n",
    "print(\"Consistency check between this tutorial notebook and BSSN.BSSN_in_terms_of_ADM NRPy+ module:\")\n",
    "\n",
    "compare(cf, BitoA.cf, \"cf\", \"BitoA.cf\")\n",
    "compare(trK, BitoA.trK, \"trK\", \"BitoA.trK\")\n",
    "# alpha is the only variable that remains unchanged:\n",
    "# print(\"alpha - BitoA.alpha = \" + str(alpha - BitoA.alpha))\n",
    "for i in range(DIM):\n",
    "    compare(vetU[i], BitoA.vetU[i], \"vetU\"+str(i), \"BitoA.vetU\"+str(i))\n",
    "    compare(betU[i], BitoA.betU[i], \"betU\"+str(i), \"BitoA.betU\"+str(i))\n",
    "    compare(lambdaU[i], BitoA.lambdaU[i], \"lambdaU\"+str(i), \"BitoA.lambdaU\"+str(i))\n",
    "    for j in range(DIM):\n",
    "        compare(hDD[i][j], BitoA.hDD[i][j], \"hDD\"+str(i)+str(j), \"BitoA.hDD\"+str(i)+str(j))\n",
    "        compare(aDD[i][j], BitoA.aDD[i][j], \"aDD\"+str(i)+str(j), \"BitoA.aDD\"+str(i)+str(j))\n",
    "\n",
    "print(\"ALL TESTS PASS\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: 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 [Tutorial-BSSN_in_terms_of_ADM.pdf](Tutorial-BSSN_in_terms_of_ADM.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": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:28.246581Z",
     "iopub.status.busy": "2021-03-07T17:21:28.245547Z",
     "iopub.status.idle": "2021-03-07T17:21:31.738562Z",
     "shell.execute_reply": "2021-03-07T17:21:31.739462Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-BSSN_in_terms_of_ADM.tex, and compiled LaTeX file to PDF\n",
      "    file Tutorial-BSSN_in_terms_of_ADM.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-BSSN_in_terms_of_ADM\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}