{
"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",
"# `BaikalETK`: NRPy+-Based BSSN Solver for the Einstein Toolkit\n",
"\n",
"## Author: Zach Etienne\n",
"\n",
"#### Special thanks to Roland Haas for help in answering many implementation questions\n",
"\n",
"## This module generates `BaikalETK`, an [Einstein Toolkit](https://einsteintoolkit.org) thorn for solving Einstein's equations in the BSSN formalism, in Cartesian coordinates. It features SIMD intrinsics and OpenMP support.\n",
"\n",
"**Notebook Status:** <font color='orange'><b> Validated against the Einstein Toolkit `McLachlan` BSSN thorn, both in the context of black hole binary simulations (excellent gravitational waveform agreement) as well as binary neutron star simulations (when parameter `add_stress_energy_source_terms` below is set to `True`). Once plots demonstrating this agreement are added to this tutorial notebook, the validation status will be set to</b></font> <font color='green'><b>Validated</b></font>.\n",
"\n",
"**Validation Notes:** This tutorial notebook has been validated against a trusted Einstein Toolkit thorn, but plots demonstrating its validation have yet to be included in this notebook.\n",
"\n",
"## Introduction\n",
"\n",
"```\n",
"How often did my soul cry out:\n",
"Come back to Baikal once again?\n",
"I still do not know this lake:\n",
"To see does not mean to know.\n",
"```\n",
"[Igor Severyanin](https://en.wikipedia.org/wiki/Igor_Severyanin), [[1]](https://1baikal.ru/en/istoriya/let’s-turn-to-baikal-a-poetic-view).\n",
"\n",
"[Lake Baikal](https://en.wikipedia.org/wiki/Lake_Baikal) is home to the [nerpa seal](https://en.wikipedia.org/wiki/Baikal_seal), NRPy+'s mascot.\n",
"\n",
"This thorn is meant to reproduce the functionality of the `McLachlan` thorn, generated by the [Mathematica](https://www.wolfram.com/mathematica/)-based [Kranc](http://kranccode.org/) code, but using the NRPy+ infrastructure.\n",
"\n",
"### Associated NRPy+ Source Code & Tutorial Modules for this module: \n",
"* [BSSN/ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear.py](../edit/BSSN/ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear.py); [\\[**tutorial**\\]](Tutorial-ADM_Initial_Data-Converting_Exact_ADM_Spherical_or_Cartesian_to_BSSNCurvilinear.ipynb): Spherical/Cartesian ADM$\\to$Curvilinear BSSN converter function, for which ADM quantities are assumed given at each gridpoint (i.e., exact, closed-form expressions are not given). This is used to generate BaikalETK's ADM$\\to$BSSN function, as in the ETK spacetime evolution modules are to assume that initial data are given as ADM quantities in the Cartesian basis at each gridpoint.\n",
"* [BSSN/ADM_in_terms_of_BSSN.py](../edit/BSSN/ADM_in_terms_of_BSSN.py); [\\[**tutorial**\\]](Tutorial-ADM_in_terms_of_BSSN.ipynb): Constructs ADM quantities in terms of BSSN quantities (in arbitrary curvilinear coordinates, though we use Cartesian here). This is used to generate BaikalETK's BSSN$\\to$ADM function, which make ADM variables available to diagnostic thorns within the ETK.\n",
"* [BSSN/BSSN_constraints.py](../edit/BSSN/BSSN_constraints.py); [\\[**tutorial**\\]](Tutorial-BSSN_constraints.ipynb): Hamiltonian constraint in BSSN curvilinear basis/coordinates\n",
"* [BSSN/BSSN_RHSs.py](../edit/BSSN/BSSN_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb): Generates the right-hand sides for the BSSN evolution equations in singular, curvilinear coordinates\n",
"* [BSSN/BSSN_gauge_RHSs.py](../edit/BSSN/BSSN_gauge_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb): Generates the right-hand sides for the BSSN gauge evolution equations in singular, curvilinear 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 needed Python/NRPy+ modules\n",
"1. [Step 2](#bssn): NRPy+-generated C code kernels for BSSN spacetime solve\n",
" 1. [Step 2.a](#bssnrhs): BSSN RHS expressions\n",
" 1. [Step 2.b](#hammomconstraints): Hamiltonian & momentum constraints\n",
" 1. [Step 2.c](#gamconstraint): Enforce conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint (in Cartesian coordinates, $\\det{\\hat{\\gamma}_{ij}}=1$)\n",
" 1. [Step 2.d](#parallel_codegen): Generate all the above C code kernels in parallel\n",
"1. [Step 3](#cclfiles): CCL files - Define how this module interacts and interfaces with the wider Einstein Toolkit infrastructure\n",
" 1. [Step 3.a](#paramccl): `param.ccl`: specify free parameters within `BaikalETK`\n",
" 1. [Step 3.b](#interfaceccl): `interface.ccl`: define needed gridfunctions; provide keywords denoting what this thorn provides and what it should inherit from other thorns\n",
" 1. [Step 3.c](#scheduleccl): `schedule.ccl`:schedule all functions used within `BaikalETK`, specify data dependencies within said functions, and allocate memory for gridfunctions\n",
"1. [Step 4](#cdrivers): C driver functions for ETK registration & NRPy+-generated kernels\n",
" 1. [Step 4.a](#etkfunctions): Needed ETK functions: Banner, Symmetry registration, Parameter sanity check, Method of Lines (`MoL`) registration, Boundary condition\n",
" 1. [Step 4.b](#bssnadmconversions): BSSN $\\leftrightarrow$ ADM conversions\n",
" 1. [Step 4.b.i](#admtobssn): ADM $\\to$ BSSN\n",
" 1. [Step 4.b.ii](#bssntoadm): BSSN $\\to$ ADM\n",
" 1. [Step 4.c](#bssnrhss) Evaluate BSSN right-hand-sides (RHSs)\n",
" 1. [Step 4.c.i](#ricci): Evaluate Ricci tensor\n",
" 1. [Step 4.c.ii](#bssnrhssricciinput): Evaluate BSSN RHSs, using Ricci tensor as input \n",
" 1. [Step 4.d](#enforcegammahatconstraint): Enforcing conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint (in Cartesian coordinates, $\\det{\\hat{\\gamma}_{ij}}=1$)\n",
" 1. [Step 4.e](#diagnostics): Diagnostics: Computing Hamiltonian & momentum constraints\n",
" 1. [Step 4.f](#t4uu): `driver_BSSN_T4UU()`: Compute $T^{\\mu\\nu}$ from `TmunuBase`'s $T_{\\mu\\nu}$\n",
" 1. [Step 4.g](#makecodedefn): `make.code.defn`: List of all C driver functions needed to compile `BaikalETK`\n",
"1. [Step 5](#code_validation): Code Validation\n",
" 1. [Step 5.a](#self_validation): Validation against [BaikalETK.BaikalETK_Pymodule](../edit/BaikalETK/BaikalETK_Pymodule.py) module\n",
"1. [Step 6](#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 Python/NRPy+ modules \\[Back to [top](#toc)\\]\n",
"\n",
"$$\\label{initializenrpy}$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Step 1: Import needed core NRPy+ modules\n",
"from outputC import * # 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 loop as lp # NRPy+: Generate C code loops\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 shutil, os, sys # Standard Python modules for multiplatform OS-level functions\n",
"\n",
"# Create directory for BaikalETK thorn & subdirectories in case they don't exist.\n",
"outrootdir = \"BaikalETK/\"\n",
"cmd.mkdir(os.path.join(outrootdir))\n",
"outdir = os.path.join(outrootdir,\"src\") # Main C code output directory\n",
"cmd.mkdir(outdir)\n",
"\n",
"# Set spatial dimension (must be 3 for BSSN)\n",
"DIM = 3\n",
"par.set_parval_from_str(\"grid::DIM\",DIM)\n",
"\n",
"# Enable stress-energy terms?\n",
"add_stress_energy_source_terms = False\n",
"\n",
"# Default Kreiss-Oliger dissipation strength\n",
"default_KO_strength = 0.1\n",
"\n",
"# Step 2: Set some core parameters, including CoordSystem MoL timestepping algorithm,\n",
"# FD order, floating point precision, and CFL factor:\n",
"# Choices are: Spherical, SinhSpherical, SinhSphericalv2, Cylindrical, SinhCylindrical, \n",
"# SymTP, SinhSymTP\n",
"# NOTE: Only CoordSystem == Cartesian makes sense here; new \n",
"# boundary conditions are needed within the ETK for \n",
"# Spherical, etc. coordinates.\n",
"CoordSystem = \"Cartesian\"\n",
"\n",
"par.set_parval_from_str(\"reference_metric::CoordSystem\",CoordSystem)\n",
"rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.\n",
"\n",
"# Set the standard 1+log lapse condition\n",
"LapseCondition = \"OnePlusLog\"\n",
"# Set the standard, second-order advecting-shift, Gamma-driving shift condition\n",
"ShiftCondition = \"GammaDriving2ndOrder_NoCovariant\"\n",
"\n",
"FD_order = 4 # Finite difference order: even numbers only, starting with 2. 12 is generally unstable\n",
"REAL = \"CCTK_REAL\" # Set REAL to CCTK_REAL, the ETK data type for \n",
" # floating point precision (typically `double`)\n",
"# Set finite differencing order:\n",
"par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", FD_order)\n",
"\n",
"# Copy SIMD/SIMD_intrinsics.h to $outdir/SIMD/SIMD_intrinsics.h\n",
"cmd.mkdir(os.path.join(outdir,\"SIMD\"))\n",
"shutil.copy(os.path.join(\"SIMD/\")+\"SIMD_intrinsics.h\",os.path.join(outdir,\"SIMD/\"))\n",
"\n",
"# Set the gridfunction memory access type to ETK-like, so that finite_difference\n",
"# knows how to read and write gridfunctions from/to memory.\n",
"par.set_parval_from_str(\"grid::GridFuncMemAccess\",\"ETK\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='bssn'></a>\n",
"\n",
"# Step 2: Output C code for BSSN spacetime solve \\[Back to [top](#toc)\\]\n",
"$$\\label{bssn}$$\n",
"\n",
"<a id='bssnrhs'></a>\n",
"\n",
"## Step 2.a: BSSN RHS expressions \\[Back to [top](#toc)\\]\n",
"$$\\label{bssnrhs}$$\n",
"\n",
"`BaikalETK` implements a fully covariant version of the BSSN 3+1 decomposition of Einstein's equations of general relativity, which is fully documented within NRPy+ ([start here](Tutorial-BSSN_formulation.ipynb)). However, especially if you are a newcomer to the field of numerical relativity, you may also find the following lectures and papers useful for understanding the adopted formalism:\n",
"\n",
"* Mathematical foundations of BSSN and 3+1 initial value problem decompositions of Einstein's equations:\n",
" * [Thomas Baumgarte's lectures on mathematical formulation of numerical relativity](https://www.youtube.com/watch?v=t3uo2R-yu4o&list=PLRVOWML3TL_djTd_nsTlq5aJjJET42Qke)\n",
" * [Yuichiro Sekiguchi's introduction to BSSN](http://www2.yukawa.kyoto-u.ac.jp/~yuichiro.sekiguchi/3+1.pdf) \n",
"* Extensions to the standard BSSN approach used in NRPy+\n",
" * [Brown's covariant \"Lagrangian\" formalism of BSSN](https://arxiv.org/abs/0902.3652)\n",
" * [BSSN in spherical coordinates, using the reference-metric approach of Baumgarte, Montero, Cordero-Carrión, and Müller (2012)](https://arxiv.org/abs/1211.6632)\n",
" * [BSSN in generic curvilinear coordinates, using the extended reference-metric approach of Ruchlin, Etienne, and Baumgarte (2018)](https://arxiv.org/abs/1712.07658)\n",
" \n",
"Here, we simply call the [BSSN.BSSN_RHSs](../edit/BSSN/BSSN_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_RHSs.ipynb) and [BSSN.BSSN_gauge_RHSs](../edit/BSSN/BSSN_gauge_RHSs.py); [\\[**tutorial**\\]](Tutorial-BSSN_time_evolution-BSSN_gauge_RHSs.ipynb) NRPy+ Python modules to generate the symbolic expressions, add Kreiss-Oliger dissipation, and then output the finite-difference C code form of the equations using NRPy+'s [finite_difference](../edit/finite_difference.py) ([**tutorial**](Tutorial-Finite_Difference_Derivatives.ipynb)) C code generation module."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Generating symbolic expressions for BSSN RHSs...\n",
"Finished BSSN symbolic expressions in 0.9867644309997559 seconds.\n"
]
}
],
"source": [
"import time # Standard Python module; useful for benchmarking below expression & code generation.\n",
"\n",
"import BSSN.BSSN_RHSs as rhs\n",
"import BSSN.BSSN_gauge_RHSs as gaugerhs\n",
"par.set_parval_from_str(\"BSSN.BSSN_gauge_RHSs::ShiftEvolutionOption\", ShiftCondition)\n",
"par.set_parval_from_str(\"BSSN.BSSN_gauge_RHSs::LapseEvolutionOption\", LapseCondition)\n",
"\n",
"print(\"Generating symbolic expressions for BSSN RHSs...\")\n",
"start = time.time()\n",
"# Enable rfm_precompute infrastructure, which results in \n",
"# BSSN RHSs that are free of transcendental functions,\n",
"# even in curvilinear coordinates, so long as \n",
"# ConformalFactor is set to \"W\" (default).\n",
"cmd.mkdir(os.path.join(outdir,\"rfm_files/\"))\n",
"par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\",\"True\")\n",
"par.set_parval_from_str(\"reference_metric::rfm_precompute_Ccode_outdir\",os.path.join(outdir,\"rfm_files/\"))\n",
"\n",
"# Evaluate BSSN + BSSN gauge RHSs with rfm_precompute enabled:\n",
"import BSSN.BSSN_quantities as Bq\n",
"par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\",\"True\")\n",
"\n",
"rhs.BSSN_RHSs()\n",
"\n",
"if add_stress_energy_source_terms == True:\n",
" T4UU = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"T4UU\",\"sym01\",DIM=4)\n",
" import BSSN.BSSN_stress_energy_source_terms as Bsest\n",
" Bsest.BSSN_source_terms_for_BSSN_RHSs(T4UU)\n",
" rhs.trK_rhs += Bsest.sourceterm_trK_rhs\n",
" for i in range(DIM):\n",
" # Needed for Gamma-driving shift RHSs:\n",
" rhs.Lambdabar_rhsU[i] += Bsest.sourceterm_Lambdabar_rhsU[i]\n",
" # Needed for BSSN RHSs:\n",
" rhs.lambda_rhsU[i] += Bsest.sourceterm_lambda_rhsU[i]\n",
" for j in range(DIM):\n",
" rhs.a_rhsDD[i][j] += Bsest.sourceterm_a_rhsDD[i][j]\n",
"\n",
"gaugerhs.BSSN_gauge_RHSs()\n",
"\n",
"# Add Kreiss-Oliger dissipation to the BSSN RHSs:\n",
"thismodule = \"KO_Dissipation\"\n",
"diss_strength = par.Cparameters(\"REAL\", thismodule, \"diss_strength\", default_KO_strength)\n",
"\n",
"alpha_dKOD = ixp.declarerank1(\"alpha_dKOD\")\n",
"cf_dKOD = ixp.declarerank1(\"cf_dKOD\")\n",
"trK_dKOD = ixp.declarerank1(\"trK_dKOD\")\n",
"betU_dKOD = ixp.declarerank2(\"betU_dKOD\",\"nosym\")\n",
"vetU_dKOD = ixp.declarerank2(\"vetU_dKOD\",\"nosym\")\n",
"lambdaU_dKOD = ixp.declarerank2(\"lambdaU_dKOD\",\"nosym\")\n",
"aDD_dKOD = ixp.declarerank3(\"aDD_dKOD\",\"sym01\")\n",
"hDD_dKOD = ixp.declarerank3(\"hDD_dKOD\",\"sym01\")\n",
"for k in range(DIM):\n",
" gaugerhs.alpha_rhs += diss_strength*alpha_dKOD[k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]\n",
" rhs.cf_rhs += diss_strength* cf_dKOD[k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]\n",
" rhs.trK_rhs += diss_strength* trK_dKOD[k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]\n",
" for i in range(DIM):\n",
" if \"2ndOrder\" in ShiftCondition:\n",
" gaugerhs.bet_rhsU[i] += diss_strength* betU_dKOD[i][k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]\n",
" gaugerhs.vet_rhsU[i] += diss_strength* vetU_dKOD[i][k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]\n",
" rhs.lambda_rhsU[i] += diss_strength*lambdaU_dKOD[i][k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]\n",
" for j in range(DIM):\n",
" rhs.a_rhsDD[i][j] += diss_strength*aDD_dKOD[i][j][k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]\n",
" rhs.h_rhsDD[i][j] += diss_strength*hDD_dKOD[i][j][k]*rfm.ReU[k] # ReU[k] = 1/scalefactor_orthog_funcform[k]\n",
"\n",
"# We use betaU as our upwinding control vector:\n",
"Bq.BSSN_basic_tensors()\n",
"betaU = Bq.betaU\n",
"\n",
"import BSSN.Enforce_Detgammabar_Constraint as EGC\n",
"enforce_detg_constraint_symb_expressions = EGC.Enforce_Detgammabar_Constraint_symb_expressions()\n",
"\n",
"# Next compute Ricci tensor\n",
"par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\",\"False\")\n",
"Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()\n",
"\n",
"# Now that we are finished with all the rfm hatted\n",
"# quantities in generic precomputed functional\n",
"# form, let's restore them to their closed-\n",
"# form expressions.\n",
"par.set_parval_from_str(\"reference_metric::enable_rfm_precompute\",\"False\") # Reset to False to disable rfm_precompute.\n",
"rfm.ref_metric__hatted_quantities()\n",
"end = time.time()\n",
"print(\"Finished BSSN symbolic expressions in \"+str(end-start)+\" seconds.\")\n",
"\n",
"\n",
"def BSSN_RHSs():\n",
" print(\"Generating C code for BSSN RHSs in \"+par.parval_from_str(\"reference_metric::CoordSystem\")+\" coordinates.\")\n",
" start = time.time()\n",
"\n",
" BSSN_evol_rhss = [ \\\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"aDD00\"),rhs=rhs.a_rhsDD[0][0]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"aDD01\"),rhs=rhs.a_rhsDD[0][1]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"aDD02\"),rhs=rhs.a_rhsDD[0][2]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"aDD11\"),rhs=rhs.a_rhsDD[1][1]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"aDD12\"),rhs=rhs.a_rhsDD[1][2]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"aDD22\"),rhs=rhs.a_rhsDD[2][2]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"alpha\"),rhs=gaugerhs.alpha_rhs),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"betU0\"),rhs=gaugerhs.bet_rhsU[0]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"betU1\"),rhs=gaugerhs.bet_rhsU[1]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"betU2\"),rhs=gaugerhs.bet_rhsU[2]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"cf\"), rhs=rhs.cf_rhs),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"hDD00\"),rhs=rhs.h_rhsDD[0][0]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"hDD01\"),rhs=rhs.h_rhsDD[0][1]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"hDD02\"),rhs=rhs.h_rhsDD[0][2]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"hDD11\"),rhs=rhs.h_rhsDD[1][1]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"hDD12\"),rhs=rhs.h_rhsDD[1][2]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"hDD22\"),rhs=rhs.h_rhsDD[2][2]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"lambdaU0\"),rhs=rhs.lambda_rhsU[0]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"lambdaU1\"),rhs=rhs.lambda_rhsU[1]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"lambdaU2\"),rhs=rhs.lambda_rhsU[2]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"trK\"), rhs=rhs.trK_rhs),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"vetU0\"),rhs=gaugerhs.vet_rhsU[0]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"vetU1\"),rhs=gaugerhs.vet_rhsU[1]),\n",
" lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"vetU2\"),rhs=gaugerhs.vet_rhsU[2]) ]\n",
"\n",
" BSSN_RHSs_string = fin.FD_outputC(\"returnstring\",BSSN_evol_rhss, params=\"outCverbose=False,SIMD_enable=True\",\n",
" upwindcontrolvec=betaU)\n",
"\n",
" with open(os.path.join(outdir,\"BSSN_RHSs.h\"), \"w\") as file:\n",
" file.write(lp.loop([\"i2\",\"i1\",\"i0\"],[\"cctk_nghostzones[2]\",\"cctk_nghostzones[1]\",\"cctk_nghostzones[0]\"],\n",
" [\"cctk_lsh[2]-cctk_nghostzones[2]\",\"cctk_lsh[1]-cctk_nghostzones[1]\",\"cctk_lsh[0]-cctk_nghostzones[0]\"],\n",
" [\"1\",\"1\",\"SIMD_width\"],\n",
" [\"#pragma omp parallel for\",\n",
" \"#include \\\"rfm_files/rfm_struct__SIMD_outer_read2.h\\\"\",\n",
" \"#include \\\"rfm_files/rfm_struct__SIMD_outer_read1.h\\\"\"],\"\",\n",
" \"#include \\\"rfm_files/rfm_struct__SIMD_inner_read0.h\\\"\\n\"+BSSN_RHSs_string))\n",
" end = time.time()\n",
" print(\"Finished BSSN_RHS C codegen in \" + str(end - start) + \" seconds.\")\n",
"\n",
"def Ricci():\n",
" print(\"Generating C code for Ricci tensor in \"+par.parval_from_str(\"reference_metric::CoordSystem\")+\" coordinates.\")\n",
" start = time.time()\n",
" Ricci_string = fin.FD_outputC(\"returnstring\",\n",
" [lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD00\"),rhs=Bq.RbarDD[0][0]),\n",
" lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD01\"),rhs=Bq.RbarDD[0][1]),\n",
" lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD02\"),rhs=Bq.RbarDD[0][2]),\n",
" lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD11\"),rhs=Bq.RbarDD[1][1]),\n",
" lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD12\"),rhs=Bq.RbarDD[1][2]),\n",
" lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"RbarDD22\"),rhs=Bq.RbarDD[2][2])],\n",
" params=\"outCverbose=False,SIMD_enable=True\")\n",
" with open(os.path.join(outdir,\"BSSN_Ricci.h\"), \"w\") as file:\n",
" file.write(lp.loop([\"i2\",\"i1\",\"i0\"],[\"cctk_nghostzones[2]\",\"cctk_nghostzones[1]\",\"cctk_nghostzones[0]\"],\n",
" [\"cctk_lsh[2]-cctk_nghostzones[2]\",\"cctk_lsh[1]-cctk_nghostzones[1]\",\"cctk_lsh[0]-cctk_nghostzones[0]\"],\n",
" [\"1\",\"1\",\"SIMD_width\"],\n",
" [\"#pragma omp parallel for\",\n",
" \"#include \\\"rfm_files/rfm_struct__SIMD_outer_read2.h\\\"\",\n",
" \"#include \\\"rfm_files/rfm_struct__SIMD_outer_read1.h\\\"\"],\"\",\n",
" \"#include \\\"rfm_files/rfm_struct__SIMD_inner_read0.h\\\"\\n\"+Ricci_string))\n",
" end = time.time()\n",
" print(\"Finished Ricci C codegen in \" + str(end - start) + \" seconds.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='hammomconstraints'></a>\n",
"\n",
"## Step 2.b: Hamiltonian & momentum constraints \\[Back to [top](#toc)\\]\n",
"$$\\label{hammomconstraints}$$\n",
"\n",
"Next output the C code for evaluating the Hamiltonian & momentum constraints [(**Tutorial**)](Tutorial-BSSN_constraints.ipynb). In the absence of numerical error, this constraint should evaluate to zero. However it does not due to numerical (typically truncation and roundoff) error. Therefore it is useful to measure the Hamiltonian & momentum constraint violation to gauge the accuracy of our simulation, and, ultimately determine whether errors are dominated by numerical finite differencing (truncation) error as expected."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# First register the Hamiltonian as a gridfunction.\n",
"H = gri.register_gridfunctions(\"AUX\",\"H\")\n",
"MU = ixp.register_gridfunctions_for_single_rank1(\"AUX\", \"MU\")\n",
"\n",
"# Then define the Hamiltonian constraint and output the optimized C code.\n",
"import BSSN.BSSN_constraints as bssncon\n",
"\n",
"def BSSNconstraints():\n",
" bssncon.BSSN_constraints(add_T4UUmunu_source_terms=False)\n",
" if add_stress_energy_source_terms == True:\n",
"# T4UU = gri.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"T4UU\",\"sym01\",DIM=4)\n",
" import BSSN.BSSN_stress_energy_source_terms as Bsest\n",
" Bsest.BSSN_source_terms_for_BSSN_constraints(T4UU)\n",
" bssncon.H += Bsest.sourceterm_H\n",
" for i in range(DIM):\n",
" bssncon.MU[i] += Bsest.sourceterm_MU[i]\n",
" \n",
" start = time.time()\n",
" print(\"Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.\")\n",
" Ham_mom_string = fin.FD_outputC(\"returnstring\", \n",
" [lhrh(lhs=gri.gfaccess(\"aux_gfs\", \"H\"), rhs=bssncon.H),\n",
" lhrh(lhs=gri.gfaccess(\"aux_gfs\", \"MU0\"), rhs=bssncon.MU[0]),\n",
" lhrh(lhs=gri.gfaccess(\"aux_gfs\", \"MU1\"), rhs=bssncon.MU[1]),\n",
" lhrh(lhs=gri.gfaccess(\"aux_gfs\", \"MU2\"), rhs=bssncon.MU[2])],\n",
" params=\"outCverbose=False\")\n",
"\n",
" with open(os.path.join(outdir,\"BSSN_constraints.h\"), \"w\") as file:\n",
" file.write(lp.loop([\"i2\",\"i1\",\"i0\"],[\"cctk_nghostzones[2]\",\"cctk_nghostzones[1]\",\"cctk_nghostzones[0]\"],\n",
" [\"cctk_lsh[2]-cctk_nghostzones[2]\",\"cctk_lsh[1]-cctk_nghostzones[1]\",\"cctk_lsh[0]-cctk_nghostzones[0]\"],\n",
" [\"1\",\"1\",\"1\"],[\"#pragma omp parallel for\",\"\",\"\"], \"\", Ham_mom_string))\n",
" end = time.time()\n",
" print(\"Finished Hamiltonian & momentum constraint C codegen in \" + str(end - start) + \" seconds.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='gamconstraint'></a>\n",
"\n",
"## Step 2.c: Enforce conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint (in Cartesian coordinates, $\\det{\\hat{\\gamma}_{ij}}=1$) \\[Back to [top](#toc)\\]\n",
"$$\\label{gamconstraint}$$\n",
"\n",
"Then 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)), as [documented in the corresponding NRPy+ tutorial notebook](Tutorial-BSSN-Enforcing_Determinant_gammabar_equals_gammahat_Constraint.ipynb)\n",
"\n",
"Applying curvilinear boundary conditions should affect the initial data at the outer boundary, and will in general cause the $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint to be violated there. Thus after we apply these boundary conditions, we must always call the routine for enforcing the $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def gammadet():\n",
" start = time.time()\n",
" print(\"Generating optimized C code for gamma constraint. May take a while, depending on CoordSystem.\")\n",
" enforce_gammadet_string = fin.FD_outputC(\"returnstring\", enforce_detg_constraint_symb_expressions,\n",
" params=\"outCverbose=False,preindent=0,includebraces=False\")\n",
"\n",
" with open(os.path.join(outdir,\"enforce_detgammabar_constraint.h\"), \"w\") as file:\n",
" file.write(lp.loop([\"i2\",\"i1\",\"i0\"],[\"0\", \"0\", \"0\"],\n",
" [\"cctk_lsh[2]\",\"cctk_lsh[1]\",\"cctk_lsh[0]\"],\n",
" [\"1\",\"1\",\"1\"],\n",
" [\"#pragma omp parallel for\",\n",
" \"#include \\\"rfm_files/rfm_struct__read2.h\\\"\",\n",
" \"#include \\\"rfm_files/rfm_struct__read1.h\\\"\"],\"\",\n",
" \"#include \\\"rfm_files/rfm_struct__read0.h\\\"\\n\"+enforce_gammadet_string))\n",
" end = time.time()\n",
" print(\"Finished gamma constraint C codegen in \" + str(end - start) + \" seconds.\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='parallel_codegen'></a>\n",
"\n",
"## Step 2.d: Generate all C codes in parallel \\[Back to [top](#toc)\\]\n",
"$$\\label{parallel_codegen}$$\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Generating C code for Ricci tensor in Cartesian coordinates.\n",
"Generating C code for BSSN RHSs in Cartesian coordinates.\n",
"Generating optimized C code for gamma constraint. May take a while, depending on CoordSystem.\n",
"Finished gamma constraint C codegen in 0.05950760841369629 seconds.\n",
"Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.\n",
"Finished Hamiltonian & momentum constraint C codegen in 5.987943410873413 seconds.\n",
"Finished BSSN_RHS C codegen in 16.53508186340332 seconds.\n",
"Finished Ricci C codegen in 17.892887353897095 seconds.\n"
]
}
],
"source": [
"# Step 0: Import the multiprocessing module.\n",
"import multiprocessing\n",
"\n",
"# Step 1: Create a list of functions we wish to evaluate in parallel\n",
"funcs = [BSSN_RHSs,Ricci,BSSNconstraints,gammadet]\n",
"# Step 1.a: Define master function for parallelization.\n",
"# Note that lambdifying this doesn't work in Python 3\n",
"def master_func(arg):\n",
" funcs[arg]()\n",
"\n",
"# Step 2: Evaluate list of functions in parallel if allowed; \n",
"# otherwise fallback to serial evaluation:\n",
"try:\n",
" if __name__ == '__main__':\n",
" pool = multiprocessing.Pool()\n",
" pool.map(master_func,range(len(funcs)))\n",
"except:\n",
" # If multiprocessing didn't work for whatever reason,\n",
" # evaluate functions in serial.\n",
" for func in funcs:\n",
" func()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='cclfiles'></a>\n",
"\n",
"# Step 3: ETK `ccl` file generation \\[Back to [top](#toc)\\]\n",
"$$\\label{cclfiles}$$\n",
"\n",
"<a id='paramccl'></a>\n",
"\n",
"## Step 3.a: `param.ccl`: specify free parameters within `BaikalETK` \\[Back to [top](#toc)\\]\n",
"$$\\label{paramccl}$$\n",
"\n",
"All parameters necessary for the computation of the BSSN right-hand side (RHS) expressions are registered within NRPy+; we use this information to automatically generate `param.ccl`. NRPy+ also specifies default values for each parameter. \n",
"\n",
"More information on `param.ccl` syntax can be found in the [official Einstein Toolkit documentation](http://cactuscode.org/documentation/referencemanual/ReferenceManualch8.html#x12-265000C2.3)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def keep_param__return_type(paramtuple):\n",
" keep_param = True # We'll not set some parameters in param.ccl; \n",
" # e.g., those that should be #define'd like M_PI.\n",
" typestring = \"\"\n",
" # Separate thorns within the ETK take care of grid/coordinate parameters;\n",
" # thus we ignore NRPy+ grid/coordinate parameters:\n",
" if paramtuple.module == \"grid\" or paramtuple.module == \"reference_metric\":\n",
" keep_param = False\n",
"\n",
" partype = paramtuple.type\n",
" if partype == \"bool\":\n",
" typestring += \"BOOLEAN \"\n",
" elif partype == \"REAL\":\n",
" if paramtuple.defaultval != 1e300: # 1e300 is a magic value indicating that the C parameter should be mutable\n",
" typestring += \"CCTK_REAL \"\n",
" else:\n",
" keep_param = False\n",
" elif partype == \"int\":\n",
" typestring += \"CCTK_INT \"\n",
" elif partype == \"#define\":\n",
" keep_param = False\n",
" elif partype == \"char\":\n",
" # FIXME: char/string parameter types should in principle be supported\n",
" print(\"Error: parameter \"+paramtuple.module+\"::\"+paramtuple.parname+\n",
" \" has unsupported type: \\\"\"+ paramtuple.type + \"\\\"\")\n",
" sys.exit(1)\n",
" else:\n",
" print(\"Error: parameter \"+paramtuple.module+\"::\"+paramtuple.parname+\n",
" \" has unsupported type: \\\"\"+ paramtuple.type + \"\\\"\")\n",
" sys.exit(1)\n",
" return keep_param, typestring\n",
"\n",
"with open(os.path.join(outrootdir,\"param.ccl\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"# This param.ccl file was automatically generated by NRPy+. \n",
"# You are advised against modifying it directly.\n",
"\n",
"shares: ADMBase\n",
"USES CCTK_INT lapse_timelevels # Needed to ensure ADMBase gridfunctions are allocated (see top of schedule.ccl)\n",
"USES CCTK_INT shift_timelevels # Needed to ensure ADMBase gridfunctions are allocated (see top of schedule.ccl)\n",
"USES CCTK_INT metric_timelevels # Needed to ensure ADMBase gridfunctions are allocated (see top of schedule.ccl)\n",
"\n",
"EXTENDS CCTK_KEYWORD evolution_method \"evolution_method\"\n",
"{\n",
" \"BaikalETK\" :: \"\"\n",
"} \n",
"\n",
"EXTENDS CCTK_KEYWORD lapse_evolution_method \"lapse_evolution_method\"\n",
"{\n",
" \"BaikalETK\" :: \"\"\n",
"} \n",
"\n",
"EXTENDS CCTK_KEYWORD shift_evolution_method \"shift_evolution_method\"\n",
"{\n",
" \"BaikalETK\" :: \"\"\n",
"} \n",
"\n",
"EXTENDS CCTK_KEYWORD dtshift_evolution_method \"dtshift_evolution_method\"\n",
"{\n",
" \"BaikalETK\" :: \"\"\n",
"} \n",
"\n",
"EXTENDS CCTK_KEYWORD dtlapse_evolution_method \"dtlapse_evolution_method\"\n",
"{\n",
" \"BaikalETK\" :: \"\"\n",
"} \n",
"\n",
"restricted:\n",
"\"\"\")\n",
" paramccl_str = \"\"\n",
" for i in range(len(par.glb_Cparams_list)):\n",
" # keep_param is a boolean indicating whether we should accept or reject\n",
" # the parameter. singleparstring will contain the string indicating\n",
" # the variable type.\n",
" keep_param, singleparstring = keep_param__return_type(par.glb_Cparams_list[i])\n",
"\n",
" if keep_param:\n",
" parname = par.glb_Cparams_list[i].parname\n",
" partype = par.glb_Cparams_list[i].type\n",
" singleparstring += parname + \" \\\"\"+ parname +\" (see NRPy+ for parameter definition)\\\"\\n\"\n",
" singleparstring += \"{\\n\"\n",
" if partype != \"bool\":\n",
" singleparstring += \" *:* :: \\\"All values accepted. NRPy+ does not restrict the allowed ranges of parameters yet.\\\"\\n\"\n",
" singleparstring += \"} \"+str(par.glb_Cparams_list[i].defaultval)+\"\\n\\n\"\n",
" \n",
" paramccl_str += singleparstring\n",
" file.write(paramccl_str)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='interfaceccl'></a>\n",
"\n",
"## Step 3.b: `interface.ccl`: define needed gridfunctions; provide keywords denoting what this thorn provides and what it should inherit from other thorns \\[Back to [top](#toc)\\]\n",
"$$\\label{interfaceccl}$$\n",
"\n",
"`interface.ccl` declares all gridfunctions and determines how `BaikalETK` interacts with other Einstein Toolkit thorns.\n",
"\n",
"The [official Einstein Toolkit (Cactus) documentation](http://cactuscode.org/documentation/referencemanual/ReferenceManual.html) defines what must/should be included in an `interface.ccl` file [**here**](http://cactuscode.org/documentation/referencemanual/ReferenceManualch8.html#x12-260000C2.2). "
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# First construct lists of the basic gridfunctions used in NRPy+. \n",
"# Each type will be its own group in BaikalETK.\n",
"evol_gfs_list = []\n",
"auxevol_gfs_list = []\n",
"aux_gfs_list = []\n",
"for i in range(len(gri.glb_gridfcs_list)):\n",
" if gri.glb_gridfcs_list[i].gftype == \"EVOL\":\n",
" evol_gfs_list.append( gri.glb_gridfcs_list[i].name+\"GF\")\n",
"\n",
" if gri.glb_gridfcs_list[i].gftype == \"AUX\":\n",
" aux_gfs_list.append( gri.glb_gridfcs_list[i].name+\"GF\")\n",
"\n",
" if gri.glb_gridfcs_list[i].gftype == \"AUXEVOL\":\n",
" auxevol_gfs_list.append(gri.glb_gridfcs_list[i].name+\"GF\")\n",
"\n",
"# NRPy+'s finite-difference code generator assumes gridfunctions\n",
"# are alphabetized; not sorting may result in unnecessary\n",
"# cache misses.\n",
"evol_gfs_list.sort()\n",
"aux_gfs_list.sort()\n",
"auxevol_gfs_list.sort()\n",
" \n",
"with open(os.path.join(outrootdir,\"interface.ccl\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"# With \"implements\", we give our thorn its unique name.\n",
"implements: BaikalETK\n",
"\n",
"# By \"inheriting\" other thorns, we tell the Toolkit that we \n",
"# will rely on variables/function that exist within those\n",
"# functions. \n",
"inherits: ADMBase Boundary Grid MethodofLines\\n\"\"\")\n",
" if add_stress_energy_source_terms == True:\n",
" file.write(\"inherits: TmunuBase\")\n",
" file.write(\"\"\"\n",
"# Needed functions and #include's:\n",
"USES INCLUDE: Symmetry.h\n",
"USES INCLUDE: Boundary.h\n",
"\n",
"# Needed Method of Lines function\n",
"CCTK_INT FUNCTION MoLRegisterEvolvedGroup(CCTK_INT IN EvolvedIndex, \\\n",
" CCTK_INT IN RHSIndex)\n",
"REQUIRES FUNCTION MoLRegisterEvolvedGroup\n",
"\n",
"# Needed Boundary Conditions function\n",
"CCTK_INT FUNCTION GetBoundarySpecification(CCTK_INT IN size, CCTK_INT OUT ARRAY nboundaryzones, CCTK_INT OUT ARRAY is_internal, CCTK_INT OUT ARRAY is_staggered, CCTK_INT OUT ARRAY shiftout)\n",
"USES FUNCTION GetBoundarySpecification\n",
"\n",
"CCTK_INT FUNCTION SymmetryTableHandleForGrid(CCTK_POINTER_TO_CONST IN cctkGH)\n",
"USES FUNCTION SymmetryTableHandleForGrid\n",
"\n",
"CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, CCTK_INT IN faces, CCTK_INT IN boundary_width, CCTK_INT IN table_handle, CCTK_STRING IN var_name, CCTK_STRING IN bc_name)\n",
"USES FUNCTION Boundary_SelectVarForBC\n",
"\n",
"# Needed for EinsteinEvolve/NewRad outer boundary condition driver:\n",
"CCTK_INT FUNCTION \\\\\n",
" NewRad_Apply \\\\\n",
" (CCTK_POINTER_TO_CONST IN cctkGH, \\\\\n",
" CCTK_REAL ARRAY IN var, \\\\\n",
" CCTK_REAL ARRAY INOUT rhs, \\\\\n",
" CCTK_REAL IN var0, \\\\\n",
" CCTK_REAL IN v0, \\\\\n",
" CCTK_INT IN radpower)\n",
"REQUIRES FUNCTION NewRad_Apply\n",
"\n",
"# Needed to convert ADM initial data into BSSN initial data (gamma extrapolation)\n",
"CCTK_INT FUNCTION \\\\\n",
" ExtrapolateGammas \\\\\n",
" (CCTK_POINTER_TO_CONST IN cctkGH, \\\\\n",
" CCTK_REAL ARRAY INOUT var)\n",
"REQUIRES FUNCTION ExtrapolateGammas\n",
"\n",
"# Tell the Toolkit that we want all gridfunctions\n",
"# to be visible to other thorns by using \n",
"# the keyword \"public\". Note that declaring these \n",
"# gridfunctions *does not* allocate memory for them;\n",
"# that is done by the schedule.ccl file.\n",
"\n",
"# FIXME: add info for symmetry conditions: \n",
"# https://einsteintoolkit.org/thornguide/CactusBase/SymBase/documentation.html\n",
"public:\n",
"\"\"\")\n",
" \n",
" # Next we declare gridfunctions based on their corresponding gridfunction groups as registered within NRPy+\n",
" \n",
" def output_list_of_gfs(gfs_list,description=\"User did not provide description\"):\n",
" gfsstr = \" \"\n",
" for i in range(len(gfs_list)):\n",
" gfsstr += gfs_list[i]\n",
" if i != len(gfs_list)-1:\n",
" gfsstr += \",\" # This is a comma-separated list of gridfunctions\n",
" else:\n",
" gfsstr += \"\\n} \\\"\"+description+\"\\\"\\n\\n\"\n",
" return gfsstr\n",
" # First EVOL type:\n",
" file.write(\"CCTK_REAL evol_variables type = GF Timelevels=3\\n{\\n\")\n",
" file.write(output_list_of_gfs(evol_gfs_list,\"BSSN evolved gridfunctions\"))\n",
" # Second EVOL right-hand-sides\n",
" file.write(\"CCTK_REAL evol_variables_rhs type = GF Timelevels=1 TAGS=\\'InterpNumTimelevels=1 prolongation=\\\"none\\\"\\'\\n{\\n\")\n",
" rhs_list = []\n",
" for gf in evol_gfs_list:\n",
" rhs_list.append(gf.replace(\"GF\",\"\")+\"_rhsGF\")\n",
" file.write(output_list_of_gfs(rhs_list,\"right-hand-side storage for BSSN evolved gridfunctions\"))\n",
" # Then AUX type:\n",
" file.write(\"CCTK_REAL aux_variables type = GF Timelevels=3\\n{\\n\")\n",
" file.write(output_list_of_gfs(aux_gfs_list,\"Auxiliary gridfunctions for BSSN diagnostics\"))\n",
" # Finally, AUXEVOL type:\n",
" file.write(\"CCTK_REAL auxevol_variables type = GF Timelevels=1 TAGS=\\'InterpNumTimelevels=1 prolongation=\\\"none\\\"\\'\\n{\\n\")\n",
" file.write(output_list_of_gfs(auxevol_gfs_list,\"Auxiliary gridfunctions needed for evaluating the BSSN RHSs\"))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='scheduleccl'></a>\n",
"\n",
"## Step 3.c: `schedule.ccl`: schedule all functions used within `BaikalETK`, specify data dependencies within said functions, and allocate memory for gridfunctions \\[Back to [top](#toc)\\]\n",
"$$\\label{scheduleccl}$$\n",
"\n",
"Official documentation on constructing ETK `schedule.ccl` files is found [here](http://cactuscode.org/documentation/referencemanual/ReferenceManualch8.html#x12-268000C2.4)."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"with open(os.path.join(outrootdir,\"schedule.ccl\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"# First allocate storage for all ADMBase gridfunctions, which are needed by NRPy+\n",
"STORAGE: ADMBase::metric[metric_timelevels], ADMBase::curv[metric_timelevels], ADMBase::lapse[lapse_timelevels], ADMBase::shift[shift_timelevels]\n",
"\n",
"# Next allocate storage for all 3 gridfunction groups used in BaikalETK\n",
"STORAGE: evol_variables[3] # Evolution variables\n",
"STORAGE: evol_variables_rhs[1] # Variables storing right-hand-sides\n",
"STORAGE: aux_variables[3] # Diagnostics variables\n",
"STORAGE: auxevol_variables[1] # Single-timelevel storage of variables needed for evolutions.\n",
"\n",
"# The following scheduler is based on Lean/LeanBSSNMoL/schedule.ccl\n",
"\n",
"schedule BaikalETK_Banner at STARTUP\n",
"{\n",
" LANG: C\n",
" OPTIONS: meta\n",
"} \"Output ASCII art banner\"\n",
"\n",
"schedule BaikalETK_RegisterSlicing at STARTUP after BaikalETK_Banner\n",
"{\n",
" LANG: C\n",
" OPTIONS: meta\n",
"} \"Register 3+1 slicing condition\"\n",
"\n",
"schedule BaikalETK_Symmetry_registration at BASEGRID\n",
"{\n",
" LANG: C\n",
" OPTIONS: Global\n",
"} \"Register symmetries, the CartGrid3D way.\"\n",
"\n",
"schedule BaikalETK_zero_rhss at BASEGRID after BaikalETK_Symmetry_registration\n",
"{\n",
" LANG: C\n",
"} \"Idea from Lean: set all rhs functions to zero to prevent spurious nans\"\n",
"\n",
"schedule BaikalETK_ADM_to_BSSN at CCTK_INITIAL after ADMBase_PostInitial\n",
"{\n",
" LANG: C\n",
" OPTIONS: Local\n",
" SYNC: evol_variables\n",
"} \"Convert initial data into BSSN variables\"\n",
"\n",
"schedule GROUP ApplyBCs as BaikalETK_ApplyBCs at CCTK_INITIAL after BaikalETK_ADM_to_BSSN\n",
"{\n",
"} \"Apply boundary conditions\"\n",
"\n",
"\n",
"# MoL: registration\n",
"\n",
"schedule BaikalETK_MoL_registration in MoL_Register\n",
"{\n",
" LANG: C\n",
" OPTIONS: META\n",
"} \"Register variables for MoL\"\n",
"\n",
"\n",
"# MoL: compute RHSs, etc\n",
"\"\"\")\n",
" if add_stress_energy_source_terms == True:\n",
" file.write(\"\"\"\n",
"schedule driver_BSSN_T4UU in MoL_CalcRHS as BaikalETK_T4UU before BaikalETK_BSSN_to_ADM\n",
"{\n",
" LANG: C\n",
"} \"MoL: Compute T4UU, needed for BSSN RHSs.\"\n",
"\n",
"schedule BaikalETK_BSSN_to_ADM in MoL_CalcRHS after BaikalETK_T4UU before BaikalETK_Ricci\n",
"{\n",
" LANG: C\n",
"} \"Perform BSSN-to-ADM conversion. Needed for HydroBase coupling.\"\n",
"\"\"\")\n",
" file.write(\"\"\"\n",
"schedule driver_pt1_BSSN_Ricci in MoL_CalcRHS as BaikalETK_Ricci before BaikalETK_RHS\n",
"{\n",
" LANG: C\n",
"} \"MoL: Compute Ricci tensor\"\n",
"\n",
"schedule driver_pt2_BSSN_RHSs in MoL_CalcRHS as BaikalETK_RHS after BaikalETK_Ricci\n",
"{\n",
" LANG: C\n",
"} \"MoL: Evaluate BSSN RHSs\"\n",
"\n",
"schedule BaikalETK_NewRad in MoL_CalcRHS after BaikalETK_RHS\n",
"{\n",
" LANG: C\n",
"} \"NewRad boundary conditions, scheduled right after RHS eval.\"\n",
"\n",
"schedule enforce_detgammabar_constraint in MoL_PostStep before BC_Update\n",
"{\n",
" LANG: C\n",
"} \"Enforce detgammabar = detgammahat (= 1 in Cartesian)\"\n",
"\n",
"schedule BaikalETK_BoundaryConditions_evolved_gfs in MoL_PostStep\n",
"{\n",
" LANG: C\n",
" OPTIONS: LEVEL\n",
" SYNC: evol_variables\n",
"} \"Apply boundary conditions and perform AMR+interprocessor synchronization\"\n",
"\n",
"schedule GROUP ApplyBCs as BaikalETK_ApplyBCs in MoL_PostStep after BaikalETK_BoundaryConditions_evolved_gfs\n",
"{\n",
"} \"Group for applying boundary conditions\"\n",
"\n",
"\n",
"# Next update ADM quantities\n",
"\n",
"schedule BaikalETK_BSSN_to_ADM in MoL_PostStep after BaikalETK_ApplyBCs before ADMBase_SetADMVars\n",
"{\n",
" LANG: C\n",
" OPTIONS: Local\n",
"} \"Perform BSSN-to-ADM conversion. Useful for diagnostics.\"\n",
"\n",
"# Compute Hamiltonian & momentum constraints\n",
"\"\"\")\n",
" if add_stress_energy_source_terms == True:\n",
" file.write(\"\"\"\n",
"schedule driver_BSSN_T4UU in MoL_PseudoEvolution before BaikalETK_BSSN_constraints\n",
"{\n",
" LANG: C\n",
" OPTIONS: Local\n",
"} \"MoL_PseudoEvolution: Compute T4UU, needed for BSSN constraints\"\n",
"\"\"\")\n",
" file.write(\"\"\"\n",
"\n",
"schedule BaikalETK_BSSN_constraints in MoL_PseudoEvolution\n",
"{\n",
" LANG: C\n",
" OPTIONS: Local\n",
"} \"Compute BSSN (Hamiltonian and momentum) constraints\"\n",
"\n",
"schedule BaikalETK_BoundaryConditions_aux_gfs in MoL_PseudoEvolution after BaikalETK_BSSN_constraints\n",
"{\n",
" LANG: C\n",
" OPTIONS: LOCAL # Needed so that cctk_nghostzones[0] (the number of boundary points) is defined.\n",
" # In other words, don't use LEVEL mode here, or the number of boundary points\n",
" # filled may not match the actual number of ghost zones. Weird, huh?\n",
" SYNC: aux_variables\n",
"} \"Enforce symmetry BCs in constraint computation\"\n",
"\n",
"\"\"\")\n",
" if add_stress_energy_source_terms == True:\n",
" file.write(\"\"\"\n",
"schedule BaikalETK_BSSN_to_ADM in MoL_PseudoEvolution after BaikalETK_BoundaryConditions_aux_gfs\n",
"{\n",
" LANG: C\n",
" OPTIONS: Local\n",
"} \"Perform BSSN-to-ADM conversion in MoL_PseudoEvolution. Needed for proper HydroBase integration.\"\n",
"\"\"\")\n",
" file.write(\"\"\"\n",
"schedule GROUP ApplyBCs as BaikalETK_auxgfs_ApplyBCs in MoL_PseudoEvolution after BaikalETK_BoundaryConditions_aux_gfs\n",
"{\n",
"} \"Apply boundary conditions\"\n",
"\"\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='cdrivers'></a>\n",
"\n",
"# Step 4: C driver functions for ETK registration & NRPy+-generated kernels \\[Back to [top](#toc)\\]\n",
"$$\\label{cdrivers}$$\n",
"\n",
"Now that we have constructed the basic C code kernels and the needed Einstein Toolkit `ccl` files, we next write the driver functions for registering `BaikalETK` within the Toolkit and the C code kernels. Each of these driver functions is called directly from [`schedule.ccl`](#scheduleccl).\n",
"\n",
"<a id='etkfunctions'></a>\n",
"## Step 4.a: Needed ETK functions: Banner, Symmetry registration, Parameter sanity check, Method of Lines (`MoL`) registration, Boundary condition \\[Back to [top](#toc)\\]\n",
"$$\\label{etkfunctions}$$\n",
"\n",
"### To-do: Parameter sanity check function. E.g., error should be thrown if `cctk_nghostzones[]` is set too small for the chosen finite-differencing order within NRPy+."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"make_code_defn_list = []\n",
"def append_to_make_code_defn_list(filename):\n",
" if filename not in make_code_defn_list:\n",
" make_code_defn_list.append(filename)\n",
" return os.path.join(outdir,filename)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"with open(append_to_make_code_defn_list(\"RegisterSlicing.c\"),\"w\") as file:\n",
" file.write(\"\"\" \n",
"#include \"cctk.h\"\n",
"\n",
"#include \"Slicing.h\"\n",
"\n",
"int BaikalETK_RegisterSlicing (void)\n",
"{\n",
" Einstein_RegisterSlicing (\"BaikalETK\");\n",
" return 0;\n",
"}\"\"\")"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# First the ETK banner code, proudly showing the NRPy+ banner\n",
"import NRPy_logo as logo\n",
"\n",
"with open(append_to_make_code_defn_list(\"Banner.c\"),\"w\") as file:\n",
" file.write(\"\"\"\n",
"#include <stdio.h>\n",
"\n",
"void BaikalETK_Banner() \n",
"{\n",
" \"\"\")\n",
" logostr = logo.print_logo(print_to_stdout=False)\n",
" file.write(\"printf(\\\"BaikalETK: another Einstein Toolkit thorn generated by\\\\n\\\");\\n\")\n",
" for line in logostr.splitlines():\n",
" file.write(\" printf(\\\"\"+line+\"\\\\n\\\");\\n\")\n",
" file.write(\"}\\n\")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# Next register symmetries\n",
"with open(append_to_make_code_defn_list(\"Symmetry_registration_oldCartGrid3D.c\"),\"w\") as file:\n",
" file.write(\"\"\"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"#include \"Symmetry.h\"\n",
"\n",
"void BaikalETK_Symmetry_registration(CCTK_ARGUMENTS)\n",
"{\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
" // Stores gridfunction parity across x=0, y=0, and z=0 planes, respectively\n",
" int sym[3];\n",
" \n",
" // Next register parities for each gridfunction based on its name \n",
" // (to ensure this algorithm is robust, gridfunctions with integers\n",
" // in their base names are forbidden in NRPy+).\n",
"\"\"\")\n",
" full_gfs_list = []\n",
" full_gfs_list.extend(evol_gfs_list)\n",
" full_gfs_list.extend(auxevol_gfs_list)\n",
" full_gfs_list.extend(aux_gfs_list)\n",
" for gf in full_gfs_list:\n",
" file.write(\"\"\"\n",
" // Default to scalar symmetry:\n",
" sym[0] = 1; sym[1] = 1; sym[2] = 1;\n",
" // Now modify sym[0], sym[1], and/or sym[2] as needed \n",
" // to account for gridfunction parity across \n",
" // x=0, y=0, and/or z=0 planes, respectively\n",
"\"\"\")\n",
" # If gridfunction name does not end in a digit, by NRPy+ syntax, it must be a scalar\n",
" if gf[len(gf) - 1].isdigit() == False:\n",
" pass # Scalar = default\n",
" elif len(gf) > 2:\n",
" # Rank-1 indexed expression (e.g., vector)\n",
" if gf[len(gf) - 2].isdigit() == False:\n",
" if int(gf[-1]) > 2:\n",
" print(\"Error: Found invalid gridfunction name: \"+gf)\n",
" sys.exit(1)\n",
" symidx = gf[-1]\n",
" file.write(\" sym[\"+symidx+\"] = -1;\\n\")\n",
" # Rank-2 indexed expression\n",
" elif gf[len(gf) - 2].isdigit() == True:\n",
" if len(gf) > 3 and gf[len(gf) - 3].isdigit() == True:\n",
" print(\"Error: Found a Rank-3 or above gridfunction: \"+gf+\", which is at the moment unsupported.\")\n",
" print(\"It should be easy to support this if desired.\")\n",
" sys.exit(1)\n",
" symidx0 = gf[-2]\n",
" file.write(\" sym[\"+symidx0+\"] *= -1;\\n\")\n",
" symidx1 = gf[-1]\n",
" file.write(\" sym[\"+symidx1+\"] *= -1;\\n\")\n",
" else:\n",
" print(\"Don't know how you got this far with a gridfunction named \"+gf+\", but I'll take no more of this nonsense.\")\n",
" print(\" Please follow best-practices and rename your gridfunction to be more descriptive\")\n",
" sys.exit(1)\n",
" file.write(\" SetCartSymVN(cctkGH, sym, \\\"BaikalETK::\"+gf+\"\\\");\\n\")\n",
" file.write(\"}\\n\")"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"# Next register symmetries\n",
"with open(append_to_make_code_defn_list(\"zero_rhss.c\"),\"w\") as file:\n",
" file.write(\"\"\"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"#include \"Symmetry.h\"\n",
"\n",
"void BaikalETK_zero_rhss(CCTK_ARGUMENTS)\n",
"{\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
"\"\"\")\n",
" set_rhss_to_zero = \"\"\n",
" for gf in rhs_list:\n",
" set_rhss_to_zero += gf+\"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)] = 0.0;\\n\"\n",
" \n",
" file.write(lp.loop([\"i2\",\"i1\",\"i0\"],[\"0\", \"0\", \"0\"],\n",
" [\"cctk_lsh[2]\",\"cctk_lsh[1]\",\"cctk_lsh[0]\"],\n",
" [\"1\",\"1\",\"1\"],\n",
" [\"#pragma omp parallel for\",\"\",\"\",],\"\",set_rhss_to_zero))\n",
" file.write(\"}\\n\")"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# Next registration with the Method of Lines thorn\n",
"with open(append_to_make_code_defn_list(\"MoL_registration.c\"),\"w\") as file:\n",
" file.write(\"\"\"\n",
"//--------------------------------------------------------------------------\n",
"// Register with the Method of Lines time stepper\n",
"// (MoL thorn, found in arrangements/CactusBase/MoL)\n",
"// MoL documentation located in arrangements/CactusBase/MoL/doc\n",
"//--------------------------------------------------------------------------\n",
"#include <stdio.h>\n",
"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"\n",
"#include \"Symmetry.h\"\n",
"\n",
"void BaikalETK_MoL_registration(CCTK_ARGUMENTS)\n",
"{\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
" CCTK_INT ierr = 0, group, rhs;\n",
"\n",
" // Register evolution & RHS gridfunction groups with MoL, so it knows\n",
"\n",
" group = CCTK_GroupIndex(\"BaikalETK::evol_variables\");\n",
" rhs = CCTK_GroupIndex(\"BaikalETK::evol_variables_rhs\");\n",
" ierr += MoLRegisterEvolvedGroup(group, rhs);\n",
" \n",
" if (ierr) CCTK_ERROR(\"Problems registering with MoL\");\n",
"}\n",
"\"\"\")"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"# Next register with the boundary conditions thorns.\n",
"# PART 1: Set BC type to \"none\" for all variables\n",
"# Since we choose NewRad boundary conditions, we must register all \n",
"# gridfunctions to have boundary type \"none\". This is because\n",
"# NewRad is seen by the rest of the Toolkit as a modification to the\n",
"# RHSs.\n",
"\n",
"# This code is based on Kranc's McLachlan/ML_BSSN/src/Boundaries.cc code.\n",
"with open(append_to_make_code_defn_list(\"BoundaryConditions.c\"),\"w\") as file:\n",
" file.write(\"\"\"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"#include \"cctk_Faces.h\"\n",
"#include \"util_Table.h\"\n",
"#include \"Symmetry.h\"\n",
"\n",
"// Set `none` boundary conditions on BSSN RHSs, as these are set via NewRad.\n",
"void BaikalETK_BoundaryConditions_evolved_gfs(CCTK_ARGUMENTS)\n",
"{\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
" CCTK_INT ierr CCTK_ATTRIBUTE_UNUSED = 0;\n",
"\"\"\")\n",
" for gf in evol_gfs_list:\n",
" file.write(\"\"\"\n",
" ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, \"BaikalETK::\"\"\"+gf+\"\"\"\", \"none\");\n",
" if (ierr < 0) CCTK_ERROR(\"Failed to register BC for BaikalETK::\"\"\"+gf+\"\"\"!\");\n",
"\"\"\")\n",
" file.write(\"\"\"\n",
"}\n",
"\n",
"// Set `flat` boundary conditions on BSSN constraints, similar to what Lean does.\n",
"void BaikalETK_BoundaryConditions_aux_gfs(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
" CCTK_INT ierr CCTK_ATTRIBUTE_UNUSED = 0;\n",
"\n",
"\"\"\")\n",
" for gf in aux_gfs_list:\n",
" file.write(\"\"\"\n",
" ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, cctk_nghostzones[0], -1, \"BaikalETK::\"\"\"+gf+\"\"\"\", \"flat\");\n",
" if (ierr < 0) CCTK_ERROR(\"Failed to register BC for BaikalETK::\"\"\"+gf+\"\"\"!\");\n",
"\"\"\")\n",
" file.write(\"}\\n\")\n",
" \n",
" \n",
"# PART 2: Set C code for calling NewRad BCs\n",
"# As explained in lean_public/LeanBSSNMoL/src/calc_bssn_rhs.F90,\n",
"# the function NewRad_Apply takes the following arguments:\n",
"# NewRad_Apply(cctkGH, var, rhs, var0, v0, radpower),\n",
"# which implement the boundary condition:\n",
"# var = var_at_infinite_r + u(r-var_char_speed*t)/r^var_radpower\n",
"# Obviously for var_radpower>0, var_at_infinite_r is the value of\n",
"# the variable at r->infinity. var_char_speed is the propagation\n",
"# speed at the outer boundary, and var_radpower is the radial \n",
"# falloff rate.\n",
"\n",
"with open(append_to_make_code_defn_list(\"BoundaryCondition_NewRad.c\"),\"w\") as file:\n",
" file.write(\"\"\"\n",
"#include <math.h>\n",
"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"\n",
"void BaikalETK_NewRad(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
"\"\"\")\n",
" for gf in evol_gfs_list:\n",
" var_at_infinite_r = \"0.0\"\n",
" var_char_speed = \"1.0\"\n",
" var_radpower = \"1.0\"\n",
" \n",
" if gf == \"alpha\":\n",
" var_at_infinite_r = \"1.0\"\n",
" if LapseCondition == \"OnePlusLog\":\n",
" var_char_speed = \"sqrt(2.0)\"\n",
" else:\n",
" pass # 1.0 (default) is fine\n",
" if \"aDD\" in gf or \"trK\" in gf: # consistent with Lean code.\n",
" var_radpower = \"2.0\"\n",
"\n",
" file.write(\" NewRad_Apply(cctkGH, \"+gf+\", \"+gf.replace(\"GF\",\"\")+\"_rhsGF, \"+var_at_infinite_r+\", \"+\n",
" var_char_speed+\", \"+var_radpower+\");\\n\")\n",
" file.write(\"}\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='bssnadmconversions'></a>\n",
"\n",
"## Step 4.b: BSSN $\\leftrightarrow$ ADM conversions \\[Back to [top](#toc)\\]\n",
"$$\\label{bssnadmconversions}$$\n",
"\n",
"<a id='admtobssn'></a>\n",
"\n",
"### Step 4.b.i: ADM $\\to$ BSSN \\[Back to [top](#toc)\\]\n",
"$$\\label{admtobssn}$$\n",
"\n",
"Initial data in the Einstein Toolkit are given in terms of [ADM quantities](https://en.wikipedia.org/wiki/ADM_formalism), so a conversion is necessary so the quantities are in terms of BSSN variables used for evolving the initial data forward in time."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"# First we convert from ADM to BSSN, as is required to convert initial data \n",
"# (given using) ADM quantities, to the BSSN evolved variables\n",
"import BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear as atob\n",
"IDhDD,IDaDD,IDtrK,IDvetU,IDbetU,IDalpha,IDcf,IDlambdaU = \\\n",
" atob.Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear(\"Cartesian\",\"DoNotOutputADMInputFunction\",outdir)\n",
"\n",
"alphaSphorCart = gri.register_gridfunctions( \"AUXEVOL\", \"alphaSphorCart\")\n",
"betaSphorCartU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\", \"betaSphorCartU\")\n",
"BSphorCartU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\", \"BSphorCartU\")\n",
"gammaSphorCartDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\", \"gammaSphorCartDD\", \"sym01\")\n",
"KSphorCartDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\", \"KSphorCartDD\", \"sym01\")\n",
"\n",
"# Step : Output ADM to BSSN conversion.\n",
"with open(append_to_make_code_defn_list(\"ADM_to_BSSN.c\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"#include <math.h>\n",
"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"\n",
"void BaikalETK_ADM_to_BSSN(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
" CCTK_REAL *alphaSphorCartGF = alp;\n",
"\"\"\")\n",
" # It's ugly if we output code in the following ordering, so we'll first\n",
" # output to a string and then sort the string to beautify the code a bit.\n",
" outstr = []\n",
" for i in range(DIM):\n",
" outstr.append(\" CCTK_REAL *betaSphorCartU\"+str(i)+\"GF = beta\"+chr(ord('x')+i)+\";\\n\")\n",
" outstr.append(\" CCTK_REAL *BSphorCartU\"+str(i)+\"GF = dtbeta\"+chr(ord('x')+i)+\";\\n\")\n",
" for j in range(i,DIM):\n",
" outstr.append(\" CCTK_REAL *gammaSphorCartDD\"+str(i)+str(j)+\"GF = g\"+chr(ord('x')+i)+chr(ord('x')+j)+\";\\n\")\n",
" outstr.append(\" CCTK_REAL *KSphorCartDD\"+str(i)+str(j)+\"GF = k\"+chr(ord('x')+i)+chr(ord('x')+j)+\";\\n\")\n",
" outstr.sort()\n",
" for line in outstr:\n",
" file.write(line)\n",
" file.write(\"\"\"\n",
" const CCTK_REAL invdx0 = 1.0/CCTK_DELTA_SPACE(0);\n",
" const CCTK_REAL invdx1 = 1.0/CCTK_DELTA_SPACE(1);\n",
" const CCTK_REAL invdx2 = 1.0/CCTK_DELTA_SPACE(2);\n",
"\"\"\")\n",
"\n",
" all_but_lambdaU_expressions = [\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD00\"),rhs=IDhDD[0][0]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD01\"),rhs=IDhDD[0][1]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD02\"),rhs=IDhDD[0][2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD11\"),rhs=IDhDD[1][1]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD12\"),rhs=IDhDD[1][2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD22\"),rhs=IDhDD[2][2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD00\"),rhs=IDaDD[0][0]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD01\"),rhs=IDaDD[0][1]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD02\"),rhs=IDaDD[0][2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD11\"),rhs=IDaDD[1][1]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD12\"),rhs=IDaDD[1][2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD22\"),rhs=IDaDD[2][2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"trK\"),rhs=IDtrK),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"vetU0\"),rhs=IDvetU[0]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"vetU1\"),rhs=IDvetU[1]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"vetU2\"),rhs=IDvetU[2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"betU0\"),rhs=IDbetU[0]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"betU1\"),rhs=IDbetU[1]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"betU2\"),rhs=IDbetU[2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"alpha\"),rhs=IDalpha),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"cf\"),rhs=IDcf)]\n",
" \n",
" outCparams = \"preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False\"\n",
" all_but_lambdaU_outC = fin.FD_outputC(\"returnstring\",all_but_lambdaU_expressions, outCparams)\n",
" file.write(lp.loop([\"i2\",\"i1\",\"i0\"],[\"0\",\"0\",\"0\"],[\"cctk_lsh[2]\",\"cctk_lsh[1]\",\"cctk_lsh[0]\"],\n",
" [\"1\",\"1\",\"1\"],[\"#pragma omp parallel for\",\"\",\"\"],\"\",all_but_lambdaU_outC))\n",
"\n",
" outCparams = \"preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False\"\n",
" lambdaU_expressions = [lhrh(lhs=gri.gfaccess(\"in_gfs\",\"lambdaU0\"),rhs=IDlambdaU[0]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"lambdaU1\"),rhs=IDlambdaU[1]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"lambdaU2\"),rhs=IDlambdaU[2])]\n",
" lambdaU_expressions_FDout = fin.FD_outputC(\"returnstring\",lambdaU_expressions, outCparams)\n",
"\n",
" file.write(lp.loop([\"i2\",\"i1\",\"i0\"],[\"cctk_nghostzones[2]\",\"cctk_nghostzones[1]\",\"cctk_nghostzones[0]\"],\n",
" [\"cctk_lsh[2]-cctk_nghostzones[2]\",\"cctk_lsh[1]-cctk_nghostzones[1]\",\"cctk_lsh[0]-cctk_nghostzones[0]\"],\n",
" [\"1\",\"1\",\"1\"],[\"#pragma omp parallel for\",\"\",\"\"],\"\",lambdaU_expressions_FDout))\n",
" \n",
" file.write(\"\"\"\n",
" ExtrapolateGammas(cctkGH,lambdaU0GF);\n",
" ExtrapolateGammas(cctkGH,lambdaU1GF);\n",
" ExtrapolateGammas(cctkGH,lambdaU2GF);\n",
"}\n",
"\"\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='bssntoadm'></a>\n",
"\n",
"### Step 4.b.ii: BSSN $\\to$ ADM \\[Back to [top](#toc)\\]\n",
"$$\\label{bssntoadm}$$\n",
"\n",
"All modules (thorns) in the Einstein Toolkit that deal with spacetime quantities do so via the core `ADMBase` module, which assumes variables are written in ADM form. Therefore, in order for `BaikalETK` to interface properly with the rest of the Toolkit, its native BSSN variables must be converted to ADM quantities."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"import BSSN.ADM_in_terms_of_BSSN as btoa\n",
"btoa.ADM_in_terms_of_BSSN()\n",
"Bq.BSSN_basic_tensors() # Gives us betaU & BU\n",
"\n",
"with open(append_to_make_code_defn_list(\"BSSN_to_ADM.c\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"#include <math.h>\n",
"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"\n",
"void BaikalETK_BSSN_to_ADM(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
" \n",
"\"\"\")\n",
" btoa_lhrh = []\n",
" for i in range(DIM):\n",
" for j in range(i,DIM):\n",
" btoa_lhrh.append(lhrh(lhs=\"g\"+chr(ord('x')+i)+chr(ord('x')+j)+\"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)]\",\n",
" rhs=btoa.gammaDD[i][j]))\n",
" for i in range(DIM):\n",
" for j in range(i,DIM):\n",
" btoa_lhrh.append(lhrh(lhs=\"k\"+chr(ord('x')+i)+chr(ord('x')+j)+\"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)]\",\n",
" rhs=btoa.KDD[i][j]))\n",
" btoa_lhrh.append(lhrh(lhs=\"alp[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)]\",rhs=Bq.alpha))\n",
" \n",
" for i in range(DIM):\n",
" btoa_lhrh.append(lhrh(lhs=\"beta\"+chr(ord('x')+i)+\"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)]\",\n",
" rhs=Bq.betaU[i]))\n",
"\n",
" for i in range(DIM):\n",
" btoa_lhrh.append(lhrh(lhs=\"dtbeta\"+chr(ord('x')+i)+\"[CCTK_GFINDEX3D(cctkGH,i0,i1,i2)]\",\n",
" rhs=Bq.BU[i]))\n",
" \n",
" outCparams = \"preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False\"\n",
" bssn_to_adm_Ccode = fin.FD_outputC(\"returnstring\",btoa_lhrh, outCparams)\n",
" file.write(lp.loop([\"i2\",\"i1\",\"i0\"],[\"0\",\"0\",\"0\"],[\"cctk_lsh[2]\",\"cctk_lsh[1]\",\"cctk_lsh[0]\"],\n",
" [\"1\",\"1\",\"1\"],[\"#pragma omp parallel for\",\"\",\"\"],\"\",bssn_to_adm_Ccode))\n",
"\n",
" file.write(\"}\\n\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='bssnrhss'></a>\n",
"\n",
"## Step 4.c: Evaluate BSSN right-hand-sides (RHSs) \\[Back to [top](#toc)\\]\n",
"$$\\label{bssnrhss}$$\n",
"\n",
"<a id='ricci'></a>\n",
"\n",
"### Step 4.c.i: Evaluate Ricci tensor \\[Back to [top](#toc)\\]\n",
"$$\\label{ricci}$$\n",
"\n",
"To slightly optimize the performance of `BaikalETK`'s BSSN solver, we split the computation of the [complicated expressions for the Ricci tensor $\\\\bar{R}_{ij}$](Tutorial-BSSN_quantities.ipynb#rbar) into its own function, and then use the result when evaluating the BSSN right-hand-side (RHS) expressions."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"with open(append_to_make_code_defn_list(\"driver_pt1_BSSN_Ricci.c\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"#include <math.h>\n",
"\n",
"#include \"SIMD/SIMD_intrinsics.h\"\n",
"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"\n",
"void driver_pt1_BSSN_Ricci(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" \n",
" const CCTK_REAL NOSIMDinvdx0 = 1.0/CCTK_DELTA_SPACE(0);\n",
" const REAL_SIMD_ARRAY invdx0 = ConstSIMD(NOSIMDinvdx0);\n",
" const CCTK_REAL NOSIMDinvdx1 = 1.0/CCTK_DELTA_SPACE(1);\n",
" const REAL_SIMD_ARRAY invdx1 = ConstSIMD(NOSIMDinvdx1);\n",
" const CCTK_REAL NOSIMDinvdx2 = 1.0/CCTK_DELTA_SPACE(2);\n",
" const REAL_SIMD_ARRAY invdx2 = ConstSIMD(NOSIMDinvdx2);\n",
"#include \"BSSN_Ricci.h\"\n",
"}\\n\"\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='bssnrhssricciinput'></a>\n",
"\n",
"### Step 4.c.ii: Evaluate BSSN RHSs, using Ricci tensor as input \\[Back to [top](#toc)\\]\n",
"$$\\label{bssnrhssricciinput}$$\n",
"\n",
"Next we construct the driver function for evaluating the BSSN RHSs, which make use of the Ricci tensor $\\bar{R}_{ij}$, which has just been computed."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"def SIMD_declare_C_params():\n",
" SIMD_declare_C_params_str = \"\"\n",
" for i in range(len(par.glb_Cparams_list)):\n",
" # keep_param is a boolean indicating whether we should accept or reject\n",
" # the parameter. singleparstring will contain the string indicating\n",
" # the variable type.\n",
" keep_param, singleparstring = keep_param__return_type(par.glb_Cparams_list[i])\n",
"\n",
" if (keep_param) and (\"CCTK_REAL\" in singleparstring):\n",
" parname = par.glb_Cparams_list[i].parname\n",
" SIMD_declare_C_params_str += \" const \"+singleparstring + \"*NOSIMD\"+parname+\\\n",
" \" = CCTK_ParameterGet(\\\"\"+parname+\"\\\",\\\"BaikalETK\\\",NULL);\\n\"\n",
" SIMD_declare_C_params_str += \" const REAL_SIMD_ARRAY \"+parname+\" = ConstSIMD(*NOSIMD\"+parname+\");\\n\"\n",
" return SIMD_declare_C_params_str"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"with open(append_to_make_code_defn_list(\"driver_pt2_BSSN_RHSs.c\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"#include <math.h>\n",
"\n",
"#include \"SIMD/SIMD_intrinsics.h\"\n",
"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"\n",
"//void BSSN_RHSs()\n",
"\n",
"void driver_pt2_BSSN_RHSs(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
"\n",
" const CCTK_REAL NOSIMDinvdx0 = 1.0/CCTK_DELTA_SPACE(0);\n",
" const REAL_SIMD_ARRAY invdx0 = ConstSIMD(NOSIMDinvdx0);\n",
" const CCTK_REAL NOSIMDinvdx1 = 1.0/CCTK_DELTA_SPACE(1);\n",
" const REAL_SIMD_ARRAY invdx1 = ConstSIMD(NOSIMDinvdx1);\n",
" const CCTK_REAL NOSIMDinvdx2 = 1.0/CCTK_DELTA_SPACE(2);\n",
" const REAL_SIMD_ARRAY invdx2 = ConstSIMD(NOSIMDinvdx2);\n",
"\"\"\"+SIMD_declare_C_params()+\"\"\"\n",
"#include \"BSSN_RHSs.h\"\n",
"}\\n\"\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='enforcegammahatconstraint'></a>\n",
"\n",
"## Step 4.d: Enforcing conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint (in Cartesian coordinates, $\\det{\\hat{\\gamma}_{ij}}=1$) \\[Back to [top](#toc)\\]\n",
"$$\\label{enforcegammahatconstraint}$$\n",
"\n",
"Here we construct the driver function for enforcing the conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint. The BSSN equations are not strongly hyperbolic if this condition is not set."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"with open(append_to_make_code_defn_list(\"enforce_detgammabar_constraint.c\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"#include <math.h>\n",
"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"\n",
"void enforce_detgammabar_constraint(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
"\n",
"#include \"enforce_detgammabar_constraint.h\"\n",
"}\\n\"\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='diagnostics'></a>\n",
"\n",
"## Step 4.e: Diagnostics: Computing Hamiltonian & momentum constraints \\[Back to [top](#toc)\\]\n",
"$$\\label{diagnostics}$$\n",
"\n",
"The BSSN Hamiltonian & momentum constraints are useful diagnostics of a numerical-relativity calculation's health, as both should converge to zero with increasing numerical resolution. Here we construct the driver function."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"with open(append_to_make_code_defn_list(\"BSSN_constraints.c\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"#include <math.h>\n",
"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"\n",
"void BaikalETK_BSSN_constraints(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
"\n",
" const CCTK_REAL invdx0 = 1.0/CCTK_DELTA_SPACE(0);\n",
" const CCTK_REAL invdx1 = 1.0/CCTK_DELTA_SPACE(1);\n",
" const CCTK_REAL invdx2 = 1.0/CCTK_DELTA_SPACE(2);\n",
"\n",
"#include \"BSSN_constraints.h\"\n",
"}\\n\"\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='t4uu'></a>\n",
"\n",
"## Step 4.f: `driver_BSSN_T4UU()`: Compute $T^{\\mu\\nu}$ from `TmunuBase`'s $T_{\\mu\\nu}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{t4uu}$$\n",
"\n",
"Here we implement $T^{\\mu\\nu} = g^{\\mu \\delta} g^{\\nu \\gamma} T_{\\delta \\gamma}.$"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"if add_stress_energy_source_terms == True:\n",
" # Declare T4DD as a set of gridfunctions. These won't \n",
" # actually appear in interface.ccl, as interface.ccl \n",
" # was set above. Thus before calling the code output\n",
" # by FD_outputC(), we'll have to set pointers\n",
" # to the actual gridfunctions they reference.\n",
" # (In this case the eTab's.)\n",
" T4DD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"T4DD\",\"sym01\",DIM=4)\n",
" import BSSN.ADMBSSN_tofrom_4metric as AB4m\n",
" AB4m.g4UU_ito_BSSN_or_ADM(\"BSSN\")\n",
"\n",
" T4UUraised = ixp.zerorank2(DIM=4)\n",
" for mu in range(4):\n",
" for nu in range(4):\n",
" for delta in range(4):\n",
" for gamma in range(4):\n",
" T4UUraised[mu][nu] += AB4m.g4UU[mu][delta]*AB4m.g4UU[nu][gamma]*T4DD[delta][gamma]\n",
" \n",
" T4UU_expressions = [\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"T4UU00\"),rhs=T4UUraised[0][0]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"T4UU01\"),rhs=T4UUraised[0][1]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"T4UU02\"),rhs=T4UUraised[0][2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"T4UU03\"),rhs=T4UUraised[0][3]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"T4UU11\"),rhs=T4UUraised[1][1]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"T4UU12\"),rhs=T4UUraised[1][2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"T4UU13\"),rhs=T4UUraised[1][3]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"T4UU22\"),rhs=T4UUraised[2][2]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"T4UU23\"),rhs=T4UUraised[2][3]),\n",
" lhrh(lhs=gri.gfaccess(\"in_gfs\",\"T4UU33\"),rhs=T4UUraised[3][3])]\n",
"\n",
" outCparams = \"outCverbose=False,includebraces=False,preindent=2,SIMD_enable=True\"\n",
" T4UUstr = fin.FD_outputC(\"returnstring\",T4UU_expressions, outCparams)\n",
" T4UUstr_loop = lp.loop([\"i2\",\"i1\",\"i0\"],[\"0\",\"0\",\"0\"],[\"cctk_lsh[2]\",\"cctk_lsh[1]\",\"cctk_lsh[0]\"],\n",
" [\"1\",\"1\",\"SIMD_width\"],[\"#pragma omp parallel for\",\"\",\"\"],\"\",T4UUstr)\n",
"\n",
" with open(append_to_make_code_defn_list(\"driver_BSSN_T4UU.c\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"#include <math.h>\n",
"\n",
"#include \"cctk.h\"\n",
"#include \"cctk_Arguments.h\"\n",
"#include \"cctk_Parameters.h\"\n",
"\n",
"#include \"SIMD/SIMD_intrinsics.h\"\n",
"\n",
"void driver_BSSN_T4UU(CCTK_ARGUMENTS) {\n",
" DECLARE_CCTK_ARGUMENTS;\n",
" DECLARE_CCTK_PARAMETERS;\n",
"\n",
" const CCTK_REAL *restrict T4DD00GF = eTtt;\n",
" const CCTK_REAL *restrict T4DD01GF = eTtx;\n",
" const CCTK_REAL *restrict T4DD02GF = eTty;\n",
" const CCTK_REAL *restrict T4DD03GF = eTtz;\n",
" const CCTK_REAL *restrict T4DD11GF = eTxx;\n",
" const CCTK_REAL *restrict T4DD12GF = eTxy;\n",
" const CCTK_REAL *restrict T4DD13GF = eTxz;\n",
" const CCTK_REAL *restrict T4DD22GF = eTyy;\n",
" const CCTK_REAL *restrict T4DD23GF = eTyz;\n",
" const CCTK_REAL *restrict T4DD33GF = eTzz;\n",
"\"\"\"+T4UUstr_loop+\"\"\"\n",
"}\\n\"\"\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='makecodedefn'></a>\n",
"\n",
"## Step 4.g: `make.code.defn`: List of all C driver functions needed to compile `BaikalETK` \\[Back to [top](#toc)\\]\n",
"$$\\label{makecodedefn}$$\n",
"\n",
"When constructing each C code driver function above, we called the `append_to_make_code_defn_list()` function, which built a list of each C code driver file. We'll now add each of those files to the `make.code.defn` file, used by the Einstein Toolkit's build system."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"with open(os.path.join(outdir,\"make.code.defn\"), \"w\") as file:\n",
" file.write(\"\"\"\n",
"# Main make.code.defn file for thorn BaikalETK\n",
"\n",
"# Source files in this directory\n",
"SRCS =\"\"\")\n",
" filestring = \"\"\n",
" for i in range(len(make_code_defn_list)):\n",
" filestring += \" \"+make_code_defn_list[i]\n",
" if i != len(make_code_defn_list)-1:\n",
" filestring += \" \\\\\\n\"\n",
" else:\n",
" filestring += \"\\n\"\n",
" file.write(filestring)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validation'></a>\n",
"\n",
"# Step 5: Code validation \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validation}$$\n",
"\n",
"Here we will show plots demonstrating good agreement between `BaikalETK` and, e.g., `McLachlan` (another, trusted ETK thorn)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='self_validation'></a>\n",
"\n",
"## Step 5.a: Validation against [BaikalETK.BaikalETK_Pymodule](../edit/BaikalETK/BaikalETK_Pymodule.py) module \\[Back to [top](#toc)\\]\n",
"$$\\label{self_validation}$$\n",
"\n",
"As a self-validation check, we verify agreement in all codes generated by\n",
"1. this tutorial notebook, and \n",
"2. the NRPy+ [Baikal.BaikalETK_Pymodule](../edit/BaikalETK/BaikalETK_Pymodule.py) module."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Generating symbolic expressions for BSSN RHSs...\n",
"Finished BSSN symbolic expressions in 0.96791672706604 seconds.\n",
"Generating C code for BSSN RHSs in Cartesian coordinates.\n",
"Finished BSSN_RHS C codegen in 16.42121911048889 seconds.\n",
"Generating C code for Ricci tensor in Cartesian coordinates.\n",
"Finished Ricci C codegen in 19.169066429138184 seconds.\n",
"Generating optimized C code for Ham. & mom. constraints. May take a while, depending on CoordSystem.\n",
"Finished Hamiltonian & momentum constraint C codegen in 6.465627431869507 seconds.\n",
"Generating optimized C code for gamma constraint. May take a while, depending on CoordSystem.\n",
"Finished gamma constraint C codegen in 0.04613375663757324 seconds.\n"
]
}
],
"source": [
"# First generate all codes using BaikalETK.BaikalETK_Pymodule.BaikalETK_codegen():\n",
"gri.glb_gridfcs_list = []\n",
"\n",
"import BaikalETK.BaikalETK_Pymodule as BE\n",
"outvalrootdir = \"BaikalETK-validate\"\n",
"BE.BaikalETK_codegen(outrootdir=outvalrootdir,\n",
" FD_order=FD_order, # Finite difference order: even numbers only, starting with 2. 12 is generally unstable\n",
" LapseCondition = LapseCondition,\n",
" ShiftCondition = ShiftCondition, # Set the standard, second-order advecting-shift,\n",
" # Gamma-driving shift condition\n",
" add_stress_energy_source_terms = add_stress_energy_source_terms, # Enable stress-energy terms?\n",
" default_KO_strength = default_KO_strength)"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Ignoring lines with only whitespace:\n",
"PASSED: interface.ccl\n",
"PASSED: param.ccl\n",
"PASSED: schedule.ccl\n",
"PASSED: driver_pt2_BSSN_RHSs.c\n",
"PASSED: BoundaryCondition_NewRad.c\n",
"PASSED: BSSN_constraints.h\n",
"PASSED: BoundaryConditions.c\n",
"PASSED: BSSN_RHSs.h\n",
"PASSED: MoL_registration.c\n",
"PASSED: driver_pt1_BSSN_Ricci.c\n",
"PASSED: BSSN_constraints.c\n",
"PASSED: BSSN_to_ADM.c\n",
"PASSED: RegisterSlicing.c\n",
"PASSED: ADM_to_BSSN.c\n",
"PASSED: BSSN_Ricci.h\n",
"PASSED: make.code.defn\n",
"PASSED: Symmetry_registration_oldCartGrid3D.c\n",
"PASSED: enforce_detgammabar_constraint.h\n",
"PASSED: enforce_detgammabar_constraint.c\n",
"PASSED: Banner.c\n",
"PASSED: zero_rhss.c\n"
]
}
],
"source": [
"# Then compare all files generated by this notebook & the separate Python module.\n",
"import difflib\n",
"rootfiles = [f for f in os.listdir(outrootdir) if os.path.isfile(os.path.join(outrootdir, f)) and \"ccl\" in os.path.join(outrootdir, f)]\n",
"srcfiles = [f for f in os.listdir(outdir) if os.path.isfile(os.path.join(outdir, f))]\n",
"\n",
"def compare_two_files(outdir1,outdir2,filebase):\n",
" with open(os.path.join(outdir1,file)) as file1, open(os.path.join(outdir2,file)) as file2:\n",
" # Read the lines of each file\n",
" file1_lines = file1.readlines()\n",
" file2_lines = file2.readlines()\n",
" num_diffs = 0\n",
" file1_lines_noleadingwhitespace = []\n",
" for line in file1_lines:\n",
" if line.strip() == \"\": # If the line contains only whitespace, remove all leading whitespace\n",
" file1_lines_noleadingwhitespace.append(line.lstrip())\n",
" else:\n",
" file1_lines_noleadingwhitespace.append(line)\n",
" file2_lines_noleadingwhitespace = []\n",
" for line in file2_lines:\n",
" if line.strip() == \"\": # If the line contains only whitespace, remove all leading whitespace\n",
" file2_lines_noleadingwhitespace.append(line.lstrip())\n",
" else:\n",
" file2_lines_noleadingwhitespace.append(line)\n",
" for line in difflib.unified_diff(file1_lines_noleadingwhitespace, file2_lines_noleadingwhitespace, \n",
" fromfile=os.path.join(outdir1,file), \n",
" tofile =os.path.join(outdir2,file)):\n",
" sys.stdout.writelines(line)\n",
" num_diffs = num_diffs + 1\n",
" if num_diffs == 0:\n",
" print(\"PASSED: \"+file)\n",
" else:\n",
" print(\"FAILED (see diff above): \"+file)\n",
"\n",
"print(\"Ignoring lines with only whitespace:\")\n",
"for file in rootfiles:\n",
" compare_two_files(outrootdir,outvalrootdir,file)\n",
"for file in srcfiles:\n",
" compare_two_files(outdir,os.path.join(outvalrootdir,\"src\"),file)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='latex_pdf_output'></a>\n",
"\n",
"# Step 6: 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-BaikalETK.pdf](Tutorial-BaikalETK.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": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)\r\n",
" restricted \\write18 enabled.\r\n",
"entering extended mode\r\n",
"This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)\r\n",
" restricted \\write18 enabled.\r\n",
"entering extended mode\r\n",
"This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)\r\n",
" restricted \\write18 enabled.\r\n",
"entering extended mode\r\n"
]
}
],
"source": [
"!jupyter nbconvert --to latex --template latex_nrpy_style.tplx --log-level='WARN' Tutorial-BaikalETK.ipynb\n",
"!pdflatex -interaction=batchmode Tutorial-BaikalETK.tex\n",
"!pdflatex -interaction=batchmode Tutorial-BaikalETK.tex\n",
"!pdflatex -interaction=batchmode Tutorial-BaikalETK.tex\n",
"!rm -f Tut*.out Tut*.aux Tut*.log"
]
}
],
"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.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}