{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# [BSSN](http://www2.yukawa.kyoto-u.ac.jp/~yuichiro.sekiguchi/3+1.pdf) Hamiltonian and momentum constraint equations, in ***curvilinear*** coordinates, using a covariant reference metric approach: C code generation\n", "\n", "## Authors: Ian Ruchlin & Zach Etienne\n", "### Formatting improvements courtesy Brandon Clark\n", "\n", "## This module constructs the BSSN Hamiltonian and momentum constraint equations as symbolic (SymPy) expressions, in terms of the core BSSN quantities $\\left\\{h_{i j},a_{i j},\\phi, K, \\lambda^{i}, \\alpha, \\mathcal{V}^i, \\mathcal{B}^i\\right\\}$, as defined in [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658) (see also [Baumgarte, Montero, Cordero-Carrión, and Müller (2012)](https://arxiv.org/abs/1211.6632).\n", "\n", "### This module implements a generic curvilinear coordinate reference metric approach matching that of [Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658), which is an extension of the spherical coordinate reference metric approach of [Baumgarte, Montero, Cordero-Carrión, and Müller (2012)](https://arxiv.org/abs/1211.6632), which builds upon the covariant \"Lagrangian\" BSSN formalism of [Brown (2009)](https://arxiv.org/abs/0902.3652). *See also citations within each article.*\n", "\n", "**Notebook Status:** Validated \n", "\n", "**Validation Notes:** All expressions generated in this module have been validated against a trusted code where applicable (the original NRPy+/SENR code, which itself was validated against [Baumgarte's code](https://arxiv.org/abs/1211.6632)).\n", "\n", "### NRPy+ Source Code for this module: [BSSN/BSSN_constraints.py](../edit/BSSN/BSSN_constraints.py)\n", "\n", "\n", "[comment]: <> (Introduction: TODO)\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 Python/NRPy+ modules\n", "1. [Step 2](#hamiltonianconstraint): Construct the Hamiltonian constraint $\\mathcal{H}$.\n", "1. [Step 3](#momentumconstraint): Construct the momentum constraint $\\mathcal{M}^i$.\n", "1. [Step 4](#code_validation): Code Validation against `BSSN.BSSN_constraints` NRPy+ module\n", "1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 1: Initialize needed Python/NRPy+ modules \\[Back to [top](#toc)\\]\n", "$$\\label{initializenrpy}$$\n", "\n", "We start by loading the needed modules. Notably, this module depends on several quantities defined in the [BSSN/BSSN_quantities.py](../edit/BSSN/BSSN_quantities.py) Python code, documented in the NRPy+ [BSSN quantities](Tutorial-BSSN_quantities.ipynb). In [Step 2](#hamiltonianconstraint) we call functions within [BSSN.BSSN_quantities](../edit/BSSN/BSSN_quantities.py) to define quantities needed in this module." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:21:01.823104Z", "iopub.status.busy": "2021-03-07T17:21:01.822043Z", "iopub.status.idle": "2021-03-07T17:21:02.561250Z", "shell.execute_reply": "2021-03-07T17:21:02.560504Z" } }, "outputs": [], "source": [ "# Step 1: Initialize needed Python/NRPy+ modules\n", "import sympy as sp # SymPy, Python's core symbolic algebra package on 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 grid as gri # NRPy+: Functions having to do with numerical grids\n", "import reference_metric as rfm # NRPy+: Reference metric support\n", "import BSSN.BSSN_quantities as Bq\n", "\n", "# Step 1.a: Set spatial dimension (must be 3 for BSSN, as BSSN is\n", "# a 3+1-dimensional decomposition of the general\n", "# relativistic field equations)\n", "DIM = 3\n", "\n", "# Step 1.b: Given the chosen coordinate system, set up\n", "# corresponding reference metric and needed\n", "# reference metric quantities\n", "# The following function call sets up the reference metric\n", "# and related quantities, including rescaling matrices ReDD,\n", "# ReU, and hatted quantities.\n", "rfm.reference_metric()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: $\\mathcal{H}$, the Hamiltonian constraint \\[Back to [top](#toc)\\]\n", "$$\\label{hamiltonianconstraint}$$\n", "\n", "Next we define the Hamiltonian constraint. Eq. 13 of [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf) yields:\n", "$$\n", "\\mathcal{H} = {\\underbrace {\\textstyle \\frac{2}{3} K^2}_{\\rm Term\\ 1}} - \n", "{\\underbrace {\\textstyle \\bar{A}_{ij} \\bar{A}^{ij}}_{\\rm Term\\ 2}} + \n", "{\\underbrace {\\textstyle e^{-4\\phi} \\left(\\bar{R} - 8 \\bar{D}^i \\phi \\bar{D}_i \\phi - 8 \\bar{D}^2 \\phi\\right)}_{\\rm Term\\ 3}}\n", "$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:21:02.635838Z", "iopub.status.busy": "2021-03-07T17:21:02.600086Z", "iopub.status.idle": "2021-03-07T17:21:03.413371Z", "shell.execute_reply": "2021-03-07T17:21:03.413908Z" } }, "outputs": [], "source": [ "# Step 2: The Hamiltonian constraint.\n", "# First declare all needed variables\n", "Bq.declare_BSSN_gridfunctions_if_not_declared_already() # Sets trK\n", "Bq.BSSN_basic_tensors() # Sets AbarDD\n", "Bq.gammabar__inverse_and_derivs() # Sets gammabarUU\n", "Bq.AbarUU_AbarUD_trAbar_AbarDD_dD() # Sets AbarUU and AbarDD_dD\n", "Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU() # Sets RbarDD\n", "Bq.phi_and_derivs() # Sets phi_dBarD & phi_dBarDD\n", "\n", "# Term 1: 2/3 K^2\n", "H = sp.Rational(2,3)*Bq.trK**2\n", "\n", "# Term 2: -A_{ij} A^{ij}\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " H += -Bq.AbarDD[i][j]*Bq.AbarUU[i][j]\n", "\n", "# Term 3a: trace(Rbar)\n", "Rbartrace = sp.sympify(0)\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " Rbartrace += Bq.gammabarUU[i][j]*Bq.RbarDD[i][j]\n", "\n", "# Term 3b: -8 \\bar{\\gamma}^{ij} \\bar{D}_i \\phi \\bar{D}_j \\phi = -8*phi_dBar_times_phi_dBar\n", "# Term 3c: -8 \\bar{\\gamma}^{ij} \\bar{D}_i \\bar{D}_j \\phi = -8*phi_dBarDD_contraction\n", "phi_dBar_times_phi_dBar = sp.sympify(0) # Term 3b\n", "phi_dBarDD_contraction = sp.sympify(0) # Term 3c\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " phi_dBar_times_phi_dBar += Bq.gammabarUU[i][j]*Bq.phi_dBarD[i]*Bq.phi_dBarD[j]\n", " phi_dBarDD_contraction += Bq.gammabarUU[i][j]*Bq.phi_dBarDD[i][j]\n", "\n", "# Add Term 3:\n", "H += Bq.exp_m4phi*(Rbartrace - 8*(phi_dBar_times_phi_dBar + phi_dBarDD_contraction))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 3: $\\mathcal{M}^i$, the momentum constraint \\[Back to [top](#toc)\\]\n", "$$\\label{momentumconstraint}$$\n", "\n", "***Courtesy Ian Ruchlin***\n", "\n", "The following definition of the momentum constraint is a simplification of Eq. 47 or [Ruchlin, Etienne, & Baumgarte (2018)](https://arxiv.org/pdf/1712.07658.pdf), which itself was a corrected version of the momentum constraint presented in Eq. 14 of [Baumgarte *et al*](https://arxiv.org/pdf/1211.6632.pdf).\n", "\n", "Start with the physical momentum constraint\n", "$$\n", "\\mathcal{M}^{i} \\equiv D_{j} \\left ( K^{i j} - \\gamma^{i j} K \\right ) = 0 \\; .\n", "$$\n", "Expanding and using metric compatibility with the physical covariant derivative $D_{i}$ yields\n", "$$\n", "\\mathcal{M}^{i} = D_{j} K^{i j} - \\gamma^{i j} \\partial_{j} K \\; .\n", "$$\n", "The physical extrinsic curvature $K_{i j}$ is related to the trace-free extrinsic curvature $A_{i j}$ by\n", "$$\n", "K_{i j} = A_{i j} + \\frac{1}{3} \\gamma_{i j} K \\; .\n", "$$\n", "Thus,\n", "$$\n", "\\mathcal{M}^{i} = D_{j} A^{i j} - \\frac{2}{3} \\gamma^{i j} \\partial_{j} K \\; .\n", "$$\n", "The physical metric $\\gamma_{i j}$ is related to the conformal metric $\\bar{\\gamma}_{i j}$ by the conformal rescaling\n", "$$\n", "\\gamma_{i j} = e^{4 \\phi} \\bar{\\gamma}_{i j} \\; ,\n", "$$\n", "and similarly for the trace-free extrinsic curvature\n", "$$\n", "A_{i j} = e^{4 \\phi} \\bar{A}_{i j} \\; .\n", "$$\n", "It can be shown (Eq. (3.34) in Baumgarte & Shapiro (2010) with $\\alpha = -4$ and $\\psi = e^{\\phi}$) that the physical and conformal covariant derivatives obey\n", "$$\n", "D_{j} A^{i j} = e^{-10 \\phi} \\bar{D}_{j} \\left (e^{6 \\phi} \\bar{A}^{i j} \\right ) \\; .\n", "$$\n", "Then, the constraint becomes\n", "$$\n", "\\mathcal{M}^i = e^{-4\\phi} \\left(\n", "{\\underbrace {\\textstyle \\bar{D}_j \\bar{A}^{ij}}_{\\rm Term\\ 1}} + \n", "{\\underbrace {\\textstyle 6 \\bar{A}^{ij}\\partial_j \\phi}_{\\rm Term\\ 2}} - \n", "{\\underbrace {\\textstyle \\frac{2}{3} \\bar{\\gamma}^{ij}\\partial_j K}_{\\rm Term\\ 3}}\\right) \\; .\n", "$$\n", "\n", "Let's first implement Terms 2 and 3:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:21:03.445640Z", "iopub.status.busy": "2021-03-07T17:21:03.444942Z", "iopub.status.idle": "2021-03-07T17:21:03.447661Z", "shell.execute_reply": "2021-03-07T17:21:03.447144Z" } }, "outputs": [], "source": [ "# Step 3: M^i, the momentum constraint\n", "\n", "MU = ixp.zerorank1()\n", "\n", "# Term 2: 6 A^{ij} \\partial_j \\phi:\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " MU[i] += 6*Bq.AbarUU[i][j]*Bq.phi_dD[j]\n", "\n", "# Term 3: -2/3 \\bar{\\gamma}^{ij} K_{,j}\n", "trK_dD = ixp.declarerank1(\"trK_dD\") # Not defined in BSSN_RHSs; only trK_dupD is defined there.\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " MU[i] += -sp.Rational(2,3)*Bq.gammabarUU[i][j]*trK_dD[j]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we turn our attention to Term 1. The covariant divergence involves upper indices in $\\bar{A}^{i j}$, but it would be easier for us to finite difference the rescaled $\\bar{A}_{i j}$. A simple application of the inverse conformal metric yields\n", "$$\n", "\\bar{D}_{j} \\bar{A}^{i j} = \\bar{\\gamma}^{i k} \\bar{\\gamma}^{j l} \\bar{D}_{j} \\bar{A}_{k l} \\; .\n", "$$\n", "As usual, the covariant derivative is related to the ordinary derivative using the conformal Christoffel symbols\n", "$$\n", "\\bar{D}_{k} \\bar{A}_{i j} = \\partial_{k} \\bar{A}_{i j} - \\bar{\\Gamma}^{l}_{k i} \\bar{A}_{l j} - \\bar{\\Gamma}^{l}_{k j} \\bar{A}_{i l} \\; .\n", "$$\n", "It is the ordinary derivative above that is approximated by finite difference. The BSSN formulation used here does not rely on spatial derivatives $\\partial_{k} \\bar{A}_{i j}$ in any of the right-hand-sides (except for the advection term, which uses the upwinded derivative), and so we must declare a new ordinary, centered stencil derivative field of rank 3." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:21:03.522418Z", "iopub.status.busy": "2021-03-07T17:21:03.486570Z", "iopub.status.idle": "2021-03-07T17:21:03.620282Z", "shell.execute_reply": "2021-03-07T17:21:03.619601Z" } }, "outputs": [], "source": [ "# First define aDD_dD:\n", "aDD_dD = ixp.declarerank3(\"aDD_dD\",\"sym01\")\n", "\n", "# Then evaluate the conformal covariant derivative \\bar{D}_j \\bar{A}_{lm}\n", "AbarDD_dBarD = ixp.zerorank3()\n", "for i in range(DIM):\n", " for j in range(DIM):\n", " for k in range(DIM):\n", " AbarDD_dBarD[i][j][k] = Bq.AbarDD_dD[i][j][k]\n", " for l in range(DIM):\n", " AbarDD_dBarD[i][j][k] += -Bq.GammabarUDD[l][k][i]*Bq.AbarDD[l][j]\n", " AbarDD_dBarD[i][j][k] += -Bq.GammabarUDD[l][k][j]*Bq.AbarDD[i][l]\n", "\n", "# Term 1: Contract twice with the metric to make \\bar{D}_{j} \\bar{A}^{ij}\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", " MU[i] += Bq.gammabarUU[i][k]*Bq.gammabarUU[j][l]*AbarDD_dBarD[k][l][j]\n", "\n", "# Finally, we multiply by e^{-4 phi} and rescale the momentum constraint:\n", "for i in range(DIM):\n", " MU[i] *= Bq.exp_m4phi / rfm.ReU[i]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 4: Code Validation against `BSSN.BSSN_constraints` 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 the RHSs of the BSSN equations between\n", "1. this tutorial and \n", "2. the NRPy+ [BSSN.BSSN_constraints](../edit/BSSN/BSSN_constraints.py) module.\n", "\n", "By default, we analyze these expressions in Spherical coordinates, though other coordinate systems may be chosen." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:21:03.626554Z", "iopub.status.busy": "2021-03-07T17:21:03.625853Z", "iopub.status.idle": "2021-03-07T17:21:05.588244Z", "shell.execute_reply": "2021-03-07T17:21:05.587701Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Consistency check between BSSN_constraints tutorial and NRPy+ module: ALL SHOULD BE ZERO.\n", "H - bssncon.H = 0\n", "MU[0] - bssncon.MU[0] = 0\n", "MU[1] - bssncon.MU[1] = 0\n", "MU[2] - bssncon.MU[2] = 0\n" ] } ], "source": [ "# Step 4: Code Validation against BSSN.BSSN_constraints NRPy+ module\n", "\n", "# We already have SymPy expressions for BSSN constraints\n", "# in terms of other SymPy variables. Even if we reset the\n", "# list of NRPy+ gridfunctions, these *SymPy* expressions for\n", "# BSSN constraint variables *will remain unaffected*.\n", "#\n", "# Here, we will use the above-defined BSSN constraint expressions\n", "# to validate against the same expressions in the\n", "# BSSN/BSSN_constraints.py file, to ensure consistency between\n", "# this tutorial and the module itself.\n", "#\n", "# Reset the list of gridfunctions, as registering a gridfunction\n", "# twice (in the bssnrhs.BSSN_RHSs() call) will spawn an error.\n", "gri.glb_gridfcs_list = []\n", "\n", "# Call the BSSN_RHSs() function from within the\n", "# BSSN/BSSN_RHSs.py module,\n", "# which should do exactly the same as in Steps 1-16 above.\n", "import BSSN.BSSN_constraints as bssncon\n", "bssncon.BSSN_constraints()\n", "\n", "print(\"Consistency check between BSSN_constraints tutorial and NRPy+ module: ALL SHOULD BE ZERO.\")\n", "\n", "print(\"H - bssncon.H = \" + str(H - bssncon.H))\n", "for i in range(DIM):\n", " print(\"MU[\"+str(i)+\"] - bssncon.MU[\"+str(i)+\"] = \" + str(MU[i] - bssncon.MU[i]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\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\n", "[Tutorial-BSSN_constraints.pdf](Tutorial-BSSN_constraints.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": 6, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:21:05.593434Z", "iopub.status.busy": "2021-03-07T17:21:05.592342Z", "iopub.status.idle": "2021-03-07T17:21:08.879753Z", "shell.execute_reply": "2021-03-07T17:21:08.880335Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Created Tutorial-BSSN_constraints.tex, and compiled LaTeX file to PDF file\n", " Tutorial-BSSN_constraints.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_constraints\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.8.2" } }, "nbformat": 4, "nbformat_minor": 2 }