Tutorial-ETK_thorn-Weyl_Scalars_and_Spacetime_Invariants.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",
    "# 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",
    "[comment]: <> (Abstract: TODO)\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",
    "* [WeylScal4NRPy/WeylScalars_Cartesian.py](../edit/WeylScal4NRPy/WeylScalars_Cartesian.py)\n",
    "* [WeylScal4NRPy/WeylScalarInvariants_Cartesian.py](../edit/WeylScal4NRPy/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 [WeylScal4NRPy/WeylScalars_Cartesian.py](../edit/WeylScal4NRPy/WeylScalars_Cartesian.py), and the Weyl invariant SymPy expressions to [WeylScal4NRPy/WeylScalarInvariants_Cartesian.py](../edit/WeylScal4NRPy/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": {
    "execution": {
     "iopub.execute_input": "2021-08-28T19:47:16.433188Z",
     "iopub.status.busy": "2021-08-28T19:47:16.423717Z",
     "iopub.status.idle": "2021-08-28T19:48:09.460141Z",
     "shell.execute_reply": "2021-08-28T19:48:09.459069Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function calc_psis() to file WeylScal4NRPy/src/calc_psis.h\n",
      "Output C function calc_invars() to file WeylScal4NRPy/src/calc_invars.h\n"
     ]
    }
   ],
   "source": [
    "from outputC import lhrh, outCfunction  # 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 os                        # Standard Python modules for multiplatform OS-level functions\n",
    "\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",
    "\n",
    "import WeylScal4NRPy.WeylScalars_Cartesian as weyl\n",
    "par.set_parval_from_str(\"output_scalars\",\"all_psis_and_invariants\")\n",
    "\n",
    "weyl.WeylScalars_Cartesian()\n",
    "\n",
    "output_scalars = par.parval_from_str(\"output_scalars\")\n",
    "\n",
    "!mkdir WeylScal4NRPy     2>/dev/null # 2>/dev/null: Don't throw an error or warning if the directory already exists.\n",
    "!mkdir WeylScal4NRPy/src 2>/dev/null # 2>/dev/null: Don't throw an error or warning if the directory already exists.\n",
    "\n",
    "scalars_lhrh = [lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi4r\"),rhs=weyl.psi4r),\n",
    "                lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi4i\"),rhs=weyl.psi4i)]\n",
    "\n",
    "if output_scalars in ('all_psis', 'all_psis_and_invariants'):\n",
    "    scalars_lhrh = [\n",
    "                    lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi4r\"),rhs=weyl.psi4r),\n",
    "                    lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi4i\"),rhs=weyl.psi4i),\n",
    "                    lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi3r\"),rhs=weyl.psi3r),\n",
    "                    lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi3i\"),rhs=weyl.psi3i),\n",
    "                    lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi2r\"),rhs=weyl.psi2r),\n",
    "                    lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi2i\"),rhs=weyl.psi2i),\n",
    "                    lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi1r\"),rhs=weyl.psi1r),\n",
    "                    lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi1i\"),rhs=weyl.psi1i),\n",
    "                    lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi0r\"),rhs=weyl.psi0r),\n",
    "                    lhrh(lhs=gri.gfaccess(\"out_gfs\",\"psi0i\"),rhs=weyl.psi0i)\n",
    "                   ]\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",
    "psis_CcodeKernel = fin.FD_outputC(\"returnstring\",scalars_lhrh,params=\"outCverbose=False,CSE_sorting=none\")\n",
    "\n",
    "cctk_params_readin_str = r\"\"\"\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  if(cctk_nghostzones[0] != cctk_nghostzones[1] ||\n",
    "     cctk_nghostzones[0] != cctk_nghostzones[2]) {\n",
    "    CCTK_ERROR(\"cctk_nghostzones[i] must be the same in all directions i\");\n",
    "  }\n",
    "  const CCTK_INT NGHOSTS CCTK_ATTRIBUTE_UNUSED = cctk_nghostzones[0];\n",
    "  const CCTK_INT Nxx0 CCTK_ATTRIBUTE_UNUSED = cctk_lsh[0] - 2*NGHOSTS;\n",
    "  const CCTK_INT Nxx1 CCTK_ATTRIBUTE_UNUSED = cctk_lsh[1] - 2*NGHOSTS;\n",
    "  const CCTK_INT Nxx2 CCTK_ATTRIBUTE_UNUSED = cctk_lsh[2] - 2*NGHOSTS;\n",
    "  const CCTK_INT Nxx_plus_2NGHOSTS0 CCTK_ATTRIBUTE_UNUSED = cctk_lsh[0];\n",
    "  const CCTK_INT Nxx_plus_2NGHOSTS1 CCTK_ATTRIBUTE_UNUSED = cctk_lsh[1];\n",
    "  const CCTK_INT Nxx_plus_2NGHOSTS2 CCTK_ATTRIBUTE_UNUSED = cctk_lsh[2];\n",
    "\"\"\"\n",
    "\n",
    "desc = \"Calculate the Weyl Scalars\"\n",
    "name = \"calc_psis\"\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(\"WeylScal4NRPy\",\"src\",name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"\"\"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 *gammaDD00GF,const CCTK_REAL *gammaDD01GF,const CCTK_REAL *gammaDD02GF,const CCTK_REAL *gammaDD11GF,const CCTK_REAL *gammaDD12GF,const CCTK_REAL *gammaDD22GF,\n",
    "               const CCTK_REAL     *kDD00GF,const CCTK_REAL     *kDD01GF,const CCTK_REAL     *kDD02GF,const CCTK_REAL     *kDD11GF,const CCTK_REAL     *kDD12GF,const CCTK_REAL     *kDD22GF,\n",
    "               CCTK_REAL *psi4rGF,CCTK_REAL *psi4iGF,\n",
    "               CCTK_REAL *psi3rGF,CCTK_REAL *psi3iGF,\n",
    "               CCTK_REAL *psi2rGF,CCTK_REAL *psi2iGF,\n",
    "               CCTK_REAL *psi1rGF,CCTK_REAL *psi1iGF,\n",
    "               CCTK_REAL *psi0rGF,CCTK_REAL *psi0iGF\"\"\",\n",
    "    preloop=cctk_params_readin_str,\n",
    "    body     = psis_CcodeKernel,\n",
    "    loopopts =\"InteriorPoints\",\n",
    "    enableCparameters=False)\n",
    "\n",
    "# Reset the registered gridfunctions list.\n",
    "gri.glb_gridfcs_list = []\n",
    "#par.set_parval_from_str(\"WeylScal4NRPy.WeylScalars_Cartesian::output_scalars\",\"all_psis_and_invariants\")\n",
    "output_scalars = par.parval_from_str(\"output_scalars\")\n",
    "\n",
    "import WeylScal4NRPy.WeylScalarInvariants_Cartesian as invar\n",
    "invar.WeylScalarInvariants_Cartesian()\n",
    "invars_lhrh = [\n",
    "               lhrh(lhs=gri.gfaccess(\"out_gfs\",\"curvIr\"),rhs=invar.curvIr),\n",
    "               lhrh(lhs=gri.gfaccess(\"out_gfs\",\"curvIi\"),rhs=invar.curvIi),\n",
    "               lhrh(lhs=gri.gfaccess(\"out_gfs\",\"curvJr\"),rhs=invar.curvJr),\n",
    "               lhrh(lhs=gri.gfaccess(\"out_gfs\",\"curvJi\"),rhs=invar.curvJi),\n",
    "               lhrh(lhs=gri.gfaccess(\"out_gfs\",\"J1curv\"),rhs=invar.J1curv),\n",
    "               lhrh(lhs=gri.gfaccess(\"out_gfs\",\"J2curv\"),rhs=invar.J2curv),\n",
    "               lhrh(lhs=gri.gfaccess(\"out_gfs\",\"J3curv\"),rhs=invar.J3curv),\n",
    "               lhrh(lhs=gri.gfaccess(\"out_gfs\",\"J4curv\"),rhs=invar.J4curv)\n",
    "              ]\n",
    "\n",
    "invars_CcodeKernel = fin.FD_outputC(\"returnstring\",invars_lhrh,\n",
    "                    params=\"outCverbose=False,CSE_sorting=none\")# 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",
    "\n",
    "desc = \"Calculate the Weyl Invariants\"\n",
    "name = \"calc_invars\"\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(\"WeylScal4NRPy\",\"src\",name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"\"\"const cGH* restrict const cctkGH,const int *cctk_lsh,const int *cctk_nghostzones,\n",
    "                 const CCTK_REAL *psi4rGF,const CCTK_REAL *psi4iGF,\n",
    "                 const CCTK_REAL *psi3rGF,const CCTK_REAL *psi3iGF,\n",
    "                 const CCTK_REAL *psi2rGF,const CCTK_REAL *psi2iGF,\n",
    "                 const CCTK_REAL *psi1rGF,const CCTK_REAL *psi1iGF,\n",
    "                 const CCTK_REAL *psi0rGF,const CCTK_REAL *psi0iGF,\n",
    "                 CCTK_REAL *curvIrGF,CCTK_REAL *curvIiGF,\n",
    "                 CCTK_REAL *curvJrGF,CCTK_REAL *curvJiGF,\n",
    "                 CCTK_REAL *J1curvGF,CCTK_REAL *J2curvGF,\n",
    "                 CCTK_REAL *J3curvGF,CCTK_REAL *J4curvGF\"\"\",\n",
    "    preloop=cctk_params_readin_str,\n",
    "    body     = invars_CcodeKernel,\n",
    "    loopopts =\"InteriorPoints\",\n",
    "    enableCparameters=False)\n"
   ]
  },
  {
   "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 (`WeylScal4NRPy_psis.h` and `WeylScal4NRPy_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": {
    "execution": {
     "iopub.execute_input": "2021-08-28T19:48:09.466392Z",
     "iopub.status.busy": "2021-08-28T19:48:09.465470Z",
     "iopub.status.idle": "2021-08-28T19:48:09.468436Z",
     "shell.execute_reply": "2021-08-28T19:48:09.469130Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting WeylScal4NRPy/src/WeylScal4NRPy.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile WeylScal4NRPy/src/WeylScal4NRPy.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",
    "#include \"calc_psis.h\"\n",
    "#include \"calc_invars.h\"\n",
    "\n",
    "extern void weylscal4_mainfunction(CCTK_ARGUMENTS) {\n",
    "\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "\n",
    "  if(cctk_iteration % WeylScal4NRPy_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",
    "  calc_psis(cctkGH,cctk_lsh,cctk_nghostzones,\n",
    "            invdx0,invdx1,invdx2,\n",
    "            x,y,z,\n",
    "            gxx,gxy,gxz,gyy,gyz,gzz,\n",
    "            kxx,kxy,kxz,kyy,kyz,kzz,\n",
    "            psi4r,psi4i,\n",
    "            psi3r,psi3i,\n",
    "            psi2r,psi2i,\n",
    "            psi1r,psi1i,\n",
    "            psi0r,psi0i);\n",
    "\n",
    "  if (CCTK_EQUALS(output_scalars, \"all_psis_and_invariants\")) {\n",
    "    calc_invars(cctkGH,cctk_lsh,cctk_nghostzones,\n",
    "      \t        psi4r,psi4i,\n",
    "                psi3r,psi3i,\n",
    "                psi2r,psi2i,\n",
    "                psi1r,psi1i,\n",
    "                psi0r,psi0i,\n",
    "                NRPycurvIr,NRPycurvIi,\n",
    "                NRPycurvJr,NRPycurvJi,\n",
    "                NRPyJ1curv,NRPyJ2curv,\n",
    "                NRPyJ3curv,NRPyJ4curv);\n",
    "    }\n",
    "}\n"
   ]
  },
  {
   "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": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-08-28T19:48:09.474328Z",
     "iopub.status.busy": "2021-08-28T19:48:09.473323Z",
     "iopub.status.idle": "2021-08-28T19:48:09.476413Z",
     "shell.execute_reply": "2021-08-28T19:48:09.477174Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting WeylScal4NRPy/interface.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile WeylScal4NRPy/interface.ccl\n",
    "\n",
    "# With \"implements\", we give our thorn its unique name.\n",
    "implements: WeylScal4NRPy\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",
    "  psi4r, psi4i\n",
    "} \"Psi4_group\"\n",
    "\n",
    "public:\n",
    "CCTK_REAL NRPyPsi3210_group type=GF timelevels=3 tags='tensortypealias=\"Scalar\" tensorweight=0 tensorparity=1'\n",
    "{\n",
    "  psi3r,psi3i,psi2r,psi2i,psi1r,psi1i,psi0r,psi0i\n",
    "} \"Psi3210_group\"\n",
    "\n",
    "public:\n",
    "CCTK_REAL NRPyInvars_group type=GF timelevels=3 tags='tensortypealias=\"Scalar\" tensorweight=0 tensorparity=1'\n",
    "{\n",
    "  NRPycurvIr,NRPycurvIi,NRPycurvJr,NRPycurvJi,NRPyJ1curv,NRPyJ2curv,NRPyJ3curv,NRPyJ4curv\n",
    "} \"NRPyInvars_group\""
   ]
  },
  {
   "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": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-08-28T19:48:09.483515Z",
     "iopub.status.busy": "2021-08-28T19:48:09.482807Z",
     "iopub.status.idle": "2021-08-28T19:48:09.485704Z",
     "shell.execute_reply": "2021-08-28T19:48:09.485189Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting WeylScal4NRPy/param.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile WeylScal4NRPy/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 WeylScal4NRPy_calc_every \"WeylScal4_psi4_calc_Nth_calc_every\" STEERABLE=ALWAYS\n",
    "{\n",
    "  *:* :: \"\"\n",
    "} 1\n",
    "\n",
    "private:\n",
    "CCTK_KEYWORD output_scalars \"Whether to output all Weyl scalars, just psi4, or all scalars and invariants\"\n",
    "{\n",
    "  \"all_psis\" :: \"\"\n",
    "  \"all_psis_and_invariants\" :: \"\"\n",
    "} \"all_psis\"\n"
   ]
  },
  {
   "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": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-08-28T19:48:09.490837Z",
     "iopub.status.busy": "2021-08-28T19:48:09.489933Z",
     "iopub.status.idle": "2021-08-28T19:48:09.492867Z",
     "shell.execute_reply": "2021-08-28T19:48:09.493554Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting WeylScal4NRPy/schedule.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile WeylScal4NRPy/schedule.ccl\n",
    "\n",
    "STORAGE: NRPyPsi4_group[timelevels]\n",
    "if (CCTK_EQUALS(output_scalars, \"all_psis_and_invariants\") || CCTK_EQUALS(output_scalars, \"all_psis\"))\n",
    "{\n",
    " STORAGE: NRPyPsi3210_group[timelevels]\n",
    "}\n",
    "if (CCTK_EQUALS(output_scalars, \"all_psis_and_invariants\"))\n",
    "{\n",
    " STORAGE: NRPyInvars_group[timelevels]\n",
    "}\n",
    "\n",
    "schedule group WeylScal4NRPy_group in MoL_PseudoEvolution after ADMBase_SetADMVars\n",
    "{\n",
    "} \"Schedule WeylScal4NRPy group\"\n",
    "\n",
    "schedule weylscal4_mainfunction in WeylScal4NRPy_group\n",
    "{\n",
    "  LANG: C\n",
    "   READS: admbase::kxx(Everywhere)\n",
    "   READS: admbase::kxy(Everywhere)\n",
    "   READS: admbase::kxz(Everywhere)\n",
    "   READS: admbase::kyy(Everywhere)\n",
    "   READS: admbase::kyz(Everywhere)\n",
    "   READS: admbase::kzz(Everywhere)\n",
    "   READS: admbase::gxx(Everywhere)\n",
    "   READS: admbase::gxy(Everywhere)\n",
    "   READS: admbase::gxz(Everywhere)\n",
    "   READS: admbase::gyy(Everywhere)\n",
    "   READS: admbase::gyz(Everywhere)\n",
    "   READS: admbase::gzz(Everywhere)\n",
    "   READS: grid::x(Everywhere)\n",
    "   READS: grid::y(Everywhere)\n",
    "   READS: grid::z(Everywhere)\n",
    "   WRITES: WeylScal4NRPy::psi4i(Interior)\n",
    "   WRITES: WeylScal4NRPy::psi4r(Interior)\n",
    "   WRITES: WeylScal4NRPy::psi3i(Interior)\n",
    "   WRITES: WeylScal4NRPy::psi3r(Interior)\n",
    "   WRITES: WeylScal4NRPy::psi2i(Interior)\n",
    "   WRITES: WeylScal4NRPy::psi2r(Interior)\n",
    "   WRITES: WeylScal4NRPy::psi1i(Interior)\n",
    "   WRITES: WeylScal4NRPy::psi1r(Interior)\n",
    "   WRITES: WeylScal4NRPy::psi0i(Interior)\n",
    "   WRITES: WeylScal4NRPy::psi0r(Interior)\n",
    "   WRITES: WeylScal4NRPy::NRPycurvIi(Interior)\n",
    "   WRITES: WeylScal4NRPy::NRPycurvIr(Interior)\n",
    "   WRITES: WeylScal4NRPy::NRPycurvJi(Interior)\n",
    "   WRITES: WeylScal4NRPy::NRPycurvJr(Interior)\n",
    "   WRITES: WeylScal4NRPy::NRPyJ1curv(Interior)\n",
    "   WRITES: WeylScal4NRPy::NRPyJ2curv(Interior)\n",
    "   WRITES: WeylScal4NRPy::NRPyJ3curv(Interior)\n",
    "   WRITES: WeylScal4NRPy::NRPyJ4curv(Interior)\n",
    "\n",
    "} \"Call WeylScal4NRPy main function\"\n"
   ]
  },
  {
   "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": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-08-28T19:48:09.497693Z",
     "iopub.status.busy": "2021-08-28T19:48:09.496797Z",
     "iopub.status.idle": "2021-08-28T19:48:09.499842Z",
     "shell.execute_reply": "2021-08-28T19:48:09.500556Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting WeylScal4NRPy/src/make.code.defn\n"
     ]
    }
   ],
   "source": [
    "%%writefile WeylScal4NRPy/src/make.code.defn\n",
    "\n",
    "SRCS = WeylScal4NRPy.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": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-08-28T19:48:09.505999Z",
     "iopub.status.busy": "2021-08-28T19:48:09.505245Z",
     "iopub.status.idle": "2021-08-28T19:48:13.945665Z",
     "shell.execute_reply": "2021-08-28T19:48:13.946209Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ETK_thorn-Weyl_Scalars_and_Spacetime_Invariants.tex, and\n",
      "    compiled LaTeX file to PDF file Tutorial-ETK_thorn-\n",
      "    Weyl_Scalars_and_Spacetime_Invariants.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-Weyl_Scalars_and_Spacetime_Invariants\")"
   ]
  }
 ],
 "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": 2
}