{
 "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",
    "# [Brill-Lindquist](https://journals.aps.org/pr/abstract/10.1103/PhysRev.131.471) Initial data\n",
    "## Author: Zach Etienne\n",
    "###  Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "<font color='green'>**This module has been validated to exhibit convergence to zero of the Hamiltonian constraint violation at the expected order to the exact solution (see plot at bottom of [the exact initial data validation start-to-finish tutorial notebook](Tutorial-Start_to_Finish-BSSNCurvilinear-Exact_Initial_Data.ipynb); momentum constraint is zero), and all quantities have been validated against the [original SENR code](https://bitbucket.org/zach_etienne/nrpy).**</font>\n",
    "\n",
    "### NRPy+ Source Code for this module: [BrillLindquist.py](../edit/BSSN/BrillLindquist.py)\n",
    "\n",
    "## Introduction:\n",
    "This module sets up initial data for a merging black hole system in Cartesian coordinates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# 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\n",
    "1. [Step 2](#initialdata): Setting up Brill-Lindquist initial data\n",
    "1. [Step 3](#code_validation): Code Validation against BSSN/BrillLindquist 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",
    "\n",
    "Let's start by importing all the needed modules from Python/NRPy+:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:05:12.331340Z",
     "iopub.status.busy": "2021-03-07T17:05:12.330479Z",
     "iopub.status.idle": "2021-03-07T17:05:12.647089Z",
     "shell.execute_reply": "2021-03-07T17:05:12.647596Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Initialize core Python/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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initialdata'></a>\n",
    "\n",
    "# Step 2: Setting up Brill-Lindquist initial data \\[Back to [top](#toc)\\]\n",
    "$$\\label{initialdata}$$ \n",
    "\n",
    "Here we set up Brill-Lindquist initial data ([Brill & Lindquist, Phys. Rev. 131, 471, 1963](https://journals.aps.org/pr/abstract/10.1103/PhysRev.131.471); see also Eq. 1 of [Brandt & Brügmann, arXiv:gr-qc/9711015v1](https://arxiv.org/pdf/gr-qc/9711015v1.pdf)):\n",
    "\n",
    "The conformal factor $\\psi$ for Brill-Lindquist initial data is given by\n",
    "$$\\psi = e^{\\phi} = 1 + \\sum_{i=1}^N \\frac{m_{(i)}}{2 \\left|\\vec{r}_{(i)} - \\vec{r}\\right|};\\quad K_{ij}=0,$$\n",
    "\n",
    "where $\\psi$ is written in terms of the 3-metric $\\gamma_{ij}$ as\n",
    "\n",
    "$$\n",
    "\\gamma_{ij} = \\psi^4 \\delta_{ij}.\n",
    "$$\n",
    "\n",
    "The extrinsic curvature is zero:\n",
    "$$\n",
    "K_{ij} = 0\n",
    "$$\n",
    "\n",
    "These data consist of $N$ nonspinning black holes initially at rest. This module restricts to the case of two such black holes, positioned anywhere. Here, we implement $N=2$.\n",
    "\n",
    "**Inputs for $\\psi$**:\n",
    "* The position and (bare) mass of black hole 1: $\\left(x_{(1)},y_{(1)},z_{(1)}\\right)$ and $m_{(1)}$, respectively\n",
    "* The position and (bare) mass of black hole 2: $\\left(x_{(2)},y_{(2)},z_{(2)}\\right)$ and $m_{(2)}$, respectively\n",
    "\n",
    "**Additional variables needed for spacetime evolution**:\n",
    "* Desired coordinate system\n",
    "* Desired initial lapse $\\alpha$ and shift $\\beta^i$. We will choose our gauge conditions as $\\alpha=\\psi^{-2}$ and $\\beta^i=B^i=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:05:12.659508Z",
     "iopub.status.busy": "2021-03-07T17:05:12.658571Z",
     "iopub.status.idle": "2021-03-07T17:05:12.727941Z",
     "shell.execute_reply": "2021-03-07T17:05:12.727267Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 2: Setting up Brill-Lindquist initial data\n",
    "\n",
    "thismodule = \"Brill-Lindquist\"\n",
    "BH1_posn_x,BH1_posn_y,BH1_posn_z = par.Cparameters(\"REAL\", thismodule,\n",
    "                                                   [\"BH1_posn_x\",\"BH1_posn_y\",\"BH1_posn_z\"],\n",
    "                                                   [         0.0,         0.0,        +0.5])\n",
    "BH1_mass = par.Cparameters(\"REAL\", thismodule, [\"BH1_mass\"],1.0)\n",
    "BH2_posn_x,BH2_posn_y,BH2_posn_z = par.Cparameters(\"REAL\", thismodule,\n",
    "                                                   [\"BH2_posn_x\",\"BH2_posn_y\",\"BH2_posn_z\"],\n",
    "                                                   [         0.0,         0.0,        -0.5])\n",
    "BH2_mass = par.Cparameters(\"REAL\", thismodule, [\"BH2_mass\"],1.0)\n",
    "\n",
    "# Step 2.a: Set spatial dimension (must be 3 for BSSN)\n",
    "DIM = 3\n",
    "par.set_parval_from_str(\"grid::DIM\",DIM)\n",
    "\n",
    "Cartxyz = ixp.declarerank1(\"Cartxyz\")\n",
    "\n",
    "# Step 2.b: Set psi, the conformal factor:\n",
    "psi = sp.sympify(1)\n",
    "psi += BH1_mass / ( 2 * sp.sqrt((Cartxyz[0]-BH1_posn_x)**2 + (Cartxyz[1]-BH1_posn_y)**2 + (Cartxyz[2]-BH1_posn_z)**2) )\n",
    "psi += BH2_mass / ( 2 * sp.sqrt((Cartxyz[0]-BH2_posn_x)**2 + (Cartxyz[1]-BH2_posn_y)**2 + (Cartxyz[2]-BH2_posn_z)**2) )\n",
    "\n",
    "# Step 2.c: Set all needed ADM variables in Cartesian coordinates\n",
    "gammaDD = ixp.zerorank2()\n",
    "KDD     = ixp.zerorank2() # K_{ij} = 0 for these initial data\n",
    "for i in range(DIM):\n",
    "    gammaDD[i][i] = psi**4\n",
    "\n",
    "alpha = 1/psi**2\n",
    "betaU = ixp.zerorank1() # We generally choose \\beta^i = 0 for these initial data\n",
    "BU    = ixp.zerorank1() # We generally choose B^i = 0 for these initial data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 3: Code Validation against BSSN/BrillLindquist NRPy+ 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 Brill-Lindquist initial data between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [BSSN/BrillLindquist.py](../edit/BSSN/BrillLindquist.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:05:12.768371Z",
     "iopub.status.busy": "2021-03-07T17:05:12.767381Z",
     "iopub.status.idle": "2021-03-07T17:05:25.097270Z",
     "shell.execute_reply": "2021-03-07T17:05:25.097749Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between Brill-Lindquist tutorial and NRPy+ BSSN.BrillLindquist module.\n",
      "ALL TESTS PASS\n"
     ]
    }
   ],
   "source": [
    "# Step 3: Code Validation against BSSN.BrillLindquist NRPy+ module\n",
    "\n",
    "import BSSN.BrillLindquist as bl\n",
    "bl.BrillLindquist()\n",
    "\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 Brill-Lindquist tutorial and NRPy+ BSSN.BrillLindquist module.\")\n",
    "compare(alpha, bl.alpha, \"alpha\", \"bl.alpha\")\n",
    "for i in range(DIM):\n",
    "    compare(betaU[i], bl.betaU[i], \"betaU\"+str(i), \"bl.betaU\"+str(i))\n",
    "    compare(BU[i], bl.BU[i], \"BU\"+str(i), \"bl.BU\"+str(i))\n",
    "    for j in range(DIM):\n",
    "        compare(gammaDD[i][j], bl.gammaDD[i][j], \"gammaDD\"+str(i)+str(j), \"bl.gammaDD\"+str(i)+str(j))\n",
    "        compare(KDD[i][j], bl.KDD[i][j], \"KDD\"+str(i)+str(j), \"bl.KDD\"+str(i)+str(j))\n",
    "print(\"ALL TESTS PASS\")"
   ]
  },
  {
   "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-ADM_Initial_Data-Brill-Lindquist](Tutorial-ADM_Initial_Data-Brill-Lindquist.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:05:25.102274Z",
     "iopub.status.busy": "2021-03-07T17:05:25.101538Z",
     "iopub.status.idle": "2021-03-07T17:05:28.441702Z",
     "shell.execute_reply": "2021-03-07T17:05:28.442352Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ADM_Initial_Data-Brill-Lindquist.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-ADM_Initial_Data-Brill-Lindquist.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ADM_Initial_Data-Brill-Lindquist\")"
   ]
  }
 ],
 "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": 2
}