{ "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 }