Tutorial-ADM_Initial_Data_Reader__BSSN_Converter.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 in Curvilinear Coordinates Initial Data Reader\n",
    "## Author: Zach Etienne\n",
    "\n",
    "## This notebook introduces a numerical relativity module specializing in initial data reading for curvilinear coordinates and BSSN conversion. The module enables initial data reading, sets quantities for BSSN-based temporal evolution, and converts ADM initial data from spherical or Cartesian to a Cartesian basis, before converting to BSSN variables in the desired reference-metric basis. It also facilitates BSSN vector/tensor transformations from Cartesian to reference-metric basis, including necessary rescaling.\n",
    "\n",
    "**Notebook Status:** <font color='orange'><b> Self-Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). **Additional validation tests may have been performed, but are as yet, undocumented. (TODO)**\n",
    "\n",
    "### NRPy+ Source Code for this module: [BSSN/Tutorial-ADM_Initial_Data_Reader__BSSN_Converter.py](../edit/BSSN/Tutorial-ADM_Initial_Data_Reader__BSSN_Converter.py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction:\n",
    "\n",
    "These C functions take as input the ADM variables\n",
    "\n",
    "$$\\left\\{\\gamma_{ij}, K_{ij}, \\alpha, \\beta^i, B^i\\right\\}$$\n",
    "\n",
    "and optionally stress-energy tensor\n",
    "\n",
    "$$T^{\\mu\\nu}$$\n",
    "\n",
    "in the spherical or Cartesian basis. It is up to the user to either provide the needed function for reading in these data or use one provided within NRPy+. The function takes the form\n",
    "\n",
    "```C\n",
    "void ID_function(const paramstruct *params, const REAL xCart[3],\n",
    "                 const ID_persist_struct *restrict ID_persist,\n",
    "                 initial_data_struct *restrict initial_data)\n",
    "\n",
    "```\n",
    "\n",
    "where\n",
    "\n",
    "* `params` provides NRPy+ input parameters\n",
    "* `xCart[3]` is the input 3-element array storing the Cartesian coordinate `x,y,z` at which `ID_function()` will set the initial data\n",
    "* `ID_persist` is an input `struct` storing any data needed by the initial data solver/interpolator. For example if `ID_function()` performs interpolations on pseudospectral grids, `ID_persist` might store the pseudospectral coefficients.\n",
    "* `initial_data` is the struct that `ID_function()` outputs, storing the ADM variables and optionally the stress-energy tensor in either the spherical or Cartesian basis.\n",
    "\n",
    "`ID_function()` is called from the main driver routine for reading and setting initial data:\n",
    "\n",
    "* in the spherical basis, the main driver routine is `initial_data_reader__convert_to_BSSN_from_ADM_spherical()`, and\n",
    "* in the Cartesian basis, the main driver routine is `initial_data_reader__convert_to_BSSN_from_ADM_Cartesian()`.\n",
    "\n",
    "These driver routines loop over all gridpoints, and at each gridpoint:\n",
    "\n",
    "1. `ID_function()` is first called, setting $\\left\\{\\gamma_{ij}, K_{ij}, \\alpha, \\beta^i, B^i\\right\\}$, and optionally $T^{\\mu\\nu}$ in the spherical or Cartesian basis.\n",
    "1. `ADM_SphorCart_to_Cart()` then converts ADM variables from the spherical or Cartesian basis to the Cartesian basis.\n",
    "1. Next, `ADM_Cart_to_BSSN_Cart()` converts ADM variables in the Cartesian basis to BSSN variables in the Cartesian basis.\n",
    "1. Finally, `initial_data_BSSN_basis_transform_Cartesian_to_rfm()` converts BSSN vectors/tensors from Cartesian to reference-metric basis\n",
    "\n",
    "The above steps will set all the BSSN quantities needed for full evolutions, at all gridpoints, *except* the *derived* BSSN quantity $\\lambda^i$. $\\lambda^i$ is set in terms of finite-difference derivatives of other BSSN quantities.\n",
    "\n",
    "As such, after the other BSSN quantities are set at all gridpoints, the function `initial_data_lambdaU_grid_interior()` is called to compute $\\lambda^i$ at all points in the grid *interior*.\n",
    "\n",
    "**It is up to the user** to then call the inner boundary and extrapolation-outer boundary condition routine from the [Curvilinear boundary conditions driver](Tutorial-Start_to_Finish-Curvilinear_BCs_new_way.ipynb), so that all BSSN quantities are set at all points on the numerical grid, *including* the grid boundaries."
   ]
  },
  {
   "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 core Python/NRPy+ modules\n",
    "1. [Step 2](#id_reader_driver): Overview of `initial_data_reader__convert_ADM_..._to_BSSN()`: Driver routine that computes/reads ADM variables at all points on all grids, and converts them to BSSN quantities in chosen curvilinear reference metric\n",
    "1. [Step 2.a](#example_id_func): Example `ID_function()` generator, for \"exact\" initial data (i.e., initial data in which all quantities are specified with closed-form expressions)\n",
    "1. [Step 2.b](#adm_sphorcart_to_cart): `ADM_SphorCart_to_Cart()`: Convert ADM variables from the spherical or Cartesian basis to the Cartesian basis\n",
    "1. [Step 2.c](#adm_cart_to_bssn_cart) `ADM_Cart_to_BSSN_Cart()`: Convert ADM variables in the Cartesian basis to BSSN variables in the Cartesian basis\n",
    "1. [Step 2.d](#bssn_cart_to_rescaled_bssn_rfm): `BSSN_Cart_to_rescaled_BSSN_rfm()`: Convert Cartesian-basis BSSN vectors/tensors *except* $\\lambda^i$, to the basis specified by `reference_metric::CoordSystem`, then rescale these BSSN quantities.\n",
    "1. [Step 2.e](#lambda): `initial_data_lambdaU_grid_interior()`: Compute $\\lambda^i$ from finite-difference derivatives of rescaled metric quantities\n",
    "1. [Step 2.f](#register_driver_function): Putting it all together: Register `initial_data_reader__convert_ADM_..._to_BSSN()` C function \n",
    "\n",
    "1. [Step 3](#nbd): `register_NRPy_basic_defines()`: Register `ID_data_struct` within `NRPy_basic_defines.h`\n",
    "1. [Step 4](#code_validation):Code Validation against  `BSSN.ADM_Initial_Data_Reader__BSSN_Converter` NRPy+ module\n",
    "1. [Step 5](#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 core Python/NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'sympy'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[1], line 2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# Step 1: Initialize core Python/NRPy+ modules\u001b[39;00m\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01moutputC\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m outputC,lhrh,add_to_Cfunction_dict, Cfunction \u001b[38;5;66;03m# NRPy+: Core C code output module\u001b[39;00m\n\u001b[1;32m      3\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01moutputC\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m outC_NRPy_basic_defines_h_dict\n\u001b[1;32m      4\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mNRPy_param_funcs\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpar\u001b[39;00m    \u001b[38;5;66;03m# NRPy+: Parameter interface\u001b[39;00m\n",
      "File \u001b[0;32m~/Documents/VS_Code/Etienne/nrpytutorial/outputC.py:18\u001b[0m\n\u001b[1;32m     12\u001b[0m __all__ \u001b[38;5;241m=\u001b[39m [\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlhrh\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutCparams\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mnrpyAbs\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124msuperfast_uniq\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mcheck_if_string__error_if_not\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m     13\u001b[0m            \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutputC\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mparse_outCparams_string\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m     14\u001b[0m            \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutC_NRPy_basic_defines_h_dict\u001b[39m\u001b[38;5;124m'\u001b[39m,\n\u001b[1;32m     15\u001b[0m            \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutC_function_prototype_dict\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutC_function_dict\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mCfunction\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124madd_to_Cfunction_dict\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124moutCfunction\u001b[39m\u001b[38;5;124m'\u001b[39m]\n\u001b[1;32m     17\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mloop\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mlp\u001b[39;00m                             \u001b[38;5;66;03m# NRPy+: C code loop interface\u001b[39;00m\n\u001b[0;32m---> 18\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01mNRPy_param_funcs\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01mpar\u001b[39;00m                \u001b[38;5;66;03m# NRPy+: parameter interface\u001b[39;00m\n\u001b[1;32m     19\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mSIMD\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m expr_convert_to_SIMD_intrins \u001b[38;5;66;03m# NRPy+: SymPy expression => SIMD intrinsics interface\u001b[39;00m\n\u001b[1;32m     20\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mcse_helpers\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m cse_preprocess,cse_postprocess  \u001b[38;5;66;03m# NRPy+: CSE preprocessing and postprocessing\u001b[39;00m\n",
      "File \u001b[0;32m~/Documents/VS_Code/Etienne/nrpytutorial/NRPy_param_funcs.py:10\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;66;03m# As documented in the NRPy+ tutorial module\u001b[39;00m\n\u001b[1;32m      2\u001b[0m \u001b[38;5;66;03m#   Tutorial-Coutput__Parameter_Interface.ipynb\u001b[39;00m\n\u001b[1;32m      3\u001b[0m \u001b[38;5;66;03m#   this core NRPy+ module is used for\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[38;5;66;03m# Author: Zachariah B. Etienne\u001b[39;00m\n\u001b[1;32m      8\u001b[0m \u001b[38;5;66;03m#         zachetie **at** gmail **dot* com\u001b[39;00m\n\u001b[0;32m---> 10\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01msympy\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m \u001b[38;5;21;01msp\u001b[39;00m                   \u001b[38;5;66;03m# Import SymPy\u001b[39;00m\n\u001b[1;32m     11\u001b[0m \u001b[38;5;28;01mimport\u001b[39;00m \u001b[38;5;21;01msys\u001b[39;00m                           \u001b[38;5;66;03m# Standard Python: OS-independent system functions\u001b[39;00m\n\u001b[1;32m     12\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mcollections\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m namedtuple   \u001b[38;5;66;03m# Standard Python: Enable namedtuple data type\u001b[39;00m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'sympy'"
     ]
    }
   ],
   "source": [
    "# Step 1: Initialize core Python/NRPy+ modules\n",
    "from outputC import outputC,lhrh,add_to_Cfunction_dict, Cfunction # NRPy+: Core C code output module\n",
    "from outputC import outC_NRPy_basic_defines_h_dict\n",
    "import NRPy_param_funcs as par    # NRPy+: Parameter interface\n",
    "import finite_difference as fin   # NRPy+: Finite difference C code generation module\n",
    "import grid as gri                # NRPy+: Functions having to do with numerical grids\n",
    "import indexedexp as ixp          # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm    # NRPy+: Reference metric support\n",
    "import BSSN.BSSN_quantities as Bq # NRPy+: Computes useful BSSN quantities; e.g., gammabarUU & GammabarUDD needed below\n",
    "import cmdline_helper as cmd      # NRPy+: Multi-platform Python command-line interface\n",
    "from pickling import pickle_NRPy_env # NRPy+: Pickle/unpickle NRPy+ environment, for parallel codegen\n",
    "import os, sys                    # Standard Python modules for multiplatform OS-level functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='id_reader_driver'></a>\n",
    "\n",
    "# Step 2: Overview of `initial_data_reader__convert_ADM_..._to_BSSN()`: Driver routine that computes/reads ADM variables at all points on all grids, and converts them to BSSN quantities in chosen curvilinear reference metric \\[Back to [top](#toc)\\]\n",
    "$$\\label{id_reader_driver}$$\n",
    "\n",
    "Initial data for constructing spacetimes in NRPy+ are provided in either the spherical or Cartesian basis, and the corresponding reader routines are `initial_data_reader__convert_ADM_spherical_to_BSSN()` or `initial_data_reader__convert_ADM_Cartesian_to_BSSN()`, respectively. These routines do the following:\n",
    "\n",
    "1. Call `ID_function()`: Read or compute ADM initial data at Cartesian point `xCart[3]`=$(x,y,z)$, in the spherical basis inside `initial_data_reader__convert_to_BSSN_from_ADM_Spherical()` or  Cartesian basis inside `initial_data_reader__convert_to_BSSN_from_ADM_Cartesian()`.\n",
    "1. Call `ADM_SphorCart_to_Cart()`: Convert ADM variables from the spherical or Cartesian basis to the Cartesian basis.\n",
    "1. Call `ADM_Cart_to_BSSN_Cart()`: Convert ADM variables in the Cartesian basis to BSSN variables in the Cartesian basis.\n",
    "1. Call `BSSN_Cart_to_rescaled_BSSN_rfm()`: Convert all BSSN vectors/tensors *except* $\\lambda^i$ in the Cartesian basis, to the basis specified by `reference_metric::CoordSystem`. Rescale BSSN quantities.\n",
    "1. Call `initial_data_lambdaU_grid_interior()` to compute $\\lambda^i$ in the `reference_metric::CoordSystem` basis, *in the grid interior only*.\n",
    "\n",
    "**Important Note**: After `initial_data_reader__convert_to_BSSN_rfm_from_ADM_sph_or_Cart()` is called, inner/outer boundary conditions must be applied to $\\lambda^i$ to ensure it is specified on the grid boundaries."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='example_id_func'></a>\n",
    "\n",
    "## Step 2.a: Example `ID_function()` generator, for \"exact\" initial data (i.e., initial data in which all quantities are specified with closed-form expressions) \\[Back to [top](#toc)\\]\n",
    "$$\\label{example_id_func}$$\n",
    "\n",
    "This function converts closed-form SymPy expressions for ADM quantities $\\alpha$, $\\beta^i$, $B^i$, $\\gamma_{ij}$ and $K_{ij}$ in the spherical or Cartesian basis at a given point $(x,y,z)$ to optimized C code. It is meant to be passed as the `ID_function()` argument into `initial_data_reader__convert_ADM_..._to_BSSN()`. By not calling it `ID_function()`, we can easily have multiple initial data functions within the same C executable.\n",
    "\n",
    "Regarding the input quantities, a number of initial data sets are provided within NRPy+'s `BSSN` module, including for example:\n",
    "\n",
    "* [`BSSN.UIUCBlackHole`](BSSN/UIUCBlackHole.py) for spinning single black hole initial data, in which the coordinate size of the BH does not shrink to zero in the limit of maximum spin\n",
    "* [`BSSN.ShiftedKerrSchild`](BSSN/ShiftedKerrSchild.py) for spinning single black hole initial data\n",
    "* [`BSSN.StaticTrumpet`](BSSN/StaticTrumpet.py) for static trumpet nonspinning black hole initial data\n",
    "* [`BSSN.BrillLindquist`](BSSN/BrillLindquist.py) for initial data of two nonspinning black holes at rest\n",
    "\n",
    "When any of these initial data set functions is called, it will export the ADM quantities as global variables. For example\n",
    "\n",
    "```python\n",
    "import BSSN.UIUCBlackHole as UIB\n",
    "UIB.UIUCBlackHole()\n",
    "pickled_NRPy_env = \\\n",
    "    add_to_Cfunction_dict_exact_ADM_ID_function(\"UIUCBlackHole\", \"Spherical\", UIB.alphaSph, UIB.betaSphU,\n",
    "                                                UIB.BSphU, UIB.betaSphU, UIB.gammaSphDD, UIB.KSphDD)\n",
    "```\n",
    "\n",
    "where the returned pickled NRPy+ environment `pickled_NRPy_env` can be used to run this function in parallel with other functions (parallelized codegen)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_exact_ADM_ID_function(IDtype, IDCoordSystem, alpha, betaU, BU, gammaDD, KDD):\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    desc = IDtype + \" initial data\"\n",
    "    c_type = \"void\"\n",
    "    name = IDtype\n",
    "    params = \"const paramstruct *params, const REAL xCart[3], const ID_persist_struct *restrict ID_persist, initial_data_struct *restrict initial_data\"\n",
    "    desired_rfm_coord = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\", IDCoordSystem)\n",
    "    rfm.reference_metric()\n",
    "    body = \"\"\n",
    "    if IDCoordSystem == \"Spherical\":\n",
    "        body += r\"\"\"\n",
    "  REAL xx0,xx1,xx2 __attribute__((unused));  // xx2 might be unused in the case of axisymmetric initial data.\n",
    "  {\n",
    "    int unused_Cart_to_i0i1i2[3];\n",
    "    REAL xx[3];\n",
    "    Cart_to_xx_and_nearest_i0i1i2(params, xCart, xx, unused_Cart_to_i0i1i2);\n",
    "    xx0=xx[0];  xx1=xx[1];  xx2=xx[2];\n",
    "  }\n",
    "  const REAL r  = xx0; // Some ID only specify r,th,ph.\n",
    "  const REAL th = xx1;\n",
    "  const REAL ph = xx2;\n",
    "\"\"\"\n",
    "    elif IDCoordSystem == \"Cartesian\":\n",
    "        body += r\"\"\"  const REAL Cartxyz0=xCart[0], Cartxyz1=xCart[1], Cartxyz2=xCart[2];\n",
    "\"\"\"\n",
    "    else:\n",
    "        print(\"add_to_Cfunction_dict_exact_ADM_ID_function() Error: IDCoordSystem == \" + IDCoordSystem + \" unsupported\")\n",
    "        sys.exit(1)\n",
    "    list_of_output_exprs = [alpha]\n",
    "    list_of_output_varnames = [\"initial_data->alpha\"]\n",
    "    for i in range(3):\n",
    "        list_of_output_exprs += [betaU[i]]\n",
    "        list_of_output_varnames += [\"initial_data->betaSphorCartU\" + str(i)]\n",
    "        list_of_output_exprs += [BU[i]]\n",
    "        list_of_output_varnames += [\"initial_data->BSphorCartU\" + str(i)]\n",
    "        for j in range(i, 3):\n",
    "            list_of_output_exprs += [gammaDD[i][j]]\n",
    "            list_of_output_varnames += [\"initial_data->gammaSphorCartDD\" + str(i) + str(j)]\n",
    "            list_of_output_exprs += [KDD[i][j]]\n",
    "            list_of_output_varnames += [\"initial_data->KSphorCartDD\" + str(i) + str(j)]\n",
    "    # Sort the outputs before calling outputC()\n",
    "    # https://stackoverflow.com/questions/9764298/is-it-possible-to-sort-two-listswhich-reference-each-other-in-the-exact-same-w\n",
    "    list_of_output_varnames, list_of_output_exprs = (list(t) for t in zip(*sorted(zip(list_of_output_varnames, list_of_output_exprs))))\n",
    "\n",
    "    body += outputC(list_of_output_exprs, list_of_output_varnames,\n",
    "                    filename=\"returnstring\", params=\"outCverbose=False,includebraces=False,preindent=1\")\n",
    "\n",
    "    # Restore CoordSystem:\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\", desired_rfm_coord)\n",
    "    rfm.reference_metric()\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc, c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        enableCparameters=True)\n",
    "    return pickle_NRPy_env()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='adm_sphorcart_to_cart'></a>\n",
    "\n",
    "## Step 2.b: `ADM_SphorCart_to_Cart()`: Convert ADM variables from the spherical or Cartesian basis to the Cartesian basis \\[Back to [top](#toc)\\]\n",
    "$$\\label{adm_sphorcart_to_cart}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Cfunction_ADM_SphorCart_to_Cart(input_Coord=\"Spherical\", include_T4UU=False):\n",
    "    includes = []\n",
    "\n",
    "    desired_rfm_basis = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
    "\n",
    "    desc = \"Convert ADM variables from the spherical or Cartesian basis to the Cartesian basis\"\n",
    "    c_type = \"static void\"\n",
    "    name = \"ADM_SphorCart_to_Cart\"\n",
    "    params = \"\"\"paramstruct *restrict params, const REAL xCart[3], const initial_data_struct *restrict initial_data,\n",
    "                                  ADM_Cart_basis_struct *restrict ADM_Cart_basis\"\"\"\n",
    "\n",
    "    body = r\"\"\"\n",
    "  // Unpack initial_data for ADM vectors/tensors\n",
    "\"\"\"\n",
    "    for i in [\"betaSphorCartU\", \"BSphorCartU\"]:\n",
    "        for j in range(3):\n",
    "            varname = i + str(j)\n",
    "            body += \"  const REAL \" + varname + \" = initial_data->\" + varname + \";\\n\"\n",
    "        body += \"\\n\"\n",
    "    for i in [\"gammaSphorCartDD\", \"KSphorCartDD\"]:\n",
    "        for j in range(3):\n",
    "            for k in range(j, 3):\n",
    "                varname = i + str(j) + str(k)\n",
    "                body += \"  const REAL \" + varname + \" = initial_data->\" + varname + \";\\n\"\n",
    "        body += \"\\n\"\n",
    "    # Read stress-energy tensor in spherical or Cartesian basis if desired.\n",
    "    if include_T4UU:\n",
    "        for mu in range(4):\n",
    "            for nu in range(mu, 4):\n",
    "                varname = \"T4SphorCartUU\" + str(mu) + str(nu)\n",
    "                body += \"  const REAL \" + varname + \" = initial_data->\" + varname + \";\\n\"\n",
    "        body += \"\\n\"\n",
    "\n",
    "    # Set reference_metric to the input_Coord\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\", input_Coord)\n",
    "    rfm.reference_metric()\n",
    "\n",
    "    body += r\"\"\"\n",
    "  // Perform the basis transform on ADM vectors/tensors from \"\"\" + input_Coord + \"\"\" to Cartesian:\n",
    "\n",
    "  // Set destination xx[3] based on desired xCart[3]\n",
    "  REAL xx0,xx1,xx2;\n",
    "  \"\"\" + outputC(rfm.Cart_to_xx[0:3], [\"xx0\", \"xx1\", \"xx2\"],\n",
    "             filename='returnstring', params=\"includebraces=True\").\\\n",
    "                                      replace(\"Cartx\",\"xCart[0]\").\\\n",
    "                                      replace(\"Carty\",\"xCart[1]\").\\\n",
    "                                      replace(\"Cartz\",\"xCart[2]\")\n",
    "\n",
    "    # Define the input variables:\n",
    "    gammaSphorCartDD = ixp.declarerank2(\"gammaSphorCartDD\", \"sym01\")\n",
    "    KSphorCartDD     = ixp.declarerank2(\"KSphorCartDD\", \"sym01\")\n",
    "    betaSphorCartU = ixp.declarerank1(\"betaSphorCartU\")\n",
    "    BSphorCartU    = ixp.declarerank1(\"BSphorCartU\")\n",
    "    T4SphorCartUU = ixp.declarerank2(\"T4SphorCartUU\", \"sym01\", DIM=4)\n",
    "\n",
    "    # Compute Jacobian to convert to Cartesian coordinates\n",
    "    Jac_dUCart_dDrfmUD, Jac_dUrfm_dDCartUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Cartesian()\n",
    "\n",
    "    gammaCartDD = rfm.basis_transform_tensorDD_from_rfmbasis_to_Cartesian(Jac_dUrfm_dDCartUD, gammaSphorCartDD)\n",
    "    KCartDD = rfm.basis_transform_tensorDD_from_rfmbasis_to_Cartesian(Jac_dUrfm_dDCartUD, KSphorCartDD)\n",
    "    betaCartU = rfm.basis_transform_vectorU_from_rfmbasis_to_Cartesian(Jac_dUCart_dDrfmUD, betaSphorCartU)\n",
    "    BCartU = rfm.basis_transform_vectorU_from_rfmbasis_to_Cartesian(Jac_dUCart_dDrfmUD, BSphorCartU)\n",
    "    T4CartUU = ixp.zerorank2(DIM=4)\n",
    "    if include_T4UU:\n",
    "        T4CartUU = rfm.basis_transform_4tensorUU_from_time_indep_rfmbasis_to_Cartesian(Jac_dUCart_dDrfmUD,\n",
    "                                                                                       T4SphorCartUU)\n",
    "\n",
    "    alpha = sp.symbols(\"initial_data->alpha\", real=True)\n",
    "    list_of_output_exprs    = [alpha]\n",
    "    list_of_output_varnames = [\"ADM_Cart_basis->alpha\"]\n",
    "    for i in range(3):\n",
    "        list_of_output_exprs += [betaCartU[i]]\n",
    "        list_of_output_varnames += [\"ADM_Cart_basis->betaU\" + str(i)]\n",
    "        list_of_output_exprs += [BCartU[i]]\n",
    "        list_of_output_varnames += [\"ADM_Cart_basis->BU\" + str(i)]\n",
    "        for j in range(i, 3):\n",
    "            list_of_output_exprs += [gammaCartDD[i][j]]\n",
    "            list_of_output_varnames += [\"ADM_Cart_basis->gammaDD\" + str(i) + str(j)]\n",
    "            list_of_output_exprs += [KCartDD[i][j]]\n",
    "            list_of_output_varnames += [\"ADM_Cart_basis->KDD\" + str(i) + str(j)]\n",
    "    if include_T4UU:\n",
    "        for mu in range(4):\n",
    "            for nu in range(mu, 4):\n",
    "                list_of_output_exprs += [T4CartUU[mu][nu]]\n",
    "                list_of_output_varnames += [\"ADM_Cart_basis->T4UU\" + str(mu) + str(nu)]\n",
    "\n",
    "    # Sort the outputs before calling outputC()\n",
    "    # https://stackoverflow.com/questions/9764298/is-it-possible-to-sort-two-listswhich-reference-each-other-in-the-exact-same-w\n",
    "    list_of_output_varnames, list_of_output_exprs = (list(t) for t in zip(*sorted(zip(list_of_output_varnames, list_of_output_exprs))))\n",
    "\n",
    "    body += outputC(list_of_output_exprs, list_of_output_varnames,\n",
    "                    filename=\"returnstring\", params=\"outCverbose=False,includebraces=False,preindent=1\")\n",
    "\n",
    "    # Restore reference metric globals to coordsystem on grid.\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\", desired_rfm_basis)\n",
    "    rfm.reference_metric()\n",
    "\n",
    "    _func_prototype, func = Cfunction(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        enableCparameters=False)\n",
    "    return func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='adm_cart_to_bssn_cart'></a>\n",
    "\n",
    "## Step 2.c: `ADM_Cart_to_BSSN_Cart()`: Convert ADM variables in the Cartesian basis to BSSN variables in the Cartesian basis \\[Back to [top](#toc)\\]\n",
    "$$\\label{adm_cart_to_bssn_cart}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Cfunction_ADM_Cart_to_BSSN_Cart(include_T4UU=False):\n",
    "    includes = []\n",
    "\n",
    "    desired_rfm_basis = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
    "\n",
    "    desc = \"Convert ADM variables in the Cartesian basis to BSSN variables in the Cartesian basis\"\n",
    "    c_type = \"static void\"\n",
    "    name = \"ADM_Cart_to_BSSN_Cart\"\n",
    "    params = \"\"\"paramstruct *restrict params, const REAL xCart[3], const ADM_Cart_basis_struct *restrict ADM_Cart_basis,\n",
    "                                  BSSN_Cart_basis_struct *restrict BSSN_Cart_basis\"\"\"\n",
    "\n",
    "    # Extract desired rfm basis from reference_metric::CoordSystem\n",
    "    desired_rfm_basis = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
    "\n",
    "    # Set CoordSystem to Cartesian\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\", \"Cartesian\")\n",
    "    rfm.reference_metric()\n",
    "\n",
    "    gammaCartDD = ixp.declarerank2(\"ADM_Cart_basis->gammaDD\", \"sym01\")\n",
    "    KCartDD     = ixp.declarerank2(\"ADM_Cart_basis->KDD\", \"sym01\")\n",
    "\n",
    "    import BSSN.BSSN_in_terms_of_ADM as BitoA\n",
    "    BitoA.trK_AbarDD_aDD(gammaCartDD, KCartDD)\n",
    "    BitoA.gammabarDD_hDD(gammaCartDD)\n",
    "    BitoA.cf_from_gammaDD(gammaCartDD)\n",
    "\n",
    "    body = r\"\"\"\n",
    "  // *In the Cartesian basis*, convert ADM quantities gammaDD & KDD\n",
    "  //   into BSSN gammabarDD, AbarDD, cf, and trK.\n",
    "  BSSN_Cart_basis->alpha = ADM_Cart_basis->alpha;\n",
    "  BSSN_Cart_basis->betaU0 = ADM_Cart_basis->betaU0;\n",
    "  BSSN_Cart_basis->betaU1 = ADM_Cart_basis->betaU1;\n",
    "  BSSN_Cart_basis->betaU2 = ADM_Cart_basis->betaU2;\n",
    "  BSSN_Cart_basis->BU0 = ADM_Cart_basis->BU0;\n",
    "  BSSN_Cart_basis->BU1 = ADM_Cart_basis->BU1;\n",
    "  BSSN_Cart_basis->BU2 = ADM_Cart_basis->BU2;\n",
    "\"\"\"\n",
    "    list_of_output_exprs    = [BitoA.cf, BitoA.trK]\n",
    "    list_of_output_varnames = [\"BSSN_Cart_basis->cf\", \"BSSN_Cart_basis->trK\"]\n",
    "    for i in range(3):\n",
    "        for j in range(i, 3):\n",
    "            list_of_output_exprs += [BitoA.gammabarDD[i][j]]\n",
    "            list_of_output_varnames += [\"BSSN_Cart_basis->gammabarDD\" + str(i) + str(j)]\n",
    "            list_of_output_exprs += [BitoA.AbarDD[i][j]]\n",
    "            list_of_output_varnames += [\"BSSN_Cart_basis->AbarDD\" + str(i) + str(j)]\n",
    "    if include_T4UU:\n",
    "        T4CartUU = ixp.declarerank2(\"ADM_Cart_basis->T4UU\", \"sym01\", DIM=4)\n",
    "        for mu in range(4):\n",
    "            for nu in range(mu, 4):\n",
    "                list_of_output_exprs += [T4CartUU[mu][nu]]\n",
    "                list_of_output_varnames += [\"BSSN_Cart_basis->T4UU\" + str(mu) + str(nu)]\n",
    "\n",
    "    # Sort the outputs before calling outputC()\n",
    "    # https://stackoverflow.com/questions/9764298/is-it-possible-to-sort-two-listswhich-reference-each-other-in-the-exact-same-w\n",
    "    list_of_output_varnames, list_of_output_exprs = (list(t) for t in zip(*sorted(zip(list_of_output_varnames, list_of_output_exprs))))\n",
    "    body += outputC(list_of_output_exprs, list_of_output_varnames,\n",
    "                    filename=\"returnstring\", params=\"outCverbose=False,includebraces=False,preindent=1\")\n",
    "\n",
    "    # Restore reference metric globals to desired reference metric.\n",
    "    par.set_parval_from_str(\"reference_metric::CoordSystem\", desired_rfm_basis)\n",
    "    rfm.reference_metric()\n",
    "\n",
    "    _func_prototype, func = Cfunction(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        enableCparameters=False)\n",
    "    return func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='bssn_cart_to_rescaled_bssn_rfm'></a>\n",
    "\n",
    "## Step 2.d: `BSSN_Cart_to_rescaled_BSSN_rfm()`: Convert Cartesian-basis BSSN vectors/tensors *except* $\\lambda^i$, to the basis specified by `reference_metric::CoordSystem`, then rescale these BSSN quantities. \\[Back to [top](#toc)\\]\n",
    "$$\\label{bssn_cart_to_rescaled_bssn_rfm}$$\n",
    "\n",
    "By the time this function is called, all BSSN tensors and vectors are in the Cartesian coordinate basis $x^i_{\\rm Cart} = (x,y,z)$, but we need them in the curvilinear coordinate basis $x^i_{\\rm rfm}=$`(xx0,xx1,xx2)` set by the `\"reference_metric::CoordSystem\"` variable. Empirically speaking, it is far easier to write `(x(xx0,xx1,xx2),y(xx0,xx1,xx2),z(xx0,xx1,xx2))` than the inverse, so we will compute the Jacobian matrix\n",
    "\n",
    "$$\n",
    "{\\rm Jac\\_dUCart\\_dDrfmUD[i][j]} = \\frac{\\partial x^i_{\\rm Cart}}{\\partial x^j_{\\rm rfm}},\n",
    "$$\n",
    "\n",
    "via exact differentiation (courtesy SymPy), and the inverse Jacobian\n",
    "$$\n",
    "{\\rm Jac\\_dUrfm\\_dDCartUD[i][j]} = \\frac{\\partial x^i_{\\rm rfm}}{\\partial x^j_{\\rm Cart}},\n",
    "$$\n",
    "\n",
    "using NRPy+'s `generic_matrix_inverter3x3()` function. In terms of these, the transformation of BSSN tensors from Spherical to `\"reference_metric::CoordSystem\"` coordinates may be written:\n",
    "\n",
    "\\begin{align}\n",
    "\\beta^i_{\\rm rfm} &= \\frac{\\partial x^i_{\\rm rfm}}{\\partial x^\\ell_{\\rm Cart}} \\beta^\\ell_{\\rm Cart}\\\\\n",
    "B^i_{\\rm rfm} &= \\frac{\\partial x^i_{\\rm rfm}}{\\partial x^\\ell_{\\rm Cart}} B^\\ell_{\\rm Cart}\\\\\n",
    "\\bar{\\gamma}^{\\rm rfm}_{ij} &= \n",
    "\\frac{\\partial x^\\ell_{\\rm Cart}}{\\partial x^i_{\\rm rfm}}\n",
    "\\frac{\\partial x^m_{\\rm Cart}}{\\partial x^j_{\\rm rfm}} \\bar{\\gamma}^{\\rm Cart}_{\\ell m}\\\\\n",
    "\\bar{A}^{\\rm rfm}_{ij} &= \n",
    "\\frac{\\partial x^\\ell_{\\rm Cart}}{\\partial x^i_{\\rm rfm}}\n",
    "\\frac{\\partial x^m_{\\rm Cart}}{\\partial x^j_{\\rm rfm}} \\bar{A}^{\\rm Cart}_{\\ell m}\n",
    "\\end{align}\n",
    "\n",
    "The above basis transforms are included in functions `basis_transform_tensorDD_from_Cartesian_to_rfmbasis()` and `basis_transform_vectorU_from_Cartesian_to_rfmbasis()` in `reference_metric.py`, and we use them below.\n",
    "\n",
    "After the basis transform has been performed, we perform tensor rescalings to compute the evolved variables $h_{ij}$, $a_{ij}$, $\\text{vet}^i$, and $\\text{bet}^i$:\n",
    "\n",
    "$\\bar{\\gamma}_{ij}$ is rescaled $h_{ij}$ according to the prescription described in the [the covariant BSSN formulation tutorial](Tutorial-BSSN_formulation.ipynb) (also [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):\n",
    "\n",
    "$$\n",
    "h_{ij} = (\\bar{\\gamma}_{ij} - \\hat{\\gamma}_{ij})/\\text{ReDD[i][j]}.\n",
    "$$\n",
    "\n",
    "Further $\\bar{A}_{ij}$, $\\beta^i$, $B^i$ are rescaled to $a_{ij}$, $\\text{vet}^i$, and $\\text{bet}^i$ respectively via the standard formulas (found in [the covariant BSSN formulation tutorial](Tutorial-BSSN_formulation.ipynb); also [Ruchlin *et al.*](https://arxiv.org/pdf/1712.07658.pdf)):\n",
    "\n",
    "\\begin{align}\n",
    "a_{ij} &= \\bar{A}_{ij}/\\text{ReDD[i][j]} \\\\\n",
    "\\text{vet}^i &= \\beta^i/\\text{ReU[i]} \\\\\n",
    "\\text{bet}^i &= B^i/\\text{ReU[i]}.\n",
    "\\end{align}\n",
    "\n",
    "Finally, functions depending on the stress-energy tensor $T^{\\mu\\nu}$ all assume it to be *unrescaled*. Thus we do not rescale $T^{\\mu\\nu}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Cfunction_BSSN_Cart_to_rescaled_BSSN_rfm(rel_path_to_Cparams=os.path.join(\".\"), include_T4UU=False):\n",
    "    includes = []\n",
    "\n",
    "    desc = r\"\"\"Convert Cartesian-basis BSSN vectors/tensors *except* lambda^i,\n",
    "to the basis specified by `reference_metric::CoordSystem`, then rescale these BSSN quantities\"\"\"\n",
    "    c_type = \"static void\"\n",
    "    name = \"BSSN_Cart_to_rescaled_BSSN_rfm\"\n",
    "    params = \"\"\"paramstruct *restrict params, const REAL xCart[3],\n",
    "                                           const BSSN_Cart_basis_struct *restrict BSSN_Cart_basis,\n",
    "                                           rescaled_BSSN_rfm_basis_struct *restrict rescaled_BSSN_rfm_basis\"\"\"\n",
    "\n",
    "    body = r\"\"\"\n",
    "  REAL xx0,xx1,xx2 __attribute__((unused));  // xx2 might be unused in the case of axisymmetric initial data.\n",
    "  {\n",
    "    int unused_Cart_to_i0i1i2[3];\n",
    "    REAL xx[3];\n",
    "    Cart_to_xx_and_nearest_i0i1i2(params, xCart, xx, unused_Cart_to_i0i1i2);\n",
    "    xx0=xx[0];  xx1=xx[1];  xx2=xx[2];\n",
    "  }\n",
    "\"\"\"\n",
    "\n",
    "    # Define the input variables:\n",
    "    gammabarCartDD = ixp.declarerank2(\"BSSN_Cart_basis->gammabarDD\", \"sym01\")\n",
    "    AbarCartDD     = ixp.declarerank2(\"BSSN_Cart_basis->AbarDD\", \"sym01\")\n",
    "    betaCartU = ixp.declarerank1(\"BSSN_Cart_basis->betaU\")\n",
    "    BCartU    = ixp.declarerank1(\"BSSN_Cart_basis->BU\")\n",
    "\n",
    "    # Compute Jacobian to convert to Cartesian coordinates\n",
    "    Jac_dUCart_dDrfmUD, Jac_dUrfm_dDCartUD = rfm.compute_Jacobian_and_inverseJacobian_tofrom_Cartesian()\n",
    "\n",
    "    gammabarDD = rfm.basis_transform_tensorDD_from_Cartesian_to_rfmbasis(Jac_dUCart_dDrfmUD, gammabarCartDD)\n",
    "    AbarDD = rfm.basis_transform_tensorDD_from_Cartesian_to_rfmbasis(Jac_dUCart_dDrfmUD, AbarCartDD)\n",
    "    betaU = rfm.basis_transform_vectorU_from_Cartesian_to_rfmbasis(Jac_dUrfm_dDCartUD, betaCartU)\n",
    "    BU = rfm.basis_transform_vectorU_from_Cartesian_to_rfmbasis(Jac_dUrfm_dDCartUD, BCartU)\n",
    "\n",
    "    # Next rescale:\n",
    "    vetU = ixp.zerorank1()\n",
    "    betU = ixp.zerorank1()\n",
    "    hDD  = ixp.zerorank2()\n",
    "    aDD  = ixp.zerorank2()\n",
    "    for i in range(3):\n",
    "        vetU[i] = betaU[i] / rfm.ReU[i]\n",
    "        betU[i] =    BU[i] / rfm.ReU[i]\n",
    "        for j in range(3):\n",
    "            hDD[i][j] = (gammabarDD[i][j] - rfm.ghatDD[i][j]) / rfm.ReDD[i][j]\n",
    "            aDD[i][j] = AbarDD[i][j] / rfm.ReDD[i][j]\n",
    "\n",
    "    if include_T4UU:\n",
    "        T4CartUU = ixp.declarerank2(\"BSSN_Cart_basis->T4UU\", \"sym01\", DIM=4)\n",
    "        T4UU = rfm.basis_transform_4tensorUU_from_Cartesian_to_time_indep_rfmbasis(Jac_dUrfm_dDCartUD, T4CartUU)\n",
    "\n",
    "    alpha, cf, trK = sp.symbols('BSSN_Cart_basis->alpha BSSN_Cart_basis->cf BSSN_Cart_basis->trK', real=True)\n",
    "\n",
    "    list_of_output_exprs    = [alpha, cf, trK]\n",
    "    list_of_output_varnames = [\"rescaled_BSSN_rfm_basis->alpha\",\n",
    "                               \"rescaled_BSSN_rfm_basis->cf\",\n",
    "                               \"rescaled_BSSN_rfm_basis->trK\"]\n",
    "    for i in range(3):\n",
    "        list_of_output_exprs += [vetU[i]]\n",
    "        list_of_output_varnames += [\"rescaled_BSSN_rfm_basis->vetU\" + str(i)]\n",
    "        list_of_output_exprs += [betU[i]]\n",
    "        list_of_output_varnames += [\"rescaled_BSSN_rfm_basis->betU\" + str(i)]\n",
    "        for j in range(i, 3):\n",
    "            list_of_output_exprs += [hDD[i][j]]\n",
    "            list_of_output_varnames += [\"rescaled_BSSN_rfm_basis->hDD\" + str(i) + str(j)]\n",
    "            list_of_output_exprs += [aDD[i][j]]\n",
    "            list_of_output_varnames += [\"rescaled_BSSN_rfm_basis->aDD\" + str(i) + str(j)]\n",
    "    if include_T4UU:\n",
    "        for mu in range(4):\n",
    "            for nu in range(mu, 4):\n",
    "                # T4UU IS ASSUMED NOT RESCALED; RESCALINGS ARE HANDLED WITHIN BSSN RHSs, etc.\n",
    "                list_of_output_exprs += [T4UU[mu][nu]]\n",
    "                list_of_output_varnames += [\"rescaled_BSSN_rfm_basis->T4UU\" + str(mu) + str(nu)]\n",
    "\n",
    "    # Sort the outputs before calling outputC()\n",
    "    # https://stackoverflow.com/questions/9764298/is-it-possible-to-sort-two-listswhich-reference-each-other-in-the-exact-same-w\n",
    "    list_of_output_varnames, list_of_output_exprs = (list(t) for t in zip(*sorted(zip(list_of_output_varnames, list_of_output_exprs))))\n",
    "\n",
    "    body += outputC(list_of_output_exprs, list_of_output_varnames,\n",
    "                    filename=\"returnstring\", params=\"outCverbose=False,includebraces=False,preindent=1\")\n",
    "\n",
    "    _func_prototype, func = Cfunction(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        enableCparameters=True, rel_path_to_Cparams=rel_path_to_Cparams)\n",
    "    return func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='lambda'></a>\n",
    "\n",
    "## Step 2.e: `initial_data_lambdaU_grid_interior()`: Compute $\\lambda^i$ from finite-difference derivatives of rescaled metric quantities \\[Back to [top](#toc)\\]\n",
    "$$\\label{lambda}$$\n",
    "\n",
    "We compute $\\bar{\\Lambda}^i$ (Eqs. 4 and 5 of [Baumgarte *et al.*](https://arxiv.org/pdf/1211.6632.pdf)), from finite-difference derivatives of rescaled metric quantities $h_{ij}$:\n",
    "\n",
    "$$\n",
    "\\bar{\\Lambda}^i = \\bar{\\gamma}^{jk}\\left(\\bar{\\Gamma}^i_{jk} - \\hat{\\Gamma}^i_{jk}\\right).\n",
    "$$\n",
    "\n",
    "The [reference_metric.py](../edit/reference_metric.py) module provides us with analytic expressions for $\\hat{\\Gamma}^i_{jk}$, so here we need only compute finite-difference expressions for $\\bar{\\Gamma}^i_{jk}$, based on the values for $h_{ij}$ provided in the initial data. Once $\\bar{\\Lambda}^i$ has been computed, we apply the usual rescaling procedure:\n",
    "$$\n",
    "\\lambda^i = \\bar{\\Lambda}^i/\\text{ReU[i]},\n",
    "$$\n",
    "and then output the result to a C file using the NRPy+ finite-difference C output routine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# initial_data_lambdaU_grid_interior() computes lambdaU from\n",
    "#  finite-difference derivatives of rescaled metric quantities\n",
    "def Cfunction_initial_data_lambdaU_grid_interior():\n",
    "    includes = []\n",
    "    c_type = \"static void\"\n",
    "\n",
    "    output_Coord = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
    "    desc = \"Compute lambdaU in \" + output_Coord + \" coordinates\"\n",
    "    name = \"initial_data_lambdaU_grid_interior\"\n",
    "    params = \"\"\"const paramstruct *restrict params, REAL *restrict xx[3], REAL *restrict in_gfs\"\"\"\n",
    "    # Step 7: Compute $\\bar{\\Lambda}^i$ from finite-difference derivatives of rescaled metric quantities\n",
    "\n",
    "    # We will need all BSSN gridfunctions to be defined, as well as\n",
    "    #     expressions for gammabarDD_dD in terms of exact derivatives of\n",
    "    #     the rescaling matrix and finite-difference derivatives of\n",
    "    #     hDD's. This functionality is provided by BSSN.BSSN_unrescaled_and_barred_vars,\n",
    "    #     which we call here to overwrite above definitions of gammabarDD,gammabarUU, etc.\n",
    "    Bq.gammabar__inverse_and_derivs()  # Provides gammabarUU and GammabarUDD\n",
    "    gammabarUU    = Bq.gammabarUU\n",
    "    GammabarUDD   = Bq.GammabarUDD\n",
    "\n",
    "    # Next evaluate \\bar{\\Lambda}^i, based on GammabarUDD above and GammahatUDD\n",
    "    #       (from the reference metric):\n",
    "    LambdabarU = ixp.zerorank1()\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            for k in range(3):\n",
    "                LambdabarU[i] += gammabarUU[j][k] * (GammabarUDD[i][j][k] - rfm.GammahatUDD[i][j][k])\n",
    "\n",
    "    # Finally apply rescaling:\n",
    "    # lambda^i = Lambdabar^i/\\text{ReU[i]}\n",
    "    lambdaU = ixp.zerorank1()\n",
    "    for i in range(3):\n",
    "        lambdaU[i] = LambdabarU[i] / rfm.ReU[i]\n",
    "\n",
    "    lambdaU_expressions = [lhrh(lhs=gri.gfaccess(\"in_gfs\", \"lambdaU0\"), rhs=lambdaU[0]),\n",
    "                           lhrh(lhs=gri.gfaccess(\"in_gfs\", \"lambdaU1\"), rhs=lambdaU[1]),\n",
    "                           lhrh(lhs=gri.gfaccess(\"in_gfs\", \"lambdaU2\"), rhs=lambdaU[2])]\n",
    "    body = fin.FD_outputC(\"returnstring\", lambdaU_expressions,\n",
    "                           params=\"outCverbose=False,includebraces=False,preindent=0\")\n",
    "    _func_prototype, func = Cfunction(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        loopopts=\"InteriorPoints,Read_xxs\",\n",
    "        enableCparameters=True)\n",
    "    return func"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='register_driver_function'></a>\n",
    "\n",
    "## Step 2.f: Putting it all together: Register `initial_data_reader__convert_ADM_..._to_BSSN()` C function \\[Back to [top](#toc)\\]\n",
    "$$\\label{register_driver_function}$$\n",
    "\n",
    "As discussed above, the function `initial_data_reader__convert_ADM_..._to_BSSN()` is our core driver function, generally called directly by `main()` to read/construct initial data in terms of ADM quantities in the spherical/Cartesian basis, and convert these quantities to rescaled BSSN variables. As such, it incorporates all the above functions just described."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:06:02.470552Z",
     "iopub.status.busy": "2021-03-07T17:06:02.469863Z",
     "iopub.status.idle": "2021-03-07T17:06:02.472566Z",
     "shell.execute_reply": "2021-03-07T17:06:02.471938Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_initial_data_reader__convert_ADM_Sph_or_Cart_to_BSSN(addl_includes=None,\n",
    "                                                                               rel_path_to_Cparams=os.path.join(\".\"),\n",
    "                                                                               input_Coord=\"Spherical\",\n",
    "                                                                               include_T4UU=False):\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    if par.parval_from_str(\"finite_difference::enable_FD_functions\"):\n",
    "        includes += [\"finite_difference_functions.h\"]\n",
    "    if addl_includes is not None:\n",
    "        if not isinstance(addl_includes, list):\n",
    "            print(\"Error: addl_includes must be a list.\")\n",
    "            sys.exit(1)\n",
    "        includes += addl_includes\n",
    "\n",
    "    def T4UU_prettyprint():\n",
    "        return r\"\"\"\n",
    "  REAL T4UU00,T4UU01,T4UU02,T4UU03;\n",
    "  REAL        T4UU11,T4UU12,T4UU13;\n",
    "  REAL               T4UU22,T4UU23;\n",
    "  REAL                      T4UU33;\n",
    "\"\"\"\n",
    "    prefunc = \"\"\"\n",
    "// ADM variables in the Cartesian basis:\n",
    "typedef struct __ADM_Cart_basis_struct__ {\n",
    "  REAL alpha, betaU0,betaU1,betaU2, BU0,BU1,BU2;\n",
    "  REAL gammaDD00,gammaDD01,gammaDD02,gammaDD11,gammaDD12,gammaDD22;\n",
    "  REAL KDD00,KDD01,KDD02,KDD11,KDD12,KDD22;\n",
    "\"\"\"\n",
    "    if include_T4UU:\n",
    "        prefunc += T4UU_prettyprint()\n",
    "    prefunc += \"} ADM_Cart_basis_struct;\\n\"\n",
    "    ##############\n",
    "    prefunc += \"\"\"\n",
    "// BSSN variables in the Cartesian basis:\n",
    "typedef struct __BSSN_Cart_basis_struct__ {\n",
    "  REAL alpha, betaU0,betaU1,betaU2, BU0,BU1,BU2;\n",
    "  REAL cf, trK;\n",
    "  REAL gammabarDD00,gammabarDD01,gammabarDD02,gammabarDD11,gammabarDD12,gammabarDD22;\n",
    "  REAL AbarDD00,AbarDD01,AbarDD02,AbarDD11,AbarDD12,AbarDD22;\n",
    "\"\"\"\n",
    "    if include_T4UU:\n",
    "        prefunc += T4UU_prettyprint()\n",
    "    prefunc += \"} BSSN_Cart_basis_struct;\\n\"\n",
    "    ##############\n",
    "    prefunc += \"\"\"\n",
    "// Rescaled BSSN variables in the rfm basis:\n",
    "typedef struct __rescaled_BSSN_rfm_basis_struct__ {\n",
    "  REAL alpha, vetU0,vetU1,vetU2, betU0,betU1,betU2;\n",
    "  REAL cf, trK;\n",
    "  REAL hDD00,hDD01,hDD02,hDD11,hDD12,hDD22;\n",
    "  REAL aDD00,aDD01,aDD02,aDD11,aDD12,aDD22;\n",
    "\"\"\"\n",
    "    if include_T4UU:\n",
    "        prefunc += T4UU_prettyprint()\n",
    "    prefunc += \"} rescaled_BSSN_rfm_basis_struct;\\n\"\n",
    "    ##############\n",
    "    ##############\n",
    "    prefunc += Cfunction_ADM_SphorCart_to_Cart(input_Coord=input_Coord, include_T4UU=include_T4UU)\n",
    "    prefunc += Cfunction_ADM_Cart_to_BSSN_Cart(                         include_T4UU=include_T4UU)\n",
    "    prefunc += Cfunction_BSSN_Cart_to_rescaled_BSSN_rfm(rel_path_to_Cparams=rel_path_to_Cparams,\n",
    "                                                        include_T4UU=include_T4UU)\n",
    "    prefunc += Cfunction_initial_data_lambdaU_grid_interior()\n",
    "\n",
    "    output_Coord = par.parval_from_str(\"reference_metric::CoordSystem\")\n",
    "    desc = \"Read in ADM initial data in the \" + input_Coord + \" basis, and convert to BSSN data in \" + output_Coord + \" coordinates\"\n",
    "    c_type = \"void\"\n",
    "    name = \"initial_data_reader__convert_ADM_\" + input_Coord + \"_to_BSSN\"\n",
    "    params = \"\"\"griddata_struct *restrict griddata, ID_persist_struct *restrict ID_persist,\n",
    "                                                             void ID_function(const paramstruct *params, const REAL xCart[3],\n",
    "                                                                              const ID_persist_struct *restrict ID_persist,\n",
    "                                                                              initial_data_struct *restrict initial_data)\"\"\"\n",
    "\n",
    "    body = r\"\"\"\n",
    "  const int Nxx_plus_2NGHOSTS0 = griddata->params.Nxx_plus_2NGHOSTS0;\n",
    "  const int Nxx_plus_2NGHOSTS1 = griddata->params.Nxx_plus_2NGHOSTS1;\n",
    "  const int Nxx_plus_2NGHOSTS2 = griddata->params.Nxx_plus_2NGHOSTS2;\n",
    "\n",
    "  LOOP_OMP(\"omp parallel for\", i0,0,Nxx_plus_2NGHOSTS0, i1,0,Nxx_plus_2NGHOSTS1, i2,0,Nxx_plus_2NGHOSTS2) {\n",
    "    // xCart is the global Cartesian coordinate, which accounts for any grid offsets from the origin.\n",
    "    REAL xCart[3];  xx_to_Cart(&griddata->params, griddata->xx, i0,i1,i2, xCart);\n",
    "\n",
    "    // Read or compute initial data at destination point xCart\n",
    "    initial_data_struct initial_data;\n",
    "    ID_function(&griddata->params, xCart, ID_persist, &initial_data);\n",
    "\n",
    "    ADM_Cart_basis_struct ADM_Cart_basis;\n",
    "    ADM_SphorCart_to_Cart(&griddata->params, xCart, &initial_data, &ADM_Cart_basis);\n",
    "\n",
    "    BSSN_Cart_basis_struct BSSN_Cart_basis;\n",
    "    ADM_Cart_to_BSSN_Cart(&griddata->params, xCart, &ADM_Cart_basis, &BSSN_Cart_basis);\n",
    "\n",
    "    rescaled_BSSN_rfm_basis_struct rescaled_BSSN_rfm_basis;\n",
    "    BSSN_Cart_to_rescaled_BSSN_rfm(&griddata->params, xCart, &BSSN_Cart_basis, &rescaled_BSSN_rfm_basis);\n",
    "\n",
    "    const int idx3 = IDX3S(i0,i1,i2);\n",
    "\n",
    "    // Output data to gridfunctions\n",
    "\"\"\"\n",
    "    gf_list = [\"alpha\", \"trK\", \"cf\"]\n",
    "    for i in range(3):\n",
    "        gf_list += [\"vetU\"+str(i), \"betU\"+str(i)]\n",
    "        for j in range(i, 3):\n",
    "            gf_list += [\"hDD\"+str(i)+str(j), \"aDD\"+str(i)+str(j)]\n",
    "    for gf in sorted(gf_list):\n",
    "        body += \"    griddata->gridfuncs.y_n_gfs[IDX4ptS(\"+gf.upper()+\"GF, idx3)] = rescaled_BSSN_rfm_basis.\"+gf+\";\\n\"\n",
    "    if include_T4UU:\n",
    "        for mu in range(4):\n",
    "            for nu in range(mu, 4):\n",
    "                gf = \"T4UU\" + str(mu) + str(nu)\n",
    "                body += \"    griddata->gridfuncs.auxevol_gfs[IDX4ptS(\"+gf.upper()+\"GF, idx3)] = rescaled_BSSN_rfm_basis.\"+gf+\";\\n\"\n",
    "    body += \"\"\"\n",
    "  } // END LOOP over all gridpoints on given grid\n",
    "\n",
    "  initial_data_lambdaU_grid_interior(&griddata->params, griddata->xx, griddata->gridfuncs.y_n_gfs);\n",
    "\"\"\"\n",
    "\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        prefunc=prefunc,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=body,\n",
    "        enableCparameters=False)\n",
    "    return pickle_NRPy_env()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nbd'></a>\n",
    "\n",
    "# Step 3: `register_NRPy_basic_defines()`: Register `ID_data_struct` within `NRPy_basic_defines.h` \\[Back to [top](#toc)\\]\n",
    "$$\\label{nbd}$$\n",
    "\n",
    "Other than its core use as a means to store ADM input quantities, `ID_data_struct` is designed to be extensible. For example, it may be used to store e.g., pseudospectral coefficients for TwoPunctures, initial data gridfunctions from NRPyElliptic, pointers to TOV 1D data from the TOV solver, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Other than its core use as a means to store ADM input quantities,\n",
    "# `initial_data_struct` is designed to be extensible. For example, it may be\n",
    "# used to store e.g., pseudospectral coefficients for TwoPunctures,\n",
    "# initial data gridfunctions from NRPyElliptic, pointers to TOV 1D data\n",
    "# from the TOV solver, etc.\n",
    "def register_NRPy_basic_defines(ID_persist_struct_contents_str=\"\", include_T4UU=False):\n",
    "    Nbd = r\"\"\"typedef struct __initial_data_struct__ {\n",
    "  REAL alpha;\n",
    "\n",
    "  REAL betaSphorCartU0, betaSphorCartU1, betaSphorCartU2;\n",
    "  REAL BSphorCartU0, BSphorCartU1, BSphorCartU2;\n",
    "\n",
    "  REAL gammaSphorCartDD00, gammaSphorCartDD01, gammaSphorCartDD02;\n",
    "  REAL gammaSphorCartDD11, gammaSphorCartDD12, gammaSphorCartDD22;\n",
    "\n",
    "  REAL KSphorCartDD00, KSphorCartDD01, KSphorCartDD02;\n",
    "  REAL KSphorCartDD11, KSphorCartDD12, KSphorCartDD22;\n",
    "\"\"\"\n",
    "    if include_T4UU:\n",
    "        Nbd += \"\"\"\n",
    "  REAL T4SphorCartUU00,T4SphorCartUU01,T4SphorCartUU02,T4SphorCartUU03;\n",
    "  REAL                 T4SphorCartUU11,T4SphorCartUU12,T4SphorCartUU13;\n",
    "  REAL                                 T4SphorCartUU22,T4SphorCartUU23;\n",
    "  REAL                                                 T4SphorCartUU33;\n",
    "\"\"\"\n",
    "    Nbd += \"\"\"\n",
    "} initial_data_struct;\n",
    "\"\"\"\n",
    "\n",
    "    Nbd += \"typedef struct __ID_persist_struct__ {\\n\"\n",
    "    Nbd += ID_persist_struct_contents_str + \"\\n\"\n",
    "    Nbd += \"} ID_persist_struct;\\n\"\n",
    "    outC_NRPy_basic_defines_h_dict[\"BSSN_initial_data\"] = Nbd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 4: Code Validation against  `BSSN.ADM_Initial_Data_Reader__BSSN_Converter` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here, as a code validation check, we verify that functions within\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [BSSN.ADM_Initial_Data_Reader__BSSN_Converter](../edit/BSSN/ADM_Initial_Data_Reader__BSSN_Converter.py) module\n",
    "\n",
    "are *identical*. This notebook serves as documentation for the Python module, and the Python module is meant to be called from external routines. Thus this validation test ensures the documentation remains consistent with the Python module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:06:04.854546Z",
     "iopub.status.busy": "2021-03-07T17:06:04.828944Z",
     "iopub.status.idle": "2021-03-07T17:06:05.909998Z",
     "shell.execute_reply": "2021-03-07T17:06:05.909451Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PASS! ALL FUNCTIONS ARE IDENTICAL\n"
     ]
    }
   ],
   "source": [
    "import BSSN.ADM_Initial_Data_Reader__BSSN_Converter as AID\n",
    "\n",
    "funclist = [(add_to_Cfunction_dict_exact_ADM_ID_function, AID.add_to_Cfunction_dict_exact_ADM_ID_function),\n",
    "            (Cfunction_ADM_SphorCart_to_Cart, AID.Cfunction_ADM_SphorCart_to_Cart),\n",
    "            (Cfunction_ADM_Cart_to_BSSN_Cart, AID.Cfunction_ADM_Cart_to_BSSN_Cart),\n",
    "            (Cfunction_BSSN_Cart_to_rescaled_BSSN_rfm, AID.Cfunction_BSSN_Cart_to_rescaled_BSSN_rfm),\n",
    "            (Cfunction_initial_data_lambdaU_grid_interior, AID.Cfunction_initial_data_lambdaU_grid_interior),\n",
    "            (add_to_Cfunction_dict_initial_data_reader__convert_ADM_Sph_or_Cart_to_BSSN, AID.add_to_Cfunction_dict_initial_data_reader__convert_ADM_Sph_or_Cart_to_BSSN),\n",
    "            (register_NRPy_basic_defines, AID.register_NRPy_basic_defines)\n",
    "           ]\n",
    "\n",
    "import inspect\n",
    "for func in funclist:\n",
    "    # https://stackoverflow.com/questions/20059011/check-if-two-python-functions-are-equal\n",
    "    if inspect.getsource(func[0]) != inspect.getsource(func[1]):\n",
    "        with open(func[0].__name__ + \"_Jupyter_notebook_version.c\", \"w\") as file:\n",
    "            file.write(inspect.getsource(func[0]))\n",
    "        with open(func[1].__name__ + \"_Python_module_version.c\", \"w\") as file:\n",
    "            file.write(inspect.getsource(func[1]))\n",
    "        print(\"ERROR: function \" + func[0].__name__ + \" is not the same as the Ccodegen_library version!\")\n",
    "        print(\" For more info, try this:\")\n",
    "        print(\"diff \" + func[0].__name__ + \"_Jupyter_notebook_version.c\" + \" \" + func[1].__name__ + \"_Python_module_version.c\")\n",
    "        sys.exit(1)\n",
    "\n",
    "print(\"PASS! ALL FUNCTIONS ARE IDENTICAL\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "Once the following code finishes running, the generated PDF may be found at the following location within the directory you have the NRPy+ tutorial saved:\n",
    "[Tutorial-ADM_Initial_Data_Reader__BSSN_Converter.pdf](Tutorial-ADM_Initial_Data_Reader__BSSN_Converter.pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:06:05.916835Z",
     "iopub.status.busy": "2021-03-07T17:06:05.914038Z",
     "iopub.status.idle": "2021-03-07T17:06:09.387723Z",
     "shell.execute_reply": "2021-03-07T17:06:09.388501Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ADM_Initial_Data_Reader__BSSN_Converter.tex, and compiled\n",
      "    LaTeX file to PDF file Tutorial-\n",
      "    ADM_Initial_Data_Reader__BSSN_Converter.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ADM_Initial_Data_Reader__BSSN_Converter\")"
   ]
  }
 ],
 "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
}