{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"\n",
"\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:** Validated \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": [
"\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": [
"\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": [
"\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": [
"\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": [
"\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
}