Tutorial-BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "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",
    "# Enforce conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint (Eq. 53 of [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658))\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This notebook uses a method for enhancing numerical relativity simulation stability, based on [Brown](https://arxiv.org/abs/0902.3652)'s covariant BSSN formalism. It resolves numerical errors causing temporal divergence in the conformal 3-metric $\\bar{\\gamma}{ij}$, which disrupts the underlying PDEs' hyperbolicity. A correction mechanism aligns the determinant of $\\bar{\\gamma}{ij}$ with the reference metric $\\hat{\\gamma}_{ij}$ at each Runge-Kutta timestep, safeguarding PDE hyperbolicity. \n",
    "\n",
    "**Notebook Status:** <font color='green'><b> 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). In addition, its output has been \n",
    "\n",
    "### NRPy+ Source Code for this module: [BSSN/Enforce_Detgammahat_Constraint.py](../edit/BSSN/Enforce_Detgammahat_Constraint.py)\n",
    "\n",
    "## Introduction:\n",
    "[Brown](https://arxiv.org/abs/0902.3652)'s covariant Lagrangian formulation of BSSN, which we adopt, requires that $\\partial_t \\bar{\\gamma} = 0$, where $\\bar{\\gamma}=\\det \\bar{\\gamma}_{ij}$. Further, all initial data we choose satisfies $\\bar{\\gamma}=\\hat{\\gamma}$. \n",
    "\n",
    "However, numerical errors will cause $\\bar{\\gamma}$ to deviate from a constant in time. This actually disrupts the hyperbolicity of the PDEs, so to cure this, we adjust $\\bar{\\gamma}_{ij}$ at the end of each Runge-Kutta timestep, so that its determinant satisfies $\\bar{\\gamma}=\\hat{\\gamma}$ at all times. We adopt the following, rather standard prescription (Eq. 53 of [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658)):\n",
    "\n",
    "$$\n",
    "\\bar{\\gamma}_{ij} \\to \\left(\\frac{\\hat{\\gamma}}{\\bar{\\gamma}}\\right)^{1/3} \\bar{\\gamma}_{ij}.\n",
    "$$"
   ]
  },
  {
   "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 needed NRPy+ modules\n",
    "1. [Step 2](#enforcegammaconstraint): Enforce the $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint\n",
    "1. [Step 3](#code_validation): Code Validation against `BSSN.Enforce_Detgammahat_Constraint` 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 needed NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:11.199625Z",
     "iopub.status.busy": "2021-03-07T17:21:11.198828Z",
     "iopub.status.idle": "2021-03-07T17:21:12.881189Z",
     "shell.execute_reply": "2021-03-07T17:21:12.880566Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step P1: import all needed modules from NRPy+:\n",
    "from outputC import nrpyAbs,lhrh,outCfunction # NRPy+: Core C code output module\n",
    "import finite_difference as fin   # NRPy+: Finite difference C code generation module\n",
    "import NRPy_param_funcs as par    # NRPy+: parameter interface\n",
    "import grid as gri                # NRPy+: Functions having to do with numerical grids\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 cmdline_helper as cmd      # NRPy+: Multi-platform Python command-line interface\n",
    "import sympy as sp                # SymPy, Python's core symbolic algebra package\n",
    "import BSSN.BSSN_quantities as Bq # NRPy+: BSSN quantities\n",
    "import os,shutil,sys              # Standard Python modules for multiplatform OS-level functions\n",
    "\n",
    "# Set spatial dimension (must be 3 for BSSN)\n",
    "DIM = 3\n",
    "par.set_parval_from_str(\"grid::DIM\",DIM)\n",
    "\n",
    "# Then we set the coordinate system for the numerical grid\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"SinhSpherical\")\n",
    "rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='enforcegammaconstraint'></a>\n",
    "\n",
    "# Step 2: Enforce the $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint \\[Back to [top](#toc)\\]\n",
    "$$\\label{enforcegammaconstraint}$$\n",
    "\n",
    "Recall that we wish to make the replacement:\n",
    "$$\n",
    "\\bar{\\gamma}_{ij} \\to \\left(\\frac{\\hat{\\gamma}}{\\bar{\\gamma}}\\right)^{1/3} \\bar{\\gamma}_{ij}.\n",
    "$$\n",
    "Notice the expression on the right is guaranteed to have determinant equal to $\\hat{\\gamma}$.\n",
    "\n",
    "$\\bar{\\gamma}_{ij}$ is not a gridfunction, so we must rewrite the above in terms of $h_{ij}$:\n",
    "\\begin{align}\n",
    "\\left(\\frac{\\hat{\\gamma}}{\\bar{\\gamma}}\\right)^{1/3} \\bar{\\gamma}_{ij} &= \\bar{\\gamma}'_{ij}  \\\\\n",
    "&= \\hat{\\gamma}_{ij} + \\varepsilon'_{ij} \\\\\n",
    "&= \\hat{\\gamma}_{ij} + \\text{Re[i][j]} h'_{ij} \\\\\n",
    "\\implies h'_{ij} &= \\left[\\left(\\frac{\\hat{\\gamma}}{\\bar{\\gamma}}\\right)^{1/3} \\bar{\\gamma}_{ij} - \\hat{\\gamma}_{ij}\\right] / \\text{Re[i][j]} \\\\\n",
    "&= \\left(\\frac{\\hat{\\gamma}}{\\bar{\\gamma}}\\right)^{1/3} \\frac{\\bar{\\gamma}_{ij}}{\\text{Re[i][j]}} - \\delta_{ij}\\\\\n",
    "&= \\left(\\frac{\\hat{\\gamma}}{\\bar{\\gamma}}\\right)^{1/3} \\frac{\\hat{\\gamma}_{ij} + \\text{Re[i][j]} h_{ij}}{\\text{Re[i][j]}} - \\delta_{ij}\\\\\n",
    "&= \\left(\\frac{\\hat{\\gamma}}{\\bar{\\gamma}}\\right)^{1/3} \\left(\\delta_{ij} + h_{ij}\\right) - \\delta_{ij}\n",
    "\\end{align}\n",
    "\n",
    "Upon inspection, when expressing $\\hat{\\gamma}$ SymPy generates expressions like `(xx0)^{4/3} = pow(xx0, 4./3.)`, which can yield $\\text{NaN}$s when `xx0 < 0` (i.e., in the `xx0` ghost zones). To prevent this, we know that $\\hat{\\gamma}\\ge 0$ for all reasonable coordinate systems, so we make the replacement $\\hat{\\gamma}\\to |\\hat{\\gamma}|$ below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:12.932489Z",
     "iopub.status.busy": "2021-03-07T17:21:12.919642Z",
     "iopub.status.idle": "2021-03-07T17:21:12.934693Z",
     "shell.execute_reply": "2021-03-07T17:21:12.935133Z"
    }
   },
   "outputs": [],
   "source": [
    "# We will need the h_{ij} quantities defined within BSSN_RHSs\n",
    "#    below when we enforce the gammahat=gammabar constraint\n",
    "# Step 1: All barred quantities are defined in terms of BSSN rescaled gridfunctions,\n",
    "#         which we declare here in case they haven't yet been declared elsewhere.\n",
    "\n",
    "Bq.declare_BSSN_gridfunctions_if_not_declared_already()\n",
    "hDD = Bq.hDD\n",
    "Bq.BSSN_basic_tensors()\n",
    "gammabarDD = Bq.gammabarDD\n",
    "\n",
    "# First define the Kronecker delta:\n",
    "KroneckerDeltaDD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    KroneckerDeltaDD[i][i] = sp.sympify(1)\n",
    "\n",
    "# The detgammabar in BSSN_RHSs is set to detgammahat when BSSN_RHSs::detgbarOverdetghat_equals_one=True (default),\n",
    "#    so we manually compute it here:\n",
    "dummygammabarUU, detgammabar = ixp.symm_matrix_inverter3x3(gammabarDD)\n",
    "\n",
    "# Next apply the constraint enforcement equation above.\n",
    "hprimeDD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        # Using nrpyAbs here, as it directly translates to fabs() without additional SymPy processing.\n",
    "        #    This acts to simplify the final expression somewhat.\n",
    "        hprimeDD[i][j] = \\\n",
    "        (nrpyAbs(rfm.detgammahat)/detgammabar)**(sp.Rational(1,3)) * (KroneckerDeltaDD[i][j] + hDD[i][j]) \\\n",
    "        - KroneckerDeltaDD[i][j]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 3: Code Validation against `BSSN.Enforce_Detgammahat_Constraint` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the C code output between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [BSSN.Enforce_Detgammahat_Constraint](../edit/BSSN/Enforce_Detgammahat_Constraint.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:12.945845Z",
     "iopub.status.busy": "2021-03-07T17:21:12.945185Z",
     "iopub.status.idle": "2021-03-07T17:21:13.792757Z",
     "shell.execute_reply": "2021-03-07T17:21:13.793209Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function enforce_detgammahat_constraint() to file enforce_detgammahat_constraint/enforce_detgammahat_constraint.h-validation\n",
      "Output C function enforce_detgammahat_constraint() to file enforce_detgammahat_constraint/enforce_detgammahat_constraint.h\n",
      "Validation test PASSED on file: enforce_detgammahat_constraint/enforce_detgammahat_constraint.h\n"
     ]
    }
   ],
   "source": [
    "##########\n",
    "# Step 1: Generate enforce_detgammahat_constraint() using functions in this tutorial notebook:\n",
    "\n",
    "Ccodesdir = os.path.join(\"enforce_detgammahat_constraint\")\n",
    "# First remove C code output directory if it exists\n",
    "# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
    "shutil.rmtree(Ccodesdir, ignore_errors=True)\n",
    "# Then create a fresh directory\n",
    "cmd.mkdir(Ccodesdir)\n",
    "\n",
    "enforce_detg_constraint_vars = [lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD00\"),rhs=hprimeDD[0][0]),\n",
    "                                lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD01\"),rhs=hprimeDD[0][1]),\n",
    "                                lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD02\"),rhs=hprimeDD[0][2]),\n",
    "                                lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD11\"),rhs=hprimeDD[1][1]),\n",
    "                                lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD12\"),rhs=hprimeDD[1][2]),\n",
    "                                lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD22\"),rhs=hprimeDD[2][2]) ]\n",
    "\n",
    "enforce_gammadet_string = fin.FD_outputC(\"returnstring\",enforce_detg_constraint_vars,\n",
    "                                         params=\"outCverbose=False,preindent=1,includebraces=False\")\n",
    "\n",
    "desc = \"Enforce det(gammabar) = det(gammahat) constraint.\"\n",
    "name = \"enforce_detgammahat_constraint\"\n",
    "outCfunction(\n",
    "    outfile=os.path.join(Ccodesdir, name + \".h-validation\"), desc=desc, name=name,\n",
    "    params=\"const rfm_struct *restrict rfmstruct, const paramstruct *restrict params, REAL *restrict in_gfs\",\n",
    "    body=enforce_gammadet_string,\n",
    "    loopopts=\"AllPoints,enable_rfm_precompute\")\n",
    "\n",
    "##########\n",
    "# Step 2: Generate enforce_detgammahat_constraint() using functions in BSSN.Enforce_Detgammahat_Constraint\n",
    "\n",
    "gri.glb_gridfcs_list = []\n",
    "\n",
    "import BSSN.Enforce_Detgammahat_Constraint as EGC\n",
    "EGC.output_Enforce_Detgammahat_Constraint_Ccode(outdir=Ccodesdir,\n",
    "                                                exprs=EGC.Enforce_Detgammahat_Constraint_symb_expressions())\n",
    "\n",
    "import filecmp\n",
    "for file in [os.path.join(Ccodesdir,\"enforce_detgammahat_constraint.h\")]:\n",
    "    if filecmp.cmp(file,file+\"-validation\") == False:\n",
    "        print(\"VALIDATION TEST FAILED on file: \"+file+\".\")\n",
    "        sys.exit(1)\n",
    "    else:\n",
    "        print(\"Validation test PASSED on file: \"+file)\n",
    "##########"
   ]
  },
  {
   "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-BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.pdf](Tutorial-BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.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": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:21:13.797559Z",
     "iopub.status.busy": "2021-03-07T17:21:13.796954Z",
     "iopub.status.idle": "2021-03-07T17:21:17.095814Z",
     "shell.execute_reply": "2021-03-07T17:21:17.096639Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-\n",
      "    BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.tex, and\n",
      "    compiled LaTeX file to PDF file Tutorial-\n",
      "    BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.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_enforcing_determinant_gammabar_equals_gammahat_constraint\")"
   ]
  }
 ],
 "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
}