{
 "cells": [
  {
   "attachments": {},
   "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",
    "# Weyl Scalars and Invariants: An Introduction to Einstein Toolkit Diagnostic Thorns\n",
    "\n",
    "## Author: Patrick Nelson & Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This notebook outlines the construction of an Einstein Toolkit diagnostic module to compute the Weyl scalars and invariants using ADMBase gridfunctions as input, with expressions generated through SymPy and NRPy+. The process includes converting SymPy expressions into C-code kernels and integrating these within the larger Einstein Toolkit infrastructure, providing a methodological approach to develop and integrate custom diagnostic tools.\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** Numerical results from this module have been confirmed to agree with the trusted WeylScal4 Einstein Toolkit thorn to roundoff error.\n",
    "\n",
    "### NRPy+ Source Code for this module:  \n",
    "* [WeylScal4NRPD/WeylScalars_Cartesian.py](../edit/WeylScal4NRPD/WeylScalars_Cartesian.py)\n",
    "* [WeylScal4NRPD/WeylScalarInvariants_Cartesian.py](../edit/WeylScal4NRPD/WeylScalarInvariants_Cartesian.py)\n",
    "\n",
    "which are fully documented in the NRPy+ [Tutorial-WeylScalars-Cartesian](Tutorial-WeylScalars-Cartesian.ipynb) module on using NRPy+ to construct the Weyl scalars and invariants as SymPy expressions.\n",
    "\n",
    "## Introduction:\n",
    "In the [previous tutorial notebook](Tutorial-WeylScalars-Cartesian.ipynb), we constructed within SymPy full expressions for the real and imaginary components of all five Weyl scalars $\\psi_0$, $\\psi_1$, $\\psi_2$, $\\psi_3$, and $\\psi_4$ as well as the Weyl invariants. So that we can easily access these expressions, we have ported the Python code needed to generate the Weyl scalar SymPy expressions to [WeylScal4NRPD/WeylScalars_Cartesian.py](../edit/WeylScal4NRPD/WeylScalars_Cartesian.py), and the Weyl invariant SymPy expressions to [WeylScal4NRPD/WeylScalarInvariants_Cartesian.py](../edit/WeylScal4NRPD/WeylScalarInvariants_Cartesian.py).\n",
    "\n",
    "Here we will work through the steps necessary to construct an Einstein Toolkit diagnostic thorn (module), starting from these SymPy expressions, which computes these expressions using ADMBase gridfunctions as input. This tutorial is in two steps:\n",
    "\n",
    "1. Call on NRPy+ to convert the SymPy expressions for the Weyl Scalars and associated Invariants into one C-code kernel for each.\n",
    "1. Write the C code and build up the needed Einstein Toolkit infrastructure (i.e., the .ccl files)."
   ]
  },
  {
   "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](#nrpy): Call on NRPy+ to convert the SymPy expressions for the Weyl scalars and associated invariants into one C-code kernel for each\n",
    "1. [Step 2](#etk): Interfacing with the Einstein Toolkit\n",
    "    1. [Step 2.a](#etkc): Constructing the Einstein Toolkit C-code calling functions that include the C code kernels\n",
    "    1. [Step 2.b](#cclfiles): CCL files - Define how this module interacts and interfaces with the larger Einstein Toolkit infrastructure\n",
    "    1. [Step 2.c](#etk_list): Add the C file to Einstein Toolkit compilation list\n",
    "1. [Step 3](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpy'></a>\n",
    "\n",
    "# Step 1: Call on NRPy+ to convert the SymPy expressions for the Weyl scalars and associated invariants into one C-code kernel for each \\[Back to [top](#toc)\\]\n",
    "$$\\label{nrpy}$$\n",
    "\n",
    "<font color='red'><b>WARNING</b></font>: It takes some time to generate the CSE-optimized C code kernels for these quantities, especially the Weyl scalars... expect 5 minutes on a modern computer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating symbolic expressions for psi4...\n",
      "(BENCH) Finished psi4 symbolic expressions in 1.7516753673553467 seconds.\n",
      "Generating C code kernel for psi4r...\n",
      "(BENCH) Finished psi4r C code kernel generation in 40.00490069389343 seconds.\n"
     ]
    }
   ],
   "source": [
    "from outputC import *            # NRPy+: Core C code output module\n",
    "import finite_difference as fin  # NRPy+: Finite difference C code generation module\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import grid as gri               # NRPy+: Functions having to do with numerical grids\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface\n",
    "import loop as lp                # NRPy+: loop infrasructure\n",
    "import shutil, os, sys, time     # Standard Python modules for multiplatform OS-level functions, benchmarking\n",
    "\n",
    "# Step 1: Set the coordinate system for the numerical grid to Cartesian.\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n",
    "rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.\n",
    "\n",
    "# Step 2: Set the finite differencing order to FD_order to 4\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", 4)\n",
    "\n",
    "# Step 3: Create output directories\n",
    "!mkdir WeylScal4NRPD     2>/dev/null # 2>/dev/null: Don't throw an error or warning if the directory already exists.\n",
    "!mkdir WeylScal4NRPD/src 2>/dev/null # 2>/dev/null: Don't throw an error or warning if the directory already exists.\n",
    "\n",
    "# Step 4: Generate symbolic expressions\n",
    "# Since we are writing an Einstein Toolkit thorn, we must set our memory access style to \"ETK\".\n",
    "par.set_parval_from_str(\"grid::GridFuncMemAccess\",\"ETK\")\n",
    "import BSSN.Psi4_tetrads as BP4t\n",
    "par.set_parval_from_str(\"BSSN.Psi4_tetrads::TetradChoice\",\"QuasiKinnersley\")\n",
    "#par.set_parval_from_str(\"BSSN.Psi4_tetrads::UseCorrectUnitNormal\",\"True\")\n",
    "import BSSN.Psi4 as BP4\n",
    "print(\"Generating symbolic expressions for psi4...\")\n",
    "start = time.time()\n",
    "BP4.Psi4()\n",
    "end = time.time()\n",
    "print(\"(BENCH) Finished psi4 symbolic expressions in \"+str(end-start)+\" seconds.\")\n",
    "\n",
    "psi4r = gri.register_gridfunctions(\"AUX\",\"psi4r\")\n",
    "psi4r0pt = gri.register_gridfunctions(\"AUX\",\"psi4r0pt\")\n",
    "psi4r1pt = gri.register_gridfunctions(\"AUX\",\"psi4r1pt\")\n",
    "psi4r2pt = gri.register_gridfunctions(\"AUX\",\"psi4r2pt\")\n",
    "\n",
    "# Construct RHSs:\n",
    "psi4r_lhrh = [lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi4r\"),rhs=BP4.psi4_re_pt[0]+BP4.psi4_re_pt[1]+BP4.psi4_re_pt[2]),\n",
    "              lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi4r0pt\"),rhs=BP4.psi4_re_pt[0]),\n",
    "              lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi4r1pt\"),rhs=BP4.psi4_re_pt[1]),\n",
    "              lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi4r2pt\"),rhs=BP4.psi4_re_pt[2])]\n",
    "\n",
    "# Generating the CSE is the slowest\n",
    "# operation in this notebook, and much of the CSE\n",
    "# time is spent sorting CSE expressions. Disabling\n",
    "# this sorting makes the C codegen 3-4x faster,\n",
    "# but the tradeoff is that every time this is\n",
    "# run, the CSE patterns will be different \n",
    "# (though they should result in mathematically \n",
    "# *identical* expressions). You can expect \n",
    "# roundoff-level differences as a result.\n",
    "start = time.time()\n",
    "print(\"Generating C code kernel for psi4r...\")\n",
    "psi4r_CcodeKernel = fin.FD_outputC(\"returnstring\",psi4r_lhrh,params=\"outCverbose=False,CSE_sorting=none\")\n",
    "end = time.time()\n",
    "print(\"(BENCH) Finished psi4r C code kernel generation in \"+str(end-start)+\" seconds.\")\n",
    "psi4r_looped = lp.loop([\"i2\",\"i1\",\"i0\"],[\"2\",\"2\",\"2\"],[\"cctk_lsh[2]-2\",\"cctk_lsh[1]-2\",\"cctk_lsh[0]-2\"],\\\n",
    "                       [\"1\",\"1\",\"1\"],[\"#pragma omp parallel for\",\"\",\"\"],\"\",\"\"\"\n",
    "               const CCTK_REAL xx0 = xGF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "               const CCTK_REAL xx1 = yGF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "               const CCTK_REAL xx2 = zGF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "\"\"\"+psi4r_CcodeKernel)\n",
    "with open(\"WeylScal4NRPD/src/WeylScal4NRPD_psi4r.h\", \"w\") as file:\n",
    "    file.write(str(psi4r_looped))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='etk'></a>\n",
    "\n",
    "# Step 2: Interfacing with the Einstein Toolkit \\[Back to [top](#toc)\\]\n",
    "$$\\label{etk}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='etkc'></a>\n",
    "\n",
    "## Step 2.a: Constructing the Einstein Toolkit calling functions that include the C code kernels \\[Back to [top](#toc)\\]\n",
    "$$\\label{etkc}$$\n",
    "\n",
    "Now that we have generated the C code kernels (`WeylScal4NRPD_psis.h` and `WeylScal4NRPD_invars.h`) express the Weyl scalars and invariants as CSE-optimized finite-difference expressions, we next need to write the C code functions that incorporate these kernels and are called by the Einstein Toolkit scheduler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting WeylScal4NRPD/src/WeylScal4NRPD.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile WeylScal4NRPD/src/WeylScal4NRPD.c\n",
    "\n",
    "#include <math.h>\n",
    "#include <stdio.h>\n",
    "#include <stdlib.h>\n",
    "#include <string.h>\n",
    "#include \"cctk.h\"\n",
    "#include \"cctk_Arguments.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "\n",
    "void WeylScal4NRPD_calc_psi4r(const cGH* restrict const cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,\n",
    "const CCTK_REAL invdx0,const CCTK_REAL invdx1,const CCTK_REAL invdx2,\n",
    "const CCTK_REAL *xGF,const CCTK_REAL *yGF,const CCTK_REAL *zGF,\n",
    "const CCTK_REAL *hDD00GF,const CCTK_REAL *hDD01GF,const CCTK_REAL *hDD02GF,const CCTK_REAL *hDD11GF,const CCTK_REAL *hDD12GF,const CCTK_REAL *hDD22GF,\n",
    "const CCTK_REAL *aDD00GF,const CCTK_REAL *aDD01GF,const CCTK_REAL *aDD02GF,const CCTK_REAL *aDD11GF,const CCTK_REAL *aDD12GF,const CCTK_REAL *aDD22GF,\n",
    "const CCTK_REAL *trKGF,const CCTK_REAL *cfGF,\n",
    "CCTK_REAL *psi4rGF,\n",
    "CCTK_REAL *psi4r0ptGF,\n",
    "CCTK_REAL *psi4r1ptGF,\n",
    "CCTK_REAL *psi4r2ptGF) {\n",
    "\n",
    "    DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "#include \"WeylScal4NRPD_psi4r.h\"\n",
    "\n",
    "}\n",
    "\n",
    "extern void WeylScal4NRPD_mainfunction(CCTK_ARGUMENTS) {\n",
    "\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "\n",
    "  if(cctk_iteration % WeylScal4NRPD_calc_every != 0) { return; }\n",
    "  \n",
    "  const CCTK_REAL invdx0 = 1.0 / (CCTK_DELTA_SPACE(0));\n",
    "  const CCTK_REAL invdx1 = 1.0 / (CCTK_DELTA_SPACE(1));\n",
    "  const CCTK_REAL invdx2 = 1.0 / (CCTK_DELTA_SPACE(2));\n",
    "\n",
    "  /* Now, to calculate psi4: */\n",
    "  WeylScal4NRPD_calc_psi4r(cctkGH,cctk_lsh,cctk_nghostzones,\n",
    "                             invdx0,invdx1,invdx2,\n",
    "                             x,y,z,\n",
    "                             hDD00GF,hDD01GF,hDD02GF,hDD11GF,hDD12GF,hDD22GF,\n",
    "                             aDD00GF,aDD01GF,aDD02GF,aDD11GF,aDD12GF,aDD22GF,\n",
    "                             trKGF,cfGF,\n",
    "                             psi4rGF,\n",
    "                             psi4r0ptGF,psi4r1ptGF,psi4r2ptGF);\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First we convert from ADM to BSSN, as is required to convert initial data \n",
    "#    (given using) ADM quantities, to the BSSN evolved variables\n",
    "import BSSN.ADM_Numerical_Spherical_or_Cartesian_to_BSSNCurvilinear as atob\n",
    "IDhDD,IDaDD,IDtrK,IDvetU,IDbetU,IDalpha,IDcf,IDlambdaU = \\\n",
    "    atob.Convert_Spherical_or_Cartesian_ADM_to_BSSN_curvilinear(\"Cartesian\",\"DoNotOutputADMInputFunction\",os.path.join(\"WeylScal4NRPD\",\"src\"))\n",
    "\n",
    "# Store the original list of registered gridfunctions; we'll want to unregister\n",
    "#   all the *SphorCart* gridfunctions after we're finished with them below.\n",
    "orig_glb_gridfcs_list = []\n",
    "for gf in gri.glb_gridfcs_list:\n",
    "    orig_glb_gridfcs_list.append(gf)\n",
    "\n",
    "alphaSphorCart   = gri.register_gridfunctions(                 \"AUXEVOL\", \"alphaSphorCart\")\n",
    "betaSphorCartU   = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\", \"betaSphorCartU\")\n",
    "BSphorCartU      = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\", \"BSphorCartU\")\n",
    "gammaSphorCartDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\", \"gammaSphorCartDD\", \"sym01\")\n",
    "KSphorCartDD     = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\", \"KSphorCartDD\", \"sym01\")\n",
    "\n",
    "# ADM to BSSN conversion, used for converting ADM initial data into a form readable by this thorn.\n",
    "# ADM to BSSN, Part 1: Set up function call and pointers to ADM gridfunctions\n",
    "outstr = \"\"\"\n",
    "#include <math.h>\n",
    "\n",
    "#include \"cctk.h\"\n",
    "#include \"cctk_Arguments.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "\n",
    "void WeylScal4NRPD_ADM_to_BSSN(CCTK_ARGUMENTS) {\n",
    "    DECLARE_CCTK_ARGUMENTS;\n",
    "    DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "    CCTK_REAL *alphaSphorCartGF = alp;\n",
    "\"\"\"\n",
    "# It's ugly if we output code in the following ordering, so we'll first\n",
    "#   output to a string and then sort the string to beautify the code a bit.\n",
    "outstrtmp = []\n",
    "for i in range(3):\n",
    "    outstrtmp.append(\"    CCTK_REAL *betaSphorCartU\"+str(i)+\"GF = beta\"+chr(ord('x')+i)+\";\\n\")\n",
    "#     outstrtmp.append(\"    CCTK_REAL *BSphorCartU\"+str(i)+\"GF = dtbeta\"+chr(ord('x')+i)+\";\\n\")\n",
    "    for j in range(i,3):\n",
    "        outstrtmp.append(\"    CCTK_REAL *gammaSphorCartDD\"+str(i)+str(j)+\"GF = g\"+chr(ord('x')+i)+chr(ord('x')+j)+\";\\n\")\n",
    "        outstrtmp.append(\"    CCTK_REAL *KSphorCartDD\"+str(i)+str(j)+\"GF = k\"+chr(ord('x')+i)+chr(ord('x')+j)+\";\\n\")\n",
    "outstrtmp.sort()\n",
    "for line in outstrtmp:\n",
    "    outstr += line\n",
    "\n",
    "# ADM to BSSN, Part 2: Set up ADM to BSSN conversions for BSSN gridfunctions that do not require\n",
    "#                      finite-difference derivatives (i.e., all gridfunctions except lambda^i (=Gamma^i \n",
    "#                      in non-covariant BSSN)):\n",
    "#                      h_{ij}, a_{ij}, trK, vet^i=beta^i,bet^i=B^i, cf (conformal factor), and alpha\n",
    "all_but_lambdaU_expressions = [\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD00\"),rhs=IDhDD[0][0]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD01\"),rhs=IDhDD[0][1]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD02\"),rhs=IDhDD[0][2]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD11\"),rhs=IDhDD[1][1]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD12\"),rhs=IDhDD[1][2]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"hDD22\"),rhs=IDhDD[2][2]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD00\"),rhs=IDaDD[0][0]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD01\"),rhs=IDaDD[0][1]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD02\"),rhs=IDaDD[0][2]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD11\"),rhs=IDaDD[1][1]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD12\"),rhs=IDaDD[1][2]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"aDD22\"),rhs=IDaDD[2][2]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"trK\"),rhs=IDtrK),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"vetU0\"),rhs=IDvetU[0]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"vetU1\"),rhs=IDvetU[1]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"vetU2\"),rhs=IDvetU[2]),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"alpha\"),rhs=IDalpha),\n",
    "    lhrh(lhs=gri.gfaccess(\"in_gfs\",\"cf\"),rhs=IDcf)]\n",
    "\n",
    "outCparams = \"preindent=1,outCfileaccess=a,outCverbose=False,includebraces=False\"\n",
    "all_but_lambdaU_outC = fin.FD_outputC(\"returnstring\",all_but_lambdaU_expressions, outCparams)\n",
    "outstr += lp.loop([\"i2\",\"i1\",\"i0\"],[\"0\",\"0\",\"0\"],[\"cctk_lsh[2]\",\"cctk_lsh[1]\",\"cctk_lsh[0]\"],\n",
    "                   [\"1\",\"1\",\"1\"],[\"#pragma omp parallel for\",\"\",\"\"],\"    \",all_but_lambdaU_outC)\n",
    "\n",
    "outstr += \"} // END void WeylScal4NRPD_ADM_to_BSSN(CCTK_ARGUMENTS)\\n\"\n",
    "\n",
    "with open(\"WeylScal4NRPD/src/ADM_to_BSSN.c\", \"w\") as file:\n",
    "    file.write(str(outstr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cclfiles'></a>\n",
    "\n",
    "## Step 2.b: CCL files - Define how this module interacts and interfaces with the larger Einstein Toolkit infrastructure \\[Back to [top](#toc)\\]\n",
    "$$\\label{cclfiles}$$\n",
    "\n",
    "Writing a module (\"thorn\") within the Einstein Toolkit requires that three \"ccl\" files be constructed, all in the root directory of the thorn:\n",
    "\n",
    "1.`interface.ccl`: defines the gridfunction groups needed, and provides keywords denoting what this thorn provides and what it should inherit from other thorns.\n",
    "1. `param.ccl`: specifies free parameters within the thorn.\n",
    "1. `schedule.ccl`: allocates storage for gridfunctions, defines how the thorn's functions should be scheduled in a broader simulation, and specifies the regions of memory written to or read from gridfunctions.\n",
    "\n",
    "Let's start with `interface.ccl`. The [official Einstein Toolkit (Cactus) documentation](http://einsteintoolkit.org/usersguide/UsersGuide.html) defines what must/should be included in an `interface.ccl` file [**here**](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-178000D2.2). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting WeylScal4NRPD/interface.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile WeylScal4NRPD/interface.ccl\n",
    "\n",
    "# With \"implements\", we give our thorn its unique name.\n",
    "implements: WeylScal4NRPD\n",
    "\n",
    "# By \"inheriting\" other thorns, we tell the Toolkit that we \n",
    "#   will rely on variables/function that exist within those\n",
    "#   functions. \n",
    "inherits:   admbase Boundary Grid methodoflines\n",
    "\n",
    "# Tell the Toolkit that we want the various Weyl scalars \n",
    "#    and invariants to be visible to other thorns by using \n",
    "#    the keyword \"public\". Note that declaring these \n",
    "#    gridfunctions *does not* allocate memory for them;\n",
    "#    that is done by the schedule.ccl file.\n",
    "public:\n",
    "CCTK_REAL NRPyPsi4_group type=GF timelevels=3 tags='tensortypealias=\"Scalar\" tensorweight=0 tensorparity=1'\n",
    "{\n",
    "  psi4rGF,psi4r0ptGF,psi4r1ptGF,psi4r2ptGF, psi4iGF\n",
    "} \"Psi4_group\"\n",
    "\n",
    "CCTK_REAL evol_variables type = GF Timelevels=3\n",
    "{\n",
    "    aDD00GF,aDD01GF,aDD02GF,aDD11GF,aDD12GF,aDD22GF,alphaGF,cfGF,hDD00GF,hDD01GF,hDD02GF,hDD11GF,hDD12GF,hDD22GF,trKGF,vetU0GF,vetU1GF,vetU2GF\n",
    "} \"BSSN evolved gridfunctions, sans lambdaU and partial t beta\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now write the file `param.ccl`. This file allows the listed parameters to be set at runtime. We also give allowed ranges and default values for each parameter. More information on this file's syntax can be found in the [official Einstein Toolkit documentation](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-183000D2.3). \n",
    "\n",
    "The first parameter specifies how many time levels need to be stored. Generally when using the ETK's adaptive-mesh refinement (AMR) driver [Carpet](https://carpetcode.org/), three timelevels are needed so that the diagnostic quantities can be properly interpolated and defined across refinement boundaries. \n",
    "\n",
    "The second parameter determines how often we will calculate $\\psi_4$, and the third parameter indicates whether just $\\psi_4$, all Weyl scalars, or all Weyl scalars and invariants are going to be output. The third parameter is currently specified entirely within NRPy+, so by this point, it is *not* a free parameter. Thus it is not quite correct to include it in this list of *free* parameters (FIXME)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting WeylScal4NRPD/param.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile WeylScal4NRPD/param.ccl\n",
    "\n",
    "restricted:\n",
    "CCTK_INT timelevels \"Number of active timelevels\" STEERABLE=RECOVER\n",
    "{\n",
    "  0:3 :: \"\"\n",
    "} 3\n",
    "\n",
    "restricted:\n",
    "CCTK_INT WeylScal4NRPD_calc_every \"WeylScal4_psi4_calc_Nth_calc_every\" STEERABLE=ALWAYS\n",
    "{\n",
    "  *:* :: \"\"\n",
    "} 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we will write the file `schedule.ccl`; its official documentation is found [here](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-186000D2.4). This file dictates when the various parts of the thorn will be run. We first assign storage for both the real and imaginary components of $\\psi_4$, and then specify that we want our code run in  the `MoL_PseudoEvolution` schedule group (consistent with the original `WeylScal4` Einstein Toolkit thorn), after the ADM variables are set. At this step, we declare that we will be writing code in C. We also specify the gridfunctions that we wish to read in from memory--in our case, we need all the components of $K_{ij}$ (the spatial extrinsic curvature) and $\\gamma_{ij}$ (the physical [as opposed to conformal] 3-metric), in addition to the coordinate values. Note that the ETK adopts the widely-used convention that components of $\\gamma_{ij}$ are prefixed in the code with $\\text{g}$ and not $\\gamma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting WeylScal4NRPD/schedule.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile WeylScal4NRPD/schedule.ccl\n",
    "\n",
    "STORAGE: NRPyPsi4_group[3], evol_variables[3]\n",
    "STORAGE: ADMBase::metric[3], ADMBase::curv[3], ADMBase::lapse[3], ADMBase::shift[3]\n",
    "\n",
    "schedule group WeylScal4NRPD_group in MoL_PseudoEvolution after ADMBase_SetADMVars\n",
    "{\n",
    "} \"Schedule WeylScal4NRPD group\"\n",
    "\n",
    "schedule WeylScal4NRPD_ADM_to_BSSN in WeylScal4NRPD_group before weylscal4_mainfunction\n",
    "{\n",
    "  LANG: C\n",
    "} \"Convert ADM into BSSN variables\"\n",
    "\n",
    "\n",
    "schedule WeylScal4NRPD_mainfunction in WeylScal4NRPD_group after WeylScal4NRPD_ADM_to_BSSN\n",
    "{\n",
    "  LANG: C\n",
    "} \"Call WeylScal4NRPD main function\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='etk_list'></a>\n",
    "\n",
    "## Step 2.c: Tell the Einstein Toolkit to compile the C code \\[Back to [top](#toc)\\]\n",
    "$$\\label{etk_list}$$\n",
    "\n",
    "The `make.code.defn` lists the source files that need to be compiled. Naturally, this thorn has only the one C file $-$ written above $-$ to compile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting WeylScal4NRPD/src/make.code.defn\n"
     ]
    }
   ],
   "source": [
    "%%writefile WeylScal4NRPD/src/make.code.defn\n",
    "\n",
    "SRCS = WeylScal4NRPD.c ADM_to_BSSN.c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 3: 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-ETK_thorn-Weyl_Scalars_and_Spacetime_Invariants.pdf](Tutorial-ETK_thorn-Weyl_Scalars_and_Spacetime_Invariants.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": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ETK_thorn-WeylScal4NRPD.tex, and compiled LaTeX file to PDF file Tutorial-ETK_thorn-WeylScal4NRPD.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ETK_thorn-WeylScal4NRPD\")"
   ]
  }
 ],
 "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
}