Tutorial-BSSN_time_evolution-C_codegen_library.ipynb

{
 "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",
    "# BSSN Time-Evolution C Code Generation Library\n",
    "\n",
    "## Author: Zach Etienne\n",
    "\n",
    "## This module implements a number of helper functions for generating C-code kernels that solve Einstein's equations in the covariant BSSN formalism [described in this NRPy+ tutorial notebook](Tutorial-BSSN_formulation.ipynb). The kernel generation process includes the creation of the BSSN RHS expressions for the Method of Lines time integration, evaluation of BSSN constraints as a measure of numerical error, computation of the 3-Ricci tensor $\\mathcal{R}{ij}$, and the enforcement of the conformal 3-metric constraint $det \\bar{\\gamma}{ij} = det \\hat{\\gamma}_{ij}$. Additional functionality includes the generation of Weyl scalar $\\psi_4$ related to gravitational wave strain, and its tetrad.\n",
    "\n",
    "**Notebook Status:** <font color = red><b> Not yet validated </b></font>\n",
    "\n",
    "**Validation Notes:** This module has NOT been validated to exhibit convergence to zero of the Hamiltonian constraint violation at the expected order to the exact solution *after a short numerical evolution of the initial data* (see [plots at bottom](#convergence)), and all quantities have been validated against the [original SENR code](https://bitbucket.org/zach_etienne/nrpy).\n",
    "\n",
    "### NRPy+ modules that generate needed symbolic expressions:\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\n",
    "* [BSSN/Enforce_Detgammahat_Constraint.py](../edit/BSSN/Enforce_Detgammahat_Constraint.py); [**tutorial**](Tutorial-BSSN_enforcing_determinant_gammabar_equals_gammahat_constraint.ipynb): Generates symbolic expressions for enforcing the $\\det{\\bar{\\gamma}}=\\det{\\hat{\\gamma}}$ constraint\n",
    "\n",
    "## Introduction:\n",
    "Here we use NRPy+ to generate the C source code kernels necessary to generate C functions needed/useful for evolving forward in time the BSSN equations, including:\n",
    "1. the BSSN RHS expressions for [Method of Lines](https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html) time integration, with arbitrary gauge choice.\n",
    "1. the BSSN constraints as a check of numerical error"
   ]
  },
  {
   "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](#importmodules): Import needed Python modules\n",
    "1. [Step 2](#helperfuncs): Helper Python functions for C code generation\n",
    "1. [Step 3.a](#bssnrhs): Generate symbolic BSSN RHS expressions\n",
    "1. [Step 3.b](#bssnrhs_c_code): `rhs_eval()`: Register C function for evaluating BSSN RHS expressions\n",
    "1. [Step 3.c](#ricci): Generate symbolic expressions for 3-Ricci tensor $\\bar{R}_{ij}$\n",
    "1. [Step 3.d](#ricci_c_code): `Ricci_eval()`: Register C function for evaluating 3-Ricci tensor $\\bar{R}_{ij}$\n",
    "1. [Step 4.a](#bssnconstraints): Generate symbolic expressions for BSSN Hamiltonian & momentum constraints\n",
    "1. [Step 4.b](#bssnconstraints_c_code): `BSSN_constraints()`: Register C function for evaluating BSSN Hamiltonian & momentum constraints\n",
    "1. [Step 5](#enforce3metric): `enforce_detgammahat_constraint()`: Register C function for enforcing the conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint\n",
    "1. [Step 6.a](#psi4): `psi4_part_{0,1,2}()`: Register C function for evaluating Weyl scalar $\\psi_4$, in 3 parts (3 functions)\n",
    "1. [Step 6.b](#psi4_tetrad): `psi4_tetrad()`: Register C function for evaluating Weyl scalar $\\psi_4$ tetrad\n",
    "1. [Step 6.c](#swm2): `SpinWeight_minus2_SphHarmonics()`: Register C function for evaluating spin-weight $s=-2$ spherical harmonics\n",
    "1. [Step 7](#validation): Confirm above functions are bytecode-identical to those in `BSSN/BSSN_Ccodegen_library.py`\n",
    "1. [Step 8](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='importmodules'></a>\n",
    "\n",
    "# Step 1: Import needed Python modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{importmodules}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:28:22.699247Z",
     "iopub.status.busy": "2021-03-07T17:28:22.691709Z",
     "iopub.status.idle": "2021-03-07T17:28:23.477052Z",
     "shell.execute_reply": "2021-03-07T17:28:23.476205Z"
    }
   },
   "outputs": [],
   "source": [
    "# RULES FOR ADDING FUNCTIONS TO THIS ROUTINE:\n",
    "# 1. The function must be runnable from a multiprocessing environment,\n",
    "#    which means that the function\n",
    "# 1.a: cannot depend on previous function calls.\n",
    "# 1.b: cannot create directories (this is not multiproc friendly)\n",
    "\n",
    "\n",
    "# Step P1: Import needed NRPy+ core modules:\n",
    "from __future__ import print_function\n",
    "from outputC import lhrh, add_to_Cfunction_dict  # NRPy+: Core C code output module\n",
    "import finite_difference as fin  # NRPy+: Finite difference C code generation module\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import grid as gri               # NRPy+: Functions having to do with numerical grids\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "from pickling import pickle_NRPy_env   # NRPy+: Pickle/unpickle NRPy+ environment, for parallel codegen\n",
    "import os, time             # Standard Python modules for multiplatform OS-level functions, benchmarking\n",
    "import BSSN.BSSN_RHSs as rhs\n",
    "import BSSN.BSSN_gauge_RHSs as gaugerhs\n",
    "import loop as lp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='helperfuncs'></a>\n",
    "\n",
    "# Step 2: Helper Python functions for C code generation \\[Back to [top](#toc)\\]\n",
    "$$\\label{helperfuncs}$$\n",
    "\n",
    "* `print_msg_with_timing()` gives the user an idea of what's going on/taking so long. Also outputs timing info.\n",
    "* `get_loopopts()` sets up options for NRPy+'s `loop` module\n",
    "* `register_stress_energy_source_terms_return_T4UU()` registers gridfunctions for $T^{\\mu\\nu}$ if needed and not yet registered."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:28:22.699247Z",
     "iopub.status.busy": "2021-03-07T17:28:22.691709Z",
     "iopub.status.idle": "2021-03-07T17:28:23.477052Z",
     "shell.execute_reply": "2021-03-07T17:28:23.476205Z"
    }
   },
   "outputs": [],
   "source": [
    "###############################################\n",
    "# Helper Python functions for C code generation\n",
    "# print_msg_with_timing() gives the user an idea of what's going on/taking so long. Also outputs timing info.\n",
    "def print_msg_with_timing(desc, msg=\"Symbolic\", startstop=\"start\", starttime=0.0):\n",
    "    CoordSystem = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
    "    elapsed = time.time()-starttime\n",
    "    if msg == \"Symbolic\":\n",
    "        if startstop == \"start\":\n",
    "            print(\"Generating symbolic expressions for \" + desc + \" (\" + CoordSystem + \" coords)...\")\n",
    "            return time.time()\n",
    "        # else:\n",
    "        print(\"Finished generating symbolic expressions for \"+desc+\n",
    "              \" (\" + CoordSystem + \" coords) in \"+str(round(elapsed, 1))+\" seconds. Next up: C codegen...\")\n",
    "    elif msg == \"Ccodegen\":\n",
    "        if startstop == \"start\":\n",
    "            print(\"Generating C code for \"+desc+\" (\" + CoordSystem + \" coords)...\")\n",
    "            return time.time()\n",
    "        # else:\n",
    "        print(\"Finished generating C code for \"+desc+\" (\" + CoordSystem + \" coords) in \"+str(round(elapsed, 1))+\" seconds.\")\n",
    "\n",
    "\n",
    "# get_loopopts() sets up options for NRPy+'s loop module\n",
    "def get_loopopts(points_to_update, enable_SIMD, enable_rfm_precompute, OMP_pragma_on, enable_xxs=True):\n",
    "    loopopts = points_to_update + \",includebraces=False\"\n",
    "    if enable_SIMD:\n",
    "        loopopts += \",enable_SIMD\"\n",
    "    if enable_rfm_precompute:\n",
    "        loopopts += \",enable_rfm_precompute\"\n",
    "    elif not enable_xxs:\n",
    "        pass\n",
    "    else:\n",
    "        loopopts += \",Read_xxs\"\n",
    "    if OMP_pragma_on != \"i2\":\n",
    "        loopopts += \",pragma_on_\"+OMP_pragma_on\n",
    "    return loopopts\n",
    "\n",
    "\n",
    "# register_stress_energy_source_terms_return_T4UU() registers gridfunctions\n",
    "#        for T4UU if needed and not yet registered.\n",
    "def register_stress_energy_source_terms_return_T4UU(enable_stress_energy_source_terms):\n",
    "    if enable_stress_energy_source_terms:\n",
    "        registered_already = False\n",
    "        for i in range(len(gri.glb_gridfcs_list)):\n",
    "            if gri.glb_gridfcs_list[i].name == \"T4UU00\":\n",
    "                registered_already = True\n",
    "        if not registered_already:\n",
    "            return ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\", \"T4UU\", \"sym01\", DIM=4)\n",
    "        # else:\n",
    "        return ixp.declarerank2(\"T4UU\", \"sym01\", DIM=4)\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssnrhs'></a>\n",
    "\n",
    "# Step 3.a: Generate symbolic BSSN RHS expressions \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssnrhs}$$\n",
    "\n",
    "First, we generate the symbolic expressions. Be sure to call this function from within a `reference_metric::enable_rfm_precompute=\"True\"` environment if reference metric precomputation is desired.\n",
    "\n",
    "\n",
    "`BSSN_RHSs__generate_symbolic_expressions()` supports the following features\n",
    "\n",
    "* (`\"OnePlusLog\"` by default) Lapse gauge choice\n",
    "* (`\"GammaDriving2ndOrder_Covariant\"` by default) Shift gauge choice\n",
    "* (disabled by default) Kreiss-Oliger dissipation\n",
    "* (disabled by default) Stress-energy ($T^{\\mu\\nu}$) source terms\n",
    "* (enabled by default) \"Leave Ricci symbolic\": do not compute the 3-Ricci tensor $\\bar{R}_{ij}$ within the BSSN RHSs, which only adds to the extreme complexity of the BSSN RHS expressions. Instead, leave computation of $\\bar{R}_{ij}$=`RbarDD` to a separate function. Doing this generally increases C-code performance by about 10%.\n",
    "\n",
    "Two lists are returned by this function:\n",
    "\n",
    "1. `betaU`: the un-rescaled shift vector $\\beta^i$, which is used to perform upwinding.\n",
    "1. `BSSN_RHSs_SymbExpressions`: the BSSN RHS symbolic expressions, using the `lhrh` named-tuple to store a list of LHSs and RHSs, where each LHS and RHS is defined as follows\n",
    "    1. LHS = BSSN gridfunction whose time derivative is being computed at grid point `i0,i1,i2`, and \n",
    "    1. RHS = time derivative expression for a given variable at the given point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BSSN_RHSs__generate_symbolic_expressions(LapseCondition=\"OnePlusLog\",\n",
    "                                             ShiftCondition=\"GammaDriving2ndOrder_Covariant\",\n",
    "                                             enable_KreissOliger_dissipation=True,\n",
    "                                             enable_stress_energy_source_terms=False,\n",
    "                                             leave_Ricci_symbolic=True):\n",
    "    ######################################\n",
    "    # START: GENERATE SYMBOLIC EXPRESSIONS\n",
    "    starttime = print_msg_with_timing(\"BSSN_RHSs\", msg=\"Symbolic\", startstop=\"start\")\n",
    "\n",
    "    # Returns None if enable_stress_energy_source_terms==False; otherwise returns symb expressions for T4UU\n",
    "    T4UU = register_stress_energy_source_terms_return_T4UU(enable_stress_energy_source_terms)\n",
    "\n",
    "    # Evaluate BSSN RHSs:\n",
    "    import BSSN.BSSN_quantities as Bq\n",
    "    par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\", str(leave_Ricci_symbolic))\n",
    "    rhs.BSSN_RHSs()\n",
    "\n",
    "    if enable_stress_energy_source_terms:\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(3):\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(3):\n",
    "                rhs.a_rhsDD[i][j] += Bsest.sourceterm_a_rhsDD[i][j]\n",
    "\n",
    "    par.set_parval_from_str(\"BSSN.BSSN_gauge_RHSs::LapseEvolutionOption\", LapseCondition)\n",
    "    par.set_parval_from_str(\"BSSN.BSSN_gauge_RHSs::ShiftEvolutionOption\", ShiftCondition)\n",
    "    gaugerhs.BSSN_gauge_RHSs()  # Can depend on above RHSs\n",
    "    # Restore BSSN.BSSN_quantities::LeaveRicciSymbolic to False\n",
    "    par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\", \"False\")\n",
    "\n",
    "    # Add Kreiss-Oliger dissipation to the BSSN RHSs:\n",
    "    if enable_KreissOliger_dissipation:\n",
    "        thismodule = \"KO_Dissipation\"\n",
    "        diss_strength = par.Cparameters(\"REAL\", thismodule, \"diss_strength\", 0.1) # *Bq.cf # *Bq.cf*Bq.cf*Bq.cf # cf**1 is found better than cf**4 over the long term.\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(3):\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(3):\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(3):\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",
    "    # END: GENERATE SYMBOLIC EXPRESSIONS\n",
    "    ######################################\n",
    "\n",
    "    lhs_names = [\"alpha\", \"cf\", \"trK\"]\n",
    "    rhs_exprs = [gaugerhs.alpha_rhs, rhs.cf_rhs, rhs.trK_rhs]\n",
    "    for i in range(3):\n",
    "        lhs_names.append(\"betU\" + str(i))\n",
    "        rhs_exprs.append(gaugerhs.bet_rhsU[i])\n",
    "        lhs_names.append(\"lambdaU\" + str(i))\n",
    "        rhs_exprs.append(rhs.lambda_rhsU[i])\n",
    "        lhs_names.append(\"vetU\" + str(i))\n",
    "        rhs_exprs.append(gaugerhs.vet_rhsU[i])\n",
    "        for j in range(i, 3):\n",
    "            lhs_names.append(\"aDD\" + str(i) + str(j))\n",
    "            rhs_exprs.append(rhs.a_rhsDD[i][j])\n",
    "            lhs_names.append(\"hDD\" + str(i) + str(j))\n",
    "            rhs_exprs.append(rhs.h_rhsDD[i][j])\n",
    "\n",
    "    # Sort the lhss list alphabetically, and rhss to match.\n",
    "    #   This ensures the RHSs are evaluated in the same order\n",
    "    #   they're allocated in memory:\n",
    "    lhs_names, rhs_exprs = [list(x) for x in zip(*sorted(zip(lhs_names, rhs_exprs), key=lambda pair: pair[0]))]\n",
    "\n",
    "    # Declare the list of lhrh's\n",
    "    BSSN_RHSs_SymbExpressions = []\n",
    "    for var in range(len(lhs_names)):\n",
    "        BSSN_RHSs_SymbExpressions.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", lhs_names[var]), rhs=rhs_exprs[var]))\n",
    "\n",
    "    print_msg_with_timing(\"BSSN_RHSs\", msg=\"Symbolic\", startstop=\"stop\", starttime=starttime)\n",
    "    return [betaU, BSSN_RHSs_SymbExpressions]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssnrhs_c_code'></a>\n",
    "\n",
    "# Step 3.b: `rhs_eval()`: Register C code for BSSN RHS expressions \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssnrhs_c_code}$$\n",
    "\n",
    "`add_rhs_eval_to_Cfunction_dict()` supports the following features\n",
    "\n",
    "* (enabled by default) reference-metric precomputation\n",
    "* (disabled by default) \"golden kernels\", which greatly increases the C-code generation time in an attempt to reduce computational cost. Most often this results in no speed-up.\n",
    "* (enabled by default) SIMD output\n",
    "* (disabled by default) splitting of RHSs into smaller pieces (multiple loops) to improve performance. Doesn't help much.\n",
    "* (`\"OnePlusLog\"` by default) Lapse gauge choice\n",
    "* (`\"GammaDriving2ndOrder_Covariant\"` by default) Shift gauge choice\n",
    "* (disabled by default) enable Kreiss-Oliger dissipation\n",
    "* (disabled by default) add stress-energy ($T^{\\mu\\nu}$) source terms\n",
    "* (enabled by default) \"Leave Ricci symbolic\": do not compute the 3-Ricci tensor $\\bar{R}_{ij}$ within the BSSN RHSs, which only adds to the extreme complexity of the BSSN RHS expressions. Instead, leave computation of $\\bar{R}_{ij}$=`RbarDD` to a separate function. Doing this generally increases C-code performance by about 10%.\n",
    "* (`\"i2\"` by default) OpenMP pragma acts on which loop (assumes `i2` is outermost and `i0` is innermost loop). For axisymmetric or near-axisymmetric calculations, `\"i1\"` may be *significantly* faster.\n",
    "\n",
    "Also to enable parallel C-code kernel generation, the NRPy+ environment is pickled and returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_rhs_eval_to_Cfunction_dict(includes=None, rel_path_to_Cparams=os.path.join(\".\"),\n",
    "                                   enable_rfm_precompute=True, enable_golden_kernels=False,\n",
    "                                   enable_SIMD=True, enable_split_for_optimizations_doesnt_help=False,\n",
    "                                   LapseCondition=\"OnePlusLog\", ShiftCondition=\"GammaDriving2ndOrder_Covariant\",\n",
    "                                   enable_KreissOliger_dissipation=False, enable_stress_energy_source_terms=False,\n",
    "                                   leave_Ricci_symbolic=True, OMP_pragma_on=\"i2\",\n",
    "                                   func_name_suffix=\"\"):\n",
    "    if includes is None:\n",
    "        includes = []\n",
    "    if enable_SIMD:\n",
    "        includes += [os.path.join(\"SIMD\", \"SIMD_intrinsics.h\")]\n",
    "    enable_FD_functions = bool(par.parval_from_str(\"finite_difference::enable_FD_functions\"))\n",
    "    if enable_FD_functions:\n",
    "        includes += [\"finite_difference_functions.h\"]\n",
    "\n",
    "    # Set up the C function for the BSSN RHSs\n",
    "    desc = \"Evaluate the BSSN RHSs\"\n",
    "    func_name = \"rhs_eval\" + func_name_suffix\n",
    "    params = \"const paramstruct *restrict params, \"\n",
    "    if enable_rfm_precompute:\n",
    "        params += \"const rfm_struct *restrict rfmstruct, \"\n",
    "    else:\n",
    "        params += \"REAL *restrict xx[3], \"\n",
    "    params += \"\"\"\n",
    "              const REAL *restrict auxevol_gfs,const REAL *restrict in_gfs,REAL *restrict rhs_gfs\"\"\"\n",
    "\n",
    "    betaU, BSSN_RHSs_SymbExpressions = \\\n",
    "        BSSN_RHSs__generate_symbolic_expressions(LapseCondition=LapseCondition, ShiftCondition=ShiftCondition,\n",
    "                                                 enable_KreissOliger_dissipation=enable_KreissOliger_dissipation,\n",
    "                                                 enable_stress_energy_source_terms=enable_stress_energy_source_terms,\n",
    "                                                 leave_Ricci_symbolic=leave_Ricci_symbolic)\n",
    "\n",
    "    # Construct body:\n",
    "    preloop=\"\"\n",
    "    enableCparameters=True\n",
    "    # Set up preloop in case we're outputting code for the Einstein Toolkit (ETK)\n",
    "    if par.parval_from_str(\"grid::GridFuncMemAccess\") == \"ETK\":\n",
    "        params, preloop = set_ETK_func_params_preloop(func_name)\n",
    "        enableCparameters=False\n",
    "\n",
    "    FD_outCparams = \"outCverbose=False,enable_SIMD=\" + str(enable_SIMD)\n",
    "    FD_outCparams += \",GoldenKernelsEnable=\" + str(enable_golden_kernels)\n",
    "\n",
    "    loopopts = get_loopopts(\"InteriorPoints\", enable_SIMD, enable_rfm_precompute, OMP_pragma_on)\n",
    "    FDorder = par.parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\")\n",
    "    starttime = print_msg_with_timing(\"BSSN_RHSs (FD order=\"+str(FDorder)+\")\", msg=\"Ccodegen\", startstop=\"start\")\n",
    "    if enable_split_for_optimizations_doesnt_help and FDorder == 6:\n",
    "        loopopts += \",DisableOpenMP\"\n",
    "        BSSN_RHSs_SymbExpressions_pt1 = []\n",
    "        BSSN_RHSs_SymbExpressions_pt2 = []\n",
    "        for lhsrhs in BSSN_RHSs_SymbExpressions:\n",
    "            if \"BETU\" in lhsrhs.lhs or \"LAMBDAU\" in lhsrhs.lhs:\n",
    "                BSSN_RHSs_SymbExpressions_pt1.append(lhrh(lhs=lhsrhs.lhs, rhs=lhsrhs.rhs))\n",
    "            else:\n",
    "                BSSN_RHSs_SymbExpressions_pt2.append(lhrh(lhs=lhsrhs.lhs, rhs=lhsrhs.rhs))\n",
    "        preloop += \"\"\"#pragma omp parallel\n",
    "    {\n",
    "\"\"\"\n",
    "        preloopbody = fin.FD_outputC(\"returnstring\", BSSN_RHSs_SymbExpressions_pt1,\n",
    "                                     params=FD_outCparams,\n",
    "                                     upwindcontrolvec=betaU)\n",
    "        preloop += \"\\n#pragma omp for\\n\" + lp.simple_loop(loopopts, preloopbody)\n",
    "        preloop += \"\\n#pragma omp for\\n\"\n",
    "        body = fin.FD_outputC(\"returnstring\", BSSN_RHSs_SymbExpressions_pt2,\n",
    "                              params=FD_outCparams,\n",
    "                              upwindcontrolvec=betaU)\n",
    "        postloop = \"\\n    } // END #pragma omp parallel\\n\"\n",
    "    else:\n",
    "        preloop += \"\"\n",
    "        body = fin.FD_outputC(\"returnstring\", BSSN_RHSs_SymbExpressions,\n",
    "                              params=FD_outCparams,\n",
    "                              upwindcontrolvec=betaU)\n",
    "        postloop = \"\"\n",
    "    print_msg_with_timing(\"BSSN_RHSs (FD order=\"+str(FDorder)+\")\", msg=\"Ccodegen\", startstop=\"stop\", starttime=starttime)\n",
    "\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        name=func_name, params=params,\n",
    "        preloop=preloop, body=body, loopopts=loopopts, postloop=postloop,\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams, enableCparameters=enableCparameters)\n",
    "    return pickle_NRPy_env()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ricci'></a>\n",
    "\n",
    "# Step 3.c: Generate symbolic expressions for 3-Ricci tensor $\\bar{R}_{ij}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{ricci}$$\n",
    "\n",
    "As described above, we find a roughly 10% speedup by computing the 3-Ricci tensor $\\bar{R}_{ij}$ separately from the BSSN RHS equations and storing the 6 independent components in memory. Here we construct the symbolic expressions for all 6 independent components of $\\bar{R}_{ij}$ (which is symmetric under interchange of indices).\n",
    "\n",
    "`Ricci__generate_symbolic_expressions()` does not support any input parameters.\n",
    "\n",
    "One list is returned by `Ricci__generate_symbolic_expressions()`: `Ricci_SymbExpressions`, which contains a list of expressions for the six independent components of $\\bar{R}_{ij}$, using the `lhrh` named-tuple to store a list of LHSs and RHSs, where each LHS and RHS is defined as follows\n",
    "\n",
    "1. LHS = gridfunction representation of the component of $\\bar{R}_{ij}$, computed at grid point i0,i1,i2, and\n",
    "1. RHS = expression for given component of $\\bar{R}_{ij}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Ricci__generate_symbolic_expressions():\n",
    "    ######################################\n",
    "    # START: GENERATE SYMBOLIC EXPRESSIONS\n",
    "    starttime = print_msg_with_timing(\"3-Ricci tensor\", msg=\"Symbolic\", startstop=\"start\")\n",
    "\n",
    "    # Evaluate 3-Ricci tensor:\n",
    "    import BSSN.BSSN_quantities as Bq\n",
    "    par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\", \"False\")\n",
    "\n",
    "    # Register all BSSN gridfunctions if not registered already\n",
    "    Bq.BSSN_basic_tensors()\n",
    "    # Next compute Ricci tensor\n",
    "    Bq.RicciBar__gammabarDD_dHatD__DGammaUDD__DGammaU()\n",
    "    # END: GENERATE SYMBOLIC EXPRESSIONS\n",
    "    ######################################\n",
    "    # Must register RbarDD as gridfunctions, as we're outputting them to gridfunctions here:\n",
    "    foundit = False\n",
    "    for i in range(len(gri.glb_gridfcs_list)):\n",
    "        if \"RbarDD00\" in gri.glb_gridfcs_list[i].name:\n",
    "            foundit = True\n",
    "    if not foundit:\n",
    "        ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\", \"RbarDD\", \"sym01\")\n",
    "\n",
    "    Ricci_SymbExpressions = [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",
    "    print_msg_with_timing(\"3-Ricci tensor\", msg=\"Symbolic\", startstop=\"stop\", starttime=starttime)\n",
    "    return Ricci_SymbExpressions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ricci_c_code'></a>\n",
    "\n",
    "# Step 3.d: `Ricci_eval()`: Register C function for evaluating 3-Ricci tensor $\\bar{R}_{ij}$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{ricci_c_code}$$\n",
    "\n",
    "`add_Ricci_eval_to_Cfunction_dict()` supports the following features\n",
    "\n",
    "* (enabled by default) reference-metric precomputation\n",
    "* (disabled by default) \"golden kernels\", which greatly increases the C-code generation time in an attempt to reduce computational cost. Most often this results in no speed-up.\n",
    "* (enabled by default) SIMD output\n",
    "* (disabled by default) splitting of RHSs into smaller pieces (multiple loops) to improve performance. Doesn't help much.\n",
    "* (`\"i2\"` by default) OpenMP pragma acts on which loop (assumes `i2` is outermost and `i0` is innermost loop). For axisymmetric or near-axisymmetric calculations, `\"i1\"` may be *significantly* faster.\n",
    "\n",
    "Also to enable parallel C-code kernel generation, the NRPy+ environment is pickled and returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_Ricci_eval_to_Cfunction_dict(includes=None, rel_path_to_Cparams=os.path.join(\".\"),\n",
    "                                     enable_rfm_precompute=True, enable_golden_kernels=False, enable_SIMD=True,\n",
    "                                     enable_split_for_optimizations_doesnt_help=False, OMP_pragma_on=\"i2\",\n",
    "                                     func_name_suffix=\"\"):\n",
    "    if includes is None:\n",
    "        includes = []\n",
    "    if enable_SIMD:\n",
    "        includes += [os.path.join(\"SIMD\", \"SIMD_intrinsics.h\")]\n",
    "    enable_FD_functions = bool(par.parval_from_str(\"finite_difference::enable_FD_functions\"))\n",
    "    if enable_FD_functions:\n",
    "        includes += [\"finite_difference_functions.h\"]\n",
    "\n",
    "    # Set up the C function for the 3-Ricci tensor\n",
    "    desc = \"Evaluate the 3-Ricci tensor\"\n",
    "    func_name = \"Ricci_eval\" + func_name_suffix\n",
    "    params = \"const paramstruct *restrict params, \"\n",
    "    if enable_rfm_precompute:\n",
    "        params += \"const rfm_struct *restrict rfmstruct, \"\n",
    "    else:\n",
    "        params += \"REAL *restrict xx[3], \"\n",
    "    params += \"const REAL *restrict in_gfs, REAL *restrict auxevol_gfs\"\n",
    "\n",
    "    # Construct body:\n",
    "    Ricci_SymbExpressions = Ricci__generate_symbolic_expressions()\n",
    "    FD_outCparams = \"outCverbose=False,enable_SIMD=\" + str(enable_SIMD)\n",
    "    FD_outCparams += \",GoldenKernelsEnable=\" + str(enable_golden_kernels)\n",
    "    loopopts = get_loopopts(\"InteriorPoints\", enable_SIMD, enable_rfm_precompute, OMP_pragma_on)\n",
    "\n",
    "    FDorder = par.parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\")\n",
    "    starttime = print_msg_with_timing(\"3-Ricci tensor (FD order=\"+str(FDorder)+\")\", msg=\"Ccodegen\", startstop=\"start\")\n",
    "\n",
    "    # Construct body:\n",
    "    preloop=\"\"\n",
    "    enableCparameters=True\n",
    "    # Set up preloop in case we're outputting code for the Einstein Toolkit (ETK)\n",
    "    if par.parval_from_str(\"grid::GridFuncMemAccess\") == \"ETK\":\n",
    "        params, preloop = set_ETK_func_params_preloop(func_name)\n",
    "        enableCparameters=False\n",
    "\n",
    "    if enable_split_for_optimizations_doesnt_help and FDorder >= 8:\n",
    "        loopopts += \",DisableOpenMP\"\n",
    "        Ricci_SymbExpressions_pt1 = []\n",
    "        Ricci_SymbExpressions_pt2 = []\n",
    "        for lhsrhs in Ricci_SymbExpressions:\n",
    "            if \"RBARDD00\" in lhsrhs.lhs or \"RBARDD11\" in lhsrhs.lhs or \"RBARDD22\" in lhsrhs.lhs:\n",
    "                Ricci_SymbExpressions_pt1.append(lhrh(lhs=lhsrhs.lhs, rhs=lhsrhs.rhs))\n",
    "            else:\n",
    "                Ricci_SymbExpressions_pt2.append(lhrh(lhs=lhsrhs.lhs, rhs=lhsrhs.rhs))\n",
    "        preloop = \"\"\"#pragma omp parallel\n",
    "    {\n",
    "#pragma omp for\n",
    "\"\"\"\n",
    "        preloopbody = fin.FD_outputC(\"returnstring\", Ricci_SymbExpressions_pt1,\n",
    "                                     params=FD_outCparams)\n",
    "        preloop += lp.simple_loop(loopopts, preloopbody)\n",
    "        preloop += \"#pragma omp for\\n\"\n",
    "        body = fin.FD_outputC(\"returnstring\", Ricci_SymbExpressions_pt2,\n",
    "                              params=FD_outCparams)\n",
    "        postloop = \"\\n    } // END #pragma omp parallel\\n\"\n",
    "    else:\n",
    "        body = fin.FD_outputC(\"returnstring\", Ricci_SymbExpressions,\n",
    "                              params=FD_outCparams)\n",
    "        postloop = \"\"\n",
    "    print_msg_with_timing(\"3-Ricci tensor (FD order=\"+str(FDorder)+\")\", msg=\"Ccodegen\", startstop=\"stop\", starttime=starttime)\n",
    "\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        name=func_name, params=params,\n",
    "        preloop=preloop, body=body, loopopts=loopopts, postloop=postloop,\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams, enableCparameters=enableCparameters)\n",
    "    return pickle_NRPy_env()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssnconstraints'></a>\n",
    "\n",
    "# Step 4.a: Generate symbolic expressions for BSSN Hamiltonian & momentum constraints \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssnconstraints}$$\n",
    "\n",
    "Next output the C code for evaluating the BSSN Hamiltonian and momentum constraints [(**Tutorial**)](Tutorial-BSSN_constraints.ipynb). In the absence of numerical error, these constraints should evaluate to zero. However, it does not due to numerical (typically truncation) error.\n",
    "\n",
    "We will therefore measure the constraint violations to gauge the accuracy of our simulation, and, ultimately determine whether errors are dominated by numerical finite differencing (truncation) error as expected.\n",
    "\n",
    "`BSSN_constraints__generate_symbolic_expressions()` supports the following features:\n",
    "\n",
    "* (disabled by default) add stress-energy ($T^{\\mu\\nu}$) source terms\n",
    "* (disabled by default) output Hamiltonian constraint only\n",
    "\n",
    "One list is returned by `BSSN_constraints__generate_symbolic_expressions()`: `BSSN_constraints_SymbExpressions`, which contains a list of expressions for the Hamiltonian and momentum constraints (4 elements total), using the `lhrh` named-tuple to store a list of LHSs and RHSs, where each LHS and RHS is defined as follows\n",
    "\n",
    "1. LHS = gridfunction representation of the BSSN constraint, computed at grid point i0,i1,i2, and\n",
    "1. RHS = expression for given BSSN constraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:28:28.150451Z",
     "iopub.status.busy": "2021-03-07T17:28:28.149610Z",
     "iopub.status.idle": "2021-03-07T17:28:28.151795Z",
     "shell.execute_reply": "2021-03-07T17:28:28.152300Z"
    }
   },
   "outputs": [],
   "source": [
    "def BSSN_constraints__generate_symbolic_expressions(enable_stress_energy_source_terms=False, leave_Ricci_symbolic=True,\n",
    "                                                    output_H_only=False):\n",
    "    ######################################\n",
    "    # START: GENERATE SYMBOLIC EXPRESSIONS\n",
    "    starttime = print_msg_with_timing(\"BSSN constraints\", msg=\"Symbolic\", startstop=\"start\")\n",
    "\n",
    "    # Define the Hamiltonian constraint and output the optimized C code.\n",
    "    par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\", str(leave_Ricci_symbolic))\n",
    "    import BSSN.BSSN_constraints as bssncon\n",
    "\n",
    "    # Returns None if enable_stress_energy_source_terms==False; otherwise returns symb expressions for T4UU\n",
    "    T4UU = register_stress_energy_source_terms_return_T4UU(enable_stress_energy_source_terms)\n",
    "\n",
    "    bssncon.BSSN_constraints(add_T4UUmunu_source_terms=False, output_H_only=output_H_only)  # We'll add them below if desired.\n",
    "    if enable_stress_energy_source_terms:\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(3):\n",
    "            bssncon.MU[i] += Bsest.sourceterm_MU[i]\n",
    "\n",
    "    BSSN_constraints_SymbExpressions = [lhrh(lhs=gri.gfaccess(\"aux_gfs\", \"H\"), rhs=bssncon.H)]\n",
    "    if not output_H_only:\n",
    "        BSSN_constraints_SymbExpressions += [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",
    "    par.set_parval_from_str(\"BSSN.BSSN_quantities::LeaveRicciSymbolic\", \"False\")\n",
    "    print_msg_with_timing(\"BSSN constraints\", msg=\"Symbolic\", startstop=\"stop\", starttime=starttime)\n",
    "    # END: GENERATE SYMBOLIC EXPRESSIONS\n",
    "    ######################################\n",
    "    return BSSN_constraints_SymbExpressions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssnconstraints_c_code'></a>\n",
    "\n",
    "# Step 4.b: `BSSN_constraints()`: Register C function for evaluating BSSN Hamiltonian & momentum constraints \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssnconstraints_c_code}$$\n",
    "\n",
    "`add_BSSN_constraints_to_Cfunction_dict()` supports the following features\n",
    "\n",
    "* (enabled by default) reference-metric precomputation\n",
    "* (disabled by default) \"golden kernels\", which greatly increases the C-code generation time in an attempt to reduce computational cost. Most often this results in no speed-up.\n",
    "* (enabled by default) SIMD output\n",
    "* (disabled by default) splitting of RHSs into smaller pieces (multiple loops) to improve performance. Doesn't help much.\n",
    "* (disabled by default) add stress-energy ($T^{\\mu\\nu}$) source terms\n",
    "* (disabled by default) output Hamiltonian constraint only\n",
    "* (`\"i2\"` by default) OpenMP pragma acts on which loop (assumes `i2` is outermost and `i0` is innermost loop). For axisymmetric or near-axisymmetric calculations, `\"i1\"` may be *significantly* faster.\n",
    "\n",
    "Also to enable parallel C-code kernel generation, the NRPy+ environment is pickled and returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_BSSN_constraints_to_Cfunction_dict(includes=None, rel_path_to_Cparams=os.path.join(\".\"),\n",
    "                                           enable_rfm_precompute=True, enable_golden_kernels=False, enable_SIMD=True,\n",
    "                                           enable_stress_energy_source_terms=False, leave_Ricci_symbolic=True,\n",
    "                                           output_H_only=False, OMP_pragma_on=\"i2\", func_name_suffix=\"\"):\n",
    "    if includes is None:\n",
    "        includes = []\n",
    "    if enable_SIMD:\n",
    "        includes += [os.path.join(\"SIMD\", \"SIMD_intrinsics.h\")]\n",
    "    enable_FD_functions = bool(par.parval_from_str(\"finite_difference::enable_FD_functions\"))\n",
    "    if enable_FD_functions:\n",
    "        includes += [\"finite_difference_functions.h\"]\n",
    "\n",
    "    # Set up the C function for the BSSN constraints\n",
    "    desc = \"Evaluate the BSSN constraints\"\n",
    "    func_name = \"BSSN_constraints\" + func_name_suffix\n",
    "    params = \"const paramstruct *restrict params, \"\n",
    "    if enable_rfm_precompute:\n",
    "        params += \"const rfm_struct *restrict rfmstruct, \"\n",
    "    else:\n",
    "        params += \"REAL *restrict xx[3], \"\n",
    "    params += \"\"\"\n",
    "                 const REAL *restrict in_gfs, const REAL *restrict auxevol_gfs, REAL *restrict aux_gfs\"\"\"\n",
    "\n",
    "    # Construct body:\n",
    "\n",
    "    BSSN_constraints_SymbExpressions = BSSN_constraints__generate_symbolic_expressions(enable_stress_energy_source_terms,\n",
    "                                                                                       leave_Ricci_symbolic=leave_Ricci_symbolic,\n",
    "                                                                                       output_H_only=output_H_only)\n",
    "\n",
    "    preloop=\"\"\n",
    "    enableCparameters=True\n",
    "    # Set up preloop in case we're outputting code for the Einstein Toolkit (ETK)\n",
    "    if par.parval_from_str(\"grid::GridFuncMemAccess\") == \"ETK\":\n",
    "        params, preloop = set_ETK_func_params_preloop(func_name)\n",
    "        enableCparameters=False\n",
    "\n",
    "    FD_outCparams = \"outCverbose=False,enable_SIMD=\" + str(enable_SIMD)\n",
    "    FD_outCparams += \",GoldenKernelsEnable=\" + str(enable_golden_kernels)\n",
    "    FDorder = par.parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\")\n",
    "    starttime = print_msg_with_timing(\"BSSN constraints (FD order=\"+str(FDorder)+\")\", msg=\"Ccodegen\", startstop=\"start\")\n",
    "    body = fin.FD_outputC(\"returnstring\", BSSN_constraints_SymbExpressions,\n",
    "                          params=FD_outCparams)\n",
    "    print_msg_with_timing(\"BSSN constraints (FD order=\"+str(FDorder)+\")\", msg=\"Ccodegen\", startstop=\"stop\", starttime=starttime)\n",
    "\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        name=func_name, params=params,\n",
    "        preloop=preloop,\n",
    "        body=body,\n",
    "        loopopts=get_loopopts(\"InteriorPoints\", enable_SIMD, enable_rfm_precompute, OMP_pragma_on),\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams, enableCparameters=enableCparameters)\n",
    "    return pickle_NRPy_env()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='enforce3metric'></a>\n",
    "\n",
    "# Step 5: `enforce_detgammahat_constraint()`: Register C function for enforcing the conformal 3-metric $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint \\[Back to [top](#toc)\\]\n",
    "$$\\label{enforce3metric}$$\n",
    "\n",
    "To ensure stability when solving the BSSN equations, we must enforce the 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). This function imposes the $\\det{\\bar{\\gamma}_{ij}}=\\det{\\hat{\\gamma}_{ij}}$ constraint.\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.\n",
    "\n",
    "`add_enforce_detgammahat_constraint_to_Cfunction_dict()` supports the following features\n",
    "\n",
    "* (enabled by default) reference-metric precomputation\n",
    "* (disabled by default) \"golden kernels\", which greatly increases the C-code generation time in an attempt to reduce computational cost. Most often this results in no speed-up.\n",
    "* (`\"i2\"` by default) OpenMP pragma acts on which loop (assumes `i2` is outermost and `i0` is innermost loop). For axisymmetric or near-axisymmetric calculations, `\"i1\"` may be *significantly* faster.\n",
    "\n",
    "Also to enable parallel C-code kernel generation, the NRPy+ environment is pickled and returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_enforce_detgammahat_constraint_to_Cfunction_dict(includes=None, rel_path_to_Cparams=os.path.join(\".\"),\n",
    "                                                         enable_rfm_precompute=True, enable_golden_kernels=False,\n",
    "                                                         OMP_pragma_on=\"i2\", func_name_suffix=\"\"):\n",
    "    # This function disables SIMD, as it includes cbrt() and abs() functions.\n",
    "    if includes is None:\n",
    "        includes = []\n",
    "    # This function does not use finite differences!\n",
    "    # enable_FD_functions = bool(par.parval_from_str(\"finite_difference::enable_FD_functions\"))\n",
    "    # if enable_FD_functions:\n",
    "    #     includes += [\"finite_difference_functions.h\"]\n",
    "\n",
    "    # Set up the C function for enforcing the det(gammabar) = det(gammahat) BSSN algebraic constraint\n",
    "    desc = \"Enforce the det(gammabar) = det(gammahat) (algebraic) constraint\"\n",
    "    func_name = \"enforce_detgammahat_constraint\" + func_name_suffix\n",
    "    params = \"const paramstruct *restrict params, \"\n",
    "    if enable_rfm_precompute:\n",
    "        params += \"const rfm_struct *restrict rfmstruct, \"\n",
    "    else:\n",
    "        params += \"REAL *restrict xx[3], \"\n",
    "    params += \"REAL *restrict in_gfs\"\n",
    "\n",
    "    # Construct body:\n",
    "    enforce_detg_constraint_symb_expressions = EGC.Enforce_Detgammahat_Constraint_symb_expressions()\n",
    "\n",
    "    preloop=\"\"\n",
    "    enableCparameters=True\n",
    "    # Set up preloop in case we're outputting code for the Einstein Toolkit (ETK)\n",
    "    if par.parval_from_str(\"grid::GridFuncMemAccess\") == \"ETK\":\n",
    "        params, preloop = set_ETK_func_params_preloop(func_name, enable_SIMD=False)\n",
    "        enableCparameters=False\n",
    "\n",
    "    FD_outCparams = \"outCverbose=False,enable_SIMD=False\"\n",
    "    FD_outCparams += \",GoldenKernelsEnable=\" + str(enable_golden_kernels)\n",
    "    starttime = print_msg_with_timing(\"Enforcing det(gammabar)=det(gammahat) constraint\", msg=\"Ccodegen\", startstop=\"start\")\n",
    "    body = fin.FD_outputC(\"returnstring\", enforce_detg_constraint_symb_expressions,\n",
    "                          params=FD_outCparams)\n",
    "    print_msg_with_timing(\"Enforcing det(gammabar)=det(gammahat) constraint\", msg=\"Ccodegen\", startstop=\"stop\", starttime=starttime)\n",
    "\n",
    "    enable_SIMD = False\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        name=func_name, params=params,\n",
    "        preloop=preloop,\n",
    "        body=body,\n",
    "        loopopts=get_loopopts(\"AllPoints\", enable_SIMD, enable_rfm_precompute, OMP_pragma_on),\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams, enableCparameters=enableCparameters)\n",
    "    return pickle_NRPy_env()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='psi4'></a>\n",
    "\n",
    "# Step 6.a: `psi4_part_{0,1,2}()`: Register C function for evaluating Weyl scalar $\\psi_4$, in 3 parts (3 functions) \\[Back to [top](#toc)\\]\n",
    "$$\\label{psi4}$$\n",
    "\n",
    "$\\psi_4$ is a complex scalar related to the gravitational wave strain via\n",
    "\n",
    "$$\n",
    "\\psi_4 = \\ddot{h}_+ - i \\ddot{h}_\\times.\n",
    "$$\n",
    "\n",
    "We construct the symbolic expression for $\\psi_4$ as described in the [corresponding NRPy+ Jupyter notebook](Tutorial-Psi4.ipynb), in three parts. The below `add_psi4_part_to_Cfunction_dict()` function will construct any of these three parts `0`, `1,` or `2`, and output the part to a function `psi4_part0()`, `psi4_part1()`, or `psi4_part2()`, respectively.\n",
    "\n",
    "`add_psi4_part_to_Cfunction_dict()` supports the following features\n",
    "\n",
    "* (`\"0\"` by default) which part? (`0`, `1,` or `2`), as described above\n",
    "* (disabled by default) \"setPsi4tozero\", which effectively turns this into a dummy function -- for when $\\psi_4$ is not needed, and it's easier to just set `psi_4=0` instead of calculating it.\n",
    "* (`\"i2\"` by default) OpenMP pragma acts on which loop (assumes `i2` is outermost and `i0` is innermost loop). For axisymmetric or near-axisymmetric calculations, `\"i1\"` may be *significantly* faster.\n",
    "\n",
    "Also to enable parallel C-code kernel generation, the NRPy+ environment is pickled and returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_psi4_part_to_Cfunction_dict(includes=None, rel_path_to_Cparams=os.path.join(\".\"), whichpart=0,\n",
    "                                    setPsi4tozero=False, OMP_pragma_on=\"i2\"):\n",
    "    starttime = print_msg_with_timing(\"psi4, part \" + str(whichpart), msg=\"Ccodegen\", startstop=\"start\")\n",
    "\n",
    "    # Set up the C function for psi4\n",
    "    if includes is None:\n",
    "        includes = []\n",
    "    includes += [\"NRPy_function_prototypes.h\"]\n",
    "\n",
    "    desc = \"Compute psi4 at all interior gridpoints, part \" + str(whichpart)\n",
    "    name = \"psi4_part\" + str(whichpart)\n",
    "    params = \"\"\"const paramstruct *restrict params, const REAL *restrict in_gfs, REAL *restrict xx[3], REAL *restrict aux_gfs\"\"\"\n",
    "\n",
    "    body = \"\"\n",
    "    gri.register_gridfunctions(\"AUX\", [\"psi4_part\" + str(whichpart) + \"re\", \"psi4_part\" + str(whichpart) + \"im\"])\n",
    "    FD_outCparams = \"outCverbose=False,enable_SIMD=False,CSE_sorting=none\"\n",
    "    if not setPsi4tozero:\n",
    "        # Set the body of the function\n",
    "        # First compute the symbolic expressions\n",
    "        psi4.Psi4(specify_tetrad=False)\n",
    "\n",
    "        # We really don't want to store these \"Cparameters\" permanently; they'll be set via function call...\n",
    "        #   so we make a copy of the original par.glb_Cparams_list (sans tetrad vectors) and restore it below\n",
    "        Cparams_list_orig = par.glb_Cparams_list.copy()\n",
    "        par.Cparameters(\"REAL\", __name__, [\"mre4U0\", \"mre4U1\", \"mre4U2\", \"mre4U3\"], [0, 0, 0, 0])\n",
    "        par.Cparameters(\"REAL\", __name__, [\"mim4U0\", \"mim4U1\", \"mim4U2\", \"mim4U3\"], [0, 0, 0, 0])\n",
    "        par.Cparameters(\"REAL\", __name__, [\"n4U0\", \"n4U1\", \"n4U2\", \"n4U3\"], [0, 0, 0, 0])\n",
    "\n",
    "        body += \"\"\"\n",
    "REAL mre4U0,mre4U1,mre4U2,mre4U3,mim4U0,mim4U1,mim4U2,mim4U3,n4U0,n4U1,n4U2,n4U3;\n",
    "psi4_tetrad(params,\n",
    "              in_gfs[IDX4S(CFGF, i0,i1,i2)],\n",
    "              in_gfs[IDX4S(HDD00GF, i0,i1,i2)],\n",
    "              in_gfs[IDX4S(HDD01GF, i0,i1,i2)],\n",
    "              in_gfs[IDX4S(HDD02GF, i0,i1,i2)],\n",
    "              in_gfs[IDX4S(HDD11GF, i0,i1,i2)],\n",
    "              in_gfs[IDX4S(HDD12GF, i0,i1,i2)],\n",
    "              in_gfs[IDX4S(HDD22GF, i0,i1,i2)],\n",
    "              &mre4U0,&mre4U1,&mre4U2,&mre4U3,&mim4U0,&mim4U1,&mim4U2,&mim4U3,&n4U0,&n4U1,&n4U2,&n4U3,\n",
    "              xx, i0,i1,i2);\n",
    "\"\"\"\n",
    "        body += \"REAL xCart_rel_to_globalgrid_center[3];\\n\"\n",
    "        body += \"xx_to_Cart(params, xx, i0, i1, i2,  xCart_rel_to_globalgrid_center);\\n\"\n",
    "        body += \"int ignore_Cart_to_i0i1i2[3];  REAL xx_rel_to_globalgridorigin[3];\\n\"\n",
    "        body += \"Cart_to_xx_and_nearest_i0i1i2_global_grid_center(params, xCart_rel_to_globalgrid_center,xx_rel_to_globalgridorigin,ignore_Cart_to_i0i1i2);\\n\"\n",
    "        for i in range(3):\n",
    "            body += \"const REAL xx\" + str(i) + \"=xx_rel_to_globalgridorigin[\" + str(i) + \"];\\n\"\n",
    "        body += fin.FD_outputC(\"returnstring\",\n",
    "                               [lhrh(lhs=gri.gfaccess(\"in_gfs\", \"psi4_part\" + str(whichpart) + \"re\"),\n",
    "                                     rhs=psi4.psi4_re_pt[whichpart]),\n",
    "                                lhrh(lhs=gri.gfaccess(\"in_gfs\", \"psi4_part\" + str(whichpart) + \"im\"),\n",
    "                                     rhs=psi4.psi4_im_pt[whichpart])],\n",
    "                               params=FD_outCparams)\n",
    "        par.glb_Cparams_list = Cparams_list_orig.copy()\n",
    "\n",
    "    elif setPsi4tozero:\n",
    "        body += fin.FD_outputC(\"returnstring\",\n",
    "                               [lhrh(lhs=gri.gfaccess(\"in_gfs\", \"psi4_part\" + str(whichpart) + \"re\"),\n",
    "                                     rhs=sp.sympify(0)),\n",
    "                                lhrh(lhs=gri.gfaccess(\"in_gfs\", \"psi4_part\" + str(whichpart) + \"im\"),\n",
    "                                     rhs=sp.sympify(0))],\n",
    "                               params=FD_outCparams)\n",
    "    enable_SIMD = False\n",
    "    enable_rfm_precompute = False\n",
    "\n",
    "    print_msg_with_timing(\"psi4, part \" + str(whichpart), msg=\"Ccodegen\", startstop=\"stop\", starttime=starttime)\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        name=name, params=params,\n",
    "        body=body,\n",
    "        loopopts=get_loopopts(\"InteriorPoints\", enable_SIMD, enable_rfm_precompute, OMP_pragma_on,\n",
    "                              enable_xxs=False),\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)\n",
    "    return pickle_NRPy_env()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='psi4_tetrad'></a>\n",
    "\n",
    "# Step 6.b: `psi4_tetrad()`: Register C function for evaluating Weyl scalar $\\psi_4$ tetrad \\[Back to [top](#toc)\\]\n",
    "$$\\label{psi4_tetrad}$$\n",
    "\n",
    "Computing $\\psi_4$ requires that an observer tetrad be specified. We adopt a \"quasi-Kinnersley tetrad\" as described in [the corresponding NRPy+ tutorial notebook](Tutorial-Psi4_tetrads.ipynb).\n",
    "\n",
    "`add_psi4_tetrad_to_Cfunction_dict()` supports the following features\n",
    "\n",
    "* (disabled by default) \"setPsi4tozero\", which effectively turns this into a dummy function -- for when $\\psi_4$ is not needed, and it's easier to just set `psi_4=0` instead of calculating it.\n",
    "\n",
    "Also to enable parallel C-code kernel generation, the NRPy+ environment is pickled and returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_psi4_tetrad_to_Cfunction_dict(includes=None, rel_path_to_Cparams=os.path.join(\".\"), setPsi4tozero=False):\n",
    "    starttime = print_msg_with_timing(\"psi4 tetrads\", msg=\"Ccodegen\", startstop=\"start\")\n",
    "\n",
    "    # Set up the C function for BSSN basis transformations\n",
    "    desc = \"Compute tetrad for psi4\"\n",
    "    name = \"psi4_tetrad\"\n",
    "\n",
    "    # First set up the symbolic expressions (RHSs) and their names (LHSs)\n",
    "    psi4tet.Psi4_tetrads()\n",
    "    list_of_varnames = []\n",
    "    list_of_symbvars = []\n",
    "    for i in range(4):\n",
    "        list_of_varnames.append(\"*mre4U\" + str(i))\n",
    "        list_of_symbvars.append(psi4tet.mre4U[i])\n",
    "    for i in range(4):\n",
    "        list_of_varnames.append(\"*mim4U\" + str(i))\n",
    "        list_of_symbvars.append(psi4tet.mim4U[i])\n",
    "    for i in range(4):\n",
    "        list_of_varnames.append(\"*n4U\" + str(i))\n",
    "        list_of_symbvars.append(psi4tet.n4U[i])\n",
    "\n",
    "    paramsindent = \"                  \"\n",
    "    params = \"\"\"const paramstruct *restrict params,\\n\"\"\" + paramsindent\n",
    "    list_of_metricvarnames = [\"cf\"]\n",
    "    for i in range(3):\n",
    "        for j in range(i, 3):\n",
    "            list_of_metricvarnames.append(\"hDD\" + str(i) + str(j))\n",
    "    for var in list_of_metricvarnames:\n",
    "        params += \"const REAL \" + var + \",\"\n",
    "    params += \"\\n\" + paramsindent\n",
    "    for var in list_of_varnames:\n",
    "        params += \"REAL \" + var + \",\"\n",
    "    params += \"\\n\" + paramsindent + \"REAL *restrict xx[3], const int i0,const int i1,const int i2\"\n",
    "\n",
    "    # Set the body of the function\n",
    "    body = \"\"\n",
    "    outCparams = \"includebraces=False,outCverbose=False,enable_SIMD=False,preindent=1\"\n",
    "    if not setPsi4tozero:\n",
    "        for i in range(3):\n",
    "            body += \"  const REAL xx\" + str(i) + \" = xx[\" + str(i) + \"][i\" + str(i) + \"];\\n\"\n",
    "        body += \"  // Compute tetrads:\\n\"\n",
    "        body += \"  {\\n\"\n",
    "        # Sort the lhss list alphabetically, and rhss to match:\n",
    "        lhss, rhss = [list(x) for x in zip(*sorted(zip(list_of_varnames, list_of_symbvars), key=lambda pair: pair[0]))]\n",
    "        body += outputC(rhss, lhss, filename=\"returnstring\", params=outCparams)\n",
    "        body += \"  }\\n\"\n",
    "\n",
    "    elif setPsi4tozero:\n",
    "        body += \"return;\\n\"\n",
    "    loopopts = \"\"\n",
    "\n",
    "    print_msg_with_timing(\"psi4 tetrads\", msg=\"Ccodegen\", startstop=\"stop\", starttime=starttime)\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        name=name, params=params,\n",
    "        body=body,\n",
    "        loopopts=loopopts,\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)\n",
    "    return pickle_NRPy_env()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='swm2'></a>\n",
    "\n",
    "# Step 6.c: `SpinWeight_minus2_SphHarmonics()`: Register C function for evaluating spin-weight $s=-2$ spherical harmonics \\[Back to [top](#toc)\\]\n",
    "$$\\label{swm2}$$\n",
    "\n",
    "After evaluating $\\psi_4$ at all interior gridpoints on a numerical grid, we next decompose $\\psi_4$ into spin-weight $s=-2$ spherical harmonics, which are documented in [this NRPy+ tutorial notebook](Tutorial-SpinWeighted_Spherical_Harmonics.ipynb).\n",
    "\n",
    "`SpinWeight_minus2_SphHarmonics()` supports the following features\n",
    "\n",
    "* (`\"8\"` by default) `maximum_l`, the maximum $\\ell$ mode to output. Symbolic expressions $(\\ell,m)$ modes up to and including `maximum_l` will be output.\n",
    "\n",
    "Also to enable parallel C-code kernel generation, the NRPy+ environment is pickled and returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_SpinWeight_minus2_SphHarmonics_to_Cfunction_dict(includes=None, rel_path_to_Cparams=os.path.join(\".\"),\n",
    "                                                         maximum_l=8):\n",
    "    starttime = print_msg_with_timing(\"Spin-weight s=-2 Spherical Harmonics\", msg=\"Ccodegen\", startstop=\"start\")\n",
    "\n",
    "    # Set up the C function for computing the spin-weight -2 spherical harmonic at theta,phi: Y_{s=-2, l,m}(theta,phi)\n",
    "    prefunc = r\"\"\"// Compute at a single point (th,ph) the spin-weight -2 spherical harmonic Y_{s=-2, l,m}(th,ph)\n",
    "// Manual \"inline void\" of this function results in compilation error with clang.\n",
    "void SpinWeight_minus2_SphHarmonics(const int l, const int m, const REAL th, const REAL ph,\n",
    "                                           REAL *reYlmswm2_l_m, REAL *imYlmswm2_l_m) {\n",
    "\"\"\"\n",
    "    # Construct prefunc:\n",
    "    outCparams = \"preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False\"\n",
    "    prefunc += \"\"\"\n",
    "  switch(l) {\n",
    "\"\"\"\n",
    "    for l in range(maximum_l + 1):  # Output values up to and including l=8.\n",
    "        prefunc += \"  case \" + str(l) + \":\\n\"\n",
    "        prefunc += \"    switch(m) {\\n\"\n",
    "        for m in range(-l, l + 1):\n",
    "            prefunc += \"    case \" + str(m) + \":\\n\"\n",
    "            prefunc += \"      {\\n\"\n",
    "            Y_m2_lm = SWm2SH.Y(-2, l, m, SWm2SH.th, SWm2SH.ph)\n",
    "            prefunc += outputC([sp.re(Y_m2_lm), sp.im(Y_m2_lm)], [\"*reYlmswm2_l_m\", \"*imYlmswm2_l_m\"],\n",
    "                            \"returnstring\", outCparams)\n",
    "            prefunc += \"      }\\n\"\n",
    "            prefunc += \"      return;\\n\"\n",
    "        prefunc += \"    }  // END switch(m)\\n\"\n",
    "    prefunc += \"  } // END switch(l)\\n\"\n",
    "\n",
    "    prefunc += r\"\"\"\n",
    "  fprintf(stderr, \"ERROR: SpinWeight_minus2_SphHarmonics handles only l=[0,\"\"\"+str(maximum_l)+r\"\"\"] and only m=[-l,+l] is defined.\\n\");\n",
    "  fprintf(stderr, \"       You chose l=%d and m=%d, which is out of these bounds.\\n\",l,m);\n",
    "  exit(1);\n",
    "}\n",
    "\n",
    "void lowlevel_decompose_psi4_into_swm2_modes(const int Nxx_plus_2NGHOSTS1,const int Nxx_plus_2NGHOSTS2,\n",
    "                                             const REAL dxx1, const REAL dxx2,\n",
    "                                             const REAL curr_time, const REAL R_ext,\n",
    "                                             const REAL *restrict th_array, const REAL *restrict sinth_array, const REAL *restrict ph_array,\n",
    "                                             const REAL *restrict psi4r_at_R_ext, const REAL *restrict psi4i_at_R_ext) {\n",
    "  for(int l=2;l<=\"\"\"+str(maximum_l)+r\"\"\";l++) {  // The maximum l here is set in Python.\n",
    "    for(int m=-l;m<=l;m++) {\n",
    "      // Parallelize the integration loop:\n",
    "      REAL psi4r_l_m = 0.0;\n",
    "      REAL psi4i_l_m = 0.0;\n",
    "#pragma omp parallel for reduction(+:psi4r_l_m,psi4i_l_m)\n",
    "      for(int i1=0;i1<Nxx_plus_2NGHOSTS1-2*NGHOSTS;i1++) {\n",
    "        const REAL th    = th_array[i1];\n",
    "        const REAL sinth = sinth_array[i1];\n",
    "        for(int i2=0;i2<Nxx_plus_2NGHOSTS2-2*NGHOSTS;i2++) {\n",
    "          const REAL ph = ph_array[i2];\n",
    "          // Construct integrand for psi4 spin-weight s=-2,l=2,m=0 spherical harmonic\n",
    "          REAL ReY_sm2_l_m,ImY_sm2_l_m;\n",
    "          SpinWeight_minus2_SphHarmonics(l,m, th,ph,  &ReY_sm2_l_m,&ImY_sm2_l_m);\n",
    "\n",
    "          const int idx2d = i1*(Nxx_plus_2NGHOSTS2-2*NGHOSTS)+i2;\n",
    "          const REAL a = psi4r_at_R_ext[idx2d];\n",
    "          const REAL b = psi4i_at_R_ext[idx2d];\n",
    "          const REAL c = ReY_sm2_l_m;\n",
    "          const REAL d = ImY_sm2_l_m;\n",
    "          psi4r_l_m += (a*c + b*d) * dxx2  * sinth*dxx1;\n",
    "          psi4i_l_m += (b*c - a*d) * dxx2  * sinth*dxx1;\n",
    "        }\n",
    "      }\n",
    "      // Step 4: Output the result of the integration to file.\n",
    "      char filename[100];\n",
    "      sprintf(filename,\"outpsi4_l%d_m%d-r%.2f.txt\",l,m, (double)R_ext);\n",
    "      // If you love \"+\"'s in filenames by all means enable this (ugh):\n",
    "      //if(m>=0) sprintf(filename,\"outpsi4_l%d_m+%d-r%.2f.txt\",l,m, (double)R_ext);\n",
    "      FILE *outpsi4_l_m;\n",
    "      // 0 = n*dt when n=0 is exactly represented in double/long double precision,\n",
    "      //          so no worries about the result being ~1e-16 in double/ld precision\n",
    "      if(curr_time==0) outpsi4_l_m = fopen(filename, \"w\");\n",
    "      else             outpsi4_l_m = fopen(filename, \"a\");\n",
    "      fprintf(outpsi4_l_m,\"%e %.15e %.15e\\n\", (double)(curr_time),\n",
    "              (double)psi4r_l_m,(double)psi4i_l_m);\n",
    "      fclose(outpsi4_l_m);\n",
    "    }\n",
    "  }\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "    desc = \"\"\n",
    "    name = \"driver__spherlikegrids__psi4_spinweightm2_decomposition\"\n",
    "    params = r\"\"\"const paramstruct *restrict params,  REAL *restrict diagnostic_output_gfs,\n",
    "        const int *restrict list_of_R_ext_idxs, const int num_of_R_ext_idxs,\n",
    "        REAL *restrict xx[3],void xx_to_Cart(const paramstruct *restrict params, REAL *restrict xx[3],const int i0,const int i1,const int i2, REAL xCart[3])\"\"\"\n",
    "\n",
    "    body = r\"\"\"  // Step 1: Allocate memory for 2D arrays used to store psi4, theta, sin(theta), and phi.\n",
    "  const int sizeof_2Darray = sizeof(REAL)*(Nxx_plus_2NGHOSTS1-2*NGHOSTS)*(Nxx_plus_2NGHOSTS2-2*NGHOSTS);\n",
    "  REAL *restrict psi4r_at_R_ext = (REAL *restrict)malloc(sizeof_2Darray);\n",
    "  REAL *restrict psi4i_at_R_ext = (REAL *restrict)malloc(sizeof_2Darray);\n",
    "  //         ... also store theta, sin(theta), and phi to corresponding 1D arrays.\n",
    "  REAL *restrict sinth_array = (REAL *restrict)malloc(sizeof(REAL)*(Nxx_plus_2NGHOSTS1-2*NGHOSTS));\n",
    "  REAL *restrict th_array    = (REAL *restrict)malloc(sizeof(REAL)*(Nxx_plus_2NGHOSTS1-2*NGHOSTS));\n",
    "  REAL *restrict ph_array    = (REAL *restrict)malloc(sizeof(REAL)*(Nxx_plus_2NGHOSTS2-2*NGHOSTS));\n",
    "\n",
    "  // Step 2: Loop over all extraction indices:\n",
    "  for(int ii0=0;ii0<num_of_R_ext_idxs;ii0++) {\n",
    "    // Step 2.a: Set the extraction radius R_ext based on the radial index R_ext_idx\n",
    "    REAL R_ext;\n",
    "    {\n",
    "      REAL xCart[3];  xx_to_Cart(params,xx,list_of_R_ext_idxs[ii0],1,1,xCart); // values for itheta and iphi don't matter.\n",
    "      R_ext = sqrt(xCart[0]*xCart[0] + xCart[1]*xCart[1] + xCart[2]*xCart[2]);\n",
    "    }\n",
    "\n",
    "    // Step 2.b: Compute psi_4 at this extraction radius and store to a local 2D array.\n",
    "    const int i0=list_of_R_ext_idxs[ii0];\n",
    "#pragma omp parallel for\n",
    "    for(int i1=NGHOSTS;i1<Nxx_plus_2NGHOSTS1-NGHOSTS;i1++) {\n",
    "      th_array[i1-NGHOSTS]    =     xx[1][i1];\n",
    "      sinth_array[i1-NGHOSTS] = sin(xx[1][i1]);\n",
    "      for(int i2=NGHOSTS;i2<Nxx_plus_2NGHOSTS2-NGHOSTS;i2++) {\n",
    "        ph_array[i2-NGHOSTS] = xx[2][i2];\n",
    "\n",
    "        // Compute real & imaginary parts of psi_4, output to diagnostic_output_gfs\n",
    "        const REAL psi4r = (diagnostic_output_gfs[IDX4S(PSI4_PART0REGF, i0,i1,i2)] +\n",
    "                            diagnostic_output_gfs[IDX4S(PSI4_PART1REGF, i0,i1,i2)] +\n",
    "                            diagnostic_output_gfs[IDX4S(PSI4_PART2REGF, i0,i1,i2)]);\n",
    "        const REAL psi4i = (diagnostic_output_gfs[IDX4S(PSI4_PART0IMGF, i0,i1,i2)] +\n",
    "                            diagnostic_output_gfs[IDX4S(PSI4_PART1IMGF, i0,i1,i2)] +\n",
    "                            diagnostic_output_gfs[IDX4S(PSI4_PART2IMGF, i0,i1,i2)]);\n",
    "\n",
    "        // Store result to \"2D\" array (actually 1D array with 2D storage):\n",
    "        const int idx2d = (i1-NGHOSTS)*(Nxx_plus_2NGHOSTS2-2*NGHOSTS)+(i2-NGHOSTS);\n",
    "        psi4r_at_R_ext[idx2d] = psi4r;\n",
    "        psi4i_at_R_ext[idx2d] = psi4i;\n",
    "      }\n",
    "    }\n",
    "    // Step 3: Perform integrations across all l,m modes from l=2 up to and including L_MAX (global variable):\n",
    "    lowlevel_decompose_psi4_into_swm2_modes(Nxx_plus_2NGHOSTS1,Nxx_plus_2NGHOSTS2,\n",
    "                                            dxx1,dxx2,\n",
    "                                            time, R_ext, th_array, sinth_array, ph_array,\n",
    "                                            psi4r_at_R_ext,psi4i_at_R_ext);\n",
    "  }\n",
    "\n",
    "  // Step 4: Free all allocated memory:\n",
    "  free(psi4r_at_R_ext); free(psi4i_at_R_ext);\n",
    "  free(sinth_array); free(th_array); free(ph_array);\n",
    "\"\"\"\n",
    "\n",
    "    print_msg_with_timing(\"Spin-weight s=-2 Spherical Harmonics\", msg=\"Ccodegen\", startstop=\"stop\", starttime=starttime)\n",
    "\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        prefunc=prefunc,\n",
    "        desc=desc,\n",
    "        name=name, params=params,\n",
    "        body=body,\n",
    "        rel_path_to_Cparams=rel_path_to_Cparams)\n",
    "    return pickle_NRPy_env()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='validation'></a>\n",
    "\n",
    "# Step 7: Confirm above functions are bytecode-identical to those in `BSSN/BSSN_Ccodegen_library.py` \\[Back to [top](#toc)\\]\n",
    "$$\\label{validation}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:28:28.173100Z",
     "iopub.status.busy": "2021-03-07T17:28:28.172204Z",
     "iopub.status.idle": "2021-03-07T17:29:26.259689Z",
     "shell.execute_reply": "2021-03-07T17:29:26.258974Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PASS! ALL FUNCTIONS ARE IDENTICAL\n"
     ]
    }
   ],
   "source": [
    "import BSSN.BSSN_Ccodegen_library as BCL\n",
    "import sys\n",
    "\n",
    "funclist = [(\"print_msg_with_timing\", print_msg_with_timing, BCL.print_msg_with_timing),\n",
    "            (\"get_loopopts\", get_loopopts, BCL.get_loopopts),\n",
    "            (\"register_stress_energy_source_terms_return_T4UU\", register_stress_energy_source_terms_return_T4UU, BCL.register_stress_energy_source_terms_return_T4UU),\n",
    "            (\"BSSN_RHSs__generate_symbolic_expressions\", BSSN_RHSs__generate_symbolic_expressions, BCL.BSSN_RHSs__generate_symbolic_expressions),\n",
    "            (\"add_rhs_eval_to_Cfunction_dict\", add_rhs_eval_to_Cfunction_dict, BCL.add_rhs_eval_to_Cfunction_dict),\n",
    "            (\"Ricci__generate_symbolic_expressions\", Ricci__generate_symbolic_expressions, BCL.Ricci__generate_symbolic_expressions),\n",
    "            (\"add_Ricci_eval_to_Cfunction_dict\", add_Ricci_eval_to_Cfunction_dict, BCL.add_Ricci_eval_to_Cfunction_dict),\n",
    "            (\"BSSN_constraints__generate_symbolic_expressions\", BSSN_constraints__generate_symbolic_expressions, BCL.BSSN_constraints__generate_symbolic_expressions),\n",
    "            (\"add_BSSN_constraints_to_Cfunction_dict\", add_BSSN_constraints_to_Cfunction_dict, BCL.add_BSSN_constraints_to_Cfunction_dict),\n",
    "            (\"add_enforce_detgammahat_constraint_to_Cfunction_dict\", add_enforce_detgammahat_constraint_to_Cfunction_dict, BCL.add_enforce_detgammahat_constraint_to_Cfunction_dict),\n",
    "            (\"add_psi4_part_to_Cfunction_dict\", add_psi4_part_to_Cfunction_dict, BCL.add_psi4_part_to_Cfunction_dict),\n",
    "            (\"add_psi4_tetrad_to_Cfunction_dict\", add_psi4_tetrad_to_Cfunction_dict, BCL.add_psi4_tetrad_to_Cfunction_dict),\n",
    "            (\"add_SpinWeight_minus2_SphHarmonics_to_Cfunction_dict\", add_SpinWeight_minus2_SphHarmonics_to_Cfunction_dict, BCL.add_SpinWeight_minus2_SphHarmonics_to_Cfunction_dict)\n",
    "           ]\n",
    "\n",
    "if sys.version_info.major >= 3:\n",
    "    import inspect, re\n",
    "    for func in funclist:\n",
    "        # https://stackoverflow.com/questions/20059011/check-if-two-python-functions-are-equal\n",
    "        # remove line continuations\n",
    "        s1 = re.sub(r'\\s*\\\\\\n',' ',inspect.getsource(func[1]))\n",
    "        s2 = re.sub(r'\\s*\\\\\\n',' ',inspect.getsource(func[2]))\n",
    "        if s1 != s2:\n",
    "            print(\"inspect.getsource(func[1]):\")\n",
    "            print(s1)\n",
    "            with open('s1.txt','w') as fd:\n",
    "                print(s1,file=fd)\n",
    "            with open('s2.txt','w') as fd:\n",
    "                print(s2,file=fd)\n",
    "            print(\"inspect.getsource(func[2]):\")\n",
    "            print(s2)\n",
    "            print(\"ERROR: function \" + func[0] + \" is not the same as the Ccodegen_library version!\")\n",
    "            sys.exit(1)\n",
    "\n",
    "    print(\"PASS! ALL FUNCTIONS ARE IDENTICAL\")\n",
    "else:\n",
    "    print(\"SORRY CANNOT CHECK FUNCTION IDENTITY WITH PYTHON 2. PLEASE update your Python installation.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 8: 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_time_evolution-C_codegen_library.pdf](Tutorial-BSSN_time_evolution-C_codegen_library.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": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:29:58.971437Z",
     "iopub.status.busy": "2021-03-07T17:29:58.970705Z",
     "iopub.status.idle": "2021-03-07T17:30:04.091527Z",
     "shell.execute_reply": "2021-03-07T17:30:04.092219Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-BSSN_time_evolution-C_codegen_library.tex, and compiled\n",
      "    LaTeX file to PDF file Tutorial-BSSN_time_evolution-\n",
      "    C_codegen_library.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_time_evolution-C_codegen_library\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}