Tutorial-ETK_thorn-u0_smallb_Poynting.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",
    "# sbPoynETNRPy: An Einstein Toolkit Thorn for Computing the 4-Velocity Time-Component $u^0$, the Magnetic Field Measured by a Comoving Observer $b^{\\mu}$, and the Poynting Vector $S^i$\n",
    "\n",
    "## Author: 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:** This module has been validated against the hand-written smallbPoynET in WVUThorns_diagnostics (a trsuted code), which itself is based on expressions in IllinoisGRMHD... which was validated against the original GRMHD code of the Illinois NR group.\n",
    "\n",
    "## Introduction:\n",
    "In the [previous tutorial notebook](Tutorial-u0_smallb_Poynting-Cartesian.ipynb), we constructed within SymPy full expressions for the 4-velocity time-component $u^0$, the magnetic field (measured by a comoving observer) $b^{\\mu}$, and the Poynting vector $S^i$.\n",
    "\n",
    "Here we will work through the steps necessary to construct an Einstein Toolkit diagnostic thorn (module) that uses ADMBase and HydroBase variables as input into the NRPy+-generated SymPy expressions for $b^{\\mu}$, $b^2$, and the Poynting Vector $S^i$, outputting to gridfunctions `smallb4U[]`, `smallb2etk` (the \"etk\" suffix must be appended because base gridfunction names ending in numbers are not allowed in NRPy+), and `SPoyn[]`, respectively. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This tutorial is organized as follows\n",
    "\n",
    "1. [Step 1](#initializenrpy): Call on NRPy+ to convert the SymPy expressions for $b^{\\mu}$, $b^2$, and the Poynting Vector $S^i$ to C code kernels\n",
    "1. [Step 2](#etk): Build up the needed Einstein Toolkit infrastructure to implement the NRPy+-generated C code kernels\n",
    "    1. [Step 2.a](#etkc): Write the C code functions called by the Einstein Toolkit scheduler that incorporate the above \".h\" files\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](#etksys): Inform the Einstein Toolkit build system of the C code\n",
    "1. [Step 3](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a> \n",
    "\n",
    "# Step 1: Call on NRPy+ to convert the SymPy expressions for $b^{\\mu}$, $b^2$, and the Poynting Vector $S^i$ to C code kernels \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:08:59.683166Z",
     "iopub.status.busy": "2021-03-07T17:08:59.682158Z",
     "iopub.status.idle": "2021-03-07T17:09:00.015095Z",
     "shell.execute_reply": "2021-03-07T17:09:00.014351Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1a: import all needed modules from NRPy+:\n",
    "import indexedexp as ixp                   # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "from outputC import outputC                # NRPy+: Basic C code output functionality\n",
    "import grid as gri                         # NRPy+: Functions having to do with numerical grids\n",
    "import NRPy_param_funcs as par             # NRPy+: parameter interface\n",
    "import sympy as sp                         # SymPy: The Python computer algebra package upon which NRPy+ depends"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:09:00.019462Z",
     "iopub.status.busy": "2021-03-07T17:09:00.018541Z",
     "iopub.status.idle": "2021-03-07T17:09:00.020794Z",
     "shell.execute_reply": "2021-03-07T17:09:00.021342Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1b: Initialize  parameters (stub; there are none for this module)\n",
    "thismodule = __name__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will to disable verbose output in the NRPy+ outputC function. This is an important step in this case, because our final expressions are very large. Verbose output, when enabled, will print (in comments) the input SymPy expressions to the top of the file *without* CSE, resulting here in an *enormous* output file.\n",
    "\n",
    "We will also declare the additional gridfunctions we need for this thorn:\n",
    "\n",
    "**Inputs from ADMBase:**\n",
    "* the physical metric $\\gamma_{ij}$\n",
    "* the spacetime gauge quantities $\\alpha$ and $\\beta^i$\n",
    "\n",
    "**Inputs from HydroBase:**\n",
    "* the Valencia 3-velocity $v^i_{(n)}$\n",
    "* the densitized magnetic field of a normal observer $\\tilde{B}^i$\n",
    "\n",
    "**Output gridfunctions:**\n",
    "* the magnetic field as observed in a frame comoving with the plasma $b^\\mu$ (`smallb4U[]}`)\n",
    "* twice the magnetic pressure $2 P_{\\rm mag} = b_\\mu b^\\mu = b^2$ (`smallb2etk`)\n",
    "* the Poynting vector $S^i$ (`SPoyn[]`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:09:00.048005Z",
     "iopub.status.busy": "2021-03-07T17:09:00.039618Z",
     "iopub.status.idle": "2021-03-07T17:09:00.768019Z",
     "shell.execute_reply": "2021-03-07T17:09:00.768466Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Wrote to file \"sbPoynETNRPy/src/u0.h\"\n",
      "Wrote to file \"sbPoynETNRPy/src/smallb4U_smallb2etk_PoynSU.h\"\n"
     ]
    }
   ],
   "source": [
    "# Step 1c: Set spatial dimension (must be 3 for BSSN)\n",
    "DIM = 3\n",
    "par.set_parval_from_str(\"grid::DIM\",DIM)\n",
    "\n",
    "# Step 1d: declare the additional gridfunctions (i.e., functions whose values are declared\n",
    "#          at every grid point, either inside or outside of our SymPy expressions) needed\n",
    "#         for this thorn\n",
    "\n",
    "# INPUT GRIDFUNCTIONS:\n",
    "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUX\",\"gammaDD\", \"sym01\") # The AUX or EVOL designation is *not*\n",
    "                                                                                # used in diagnostic modules.\n",
    "betaU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"betaU\") # The AUX or EVOL designation is *not*\n",
    "                                                                   # used in diagnostic modules.\n",
    "alpha = gri.register_gridfunctions(\"AUX\",\"alpha\") # The AUX or EVOL designation is *not*\n",
    "                                                  # used in diagnostic modules.\n",
    "ValenciavU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"ValenciavU\") # The AUX or EVOL designation is *not*\n",
    "                                                                             # used in diagnostic modules.\n",
    "BU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"BU\") # The AUX or EVOL designation is *not*\n",
    "                                                             # used in diagnostic modules.\n",
    "\n",
    "# OUTPUT GRIDFUNCTIONS:\n",
    "smallb4U = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"smallb4U\",DIM=4) # The AUX or EVOL designation is *not*\n",
    "                                                                               # used in diagnostic modules.\n",
    "smallb2etk = gri.register_gridfunctions(\"AUX\",\"smallb2etk\") # The AUX or EVOL designation is *not*\n",
    "                                                            # used in diagnostic modules.\n",
    "PoynSU = ixp.register_gridfunctions_for_single_rank1(\"AUX\",\"PoynSU\") # The AUX or EVOL designation is *not*\n",
    "                                                                     # used in diagnostic modules.\n",
    "\n",
    "# Step 1f: Call the NRPy+ module to set up the SymPy expressions for the output, as well as the C code for computing u^0\n",
    "import u0_smallb_Poynting__Cartesian.u0_smallb_Poynting__Cartesian as u0etc\n",
    "u0etc.compute_u0_smallb_Poynting__Cartesian(gammaDD,betaU,alpha,ValenciavU,BU)\n",
    "\n",
    "# Step 1g: Set the gridfunction memory access type to \"ETK\":\n",
    "par.set_parval_from_str(\"GridFuncMemAccess\",\"ETK\")\n",
    "\n",
    "# Step 1h: Make output directories:\n",
    "!mkdir sbPoynETNRPy     2>/dev/null # 2>/dev/null: Don't throw an error or warning if the directory already exists.\n",
    "!mkdir sbPoynETNRPy/src 2>/dev/null # 2>/dev/null: Don't throw an error or warning if the directory already exists.\n",
    "\n",
    "# Step 1i: Output routine for computing u0:\n",
    "with open(\"sbPoynETNRPy/src/u0.h\", \"w\") as file:\n",
    "    file.write(str(u0etc.computeu0_Cfunction))\n",
    "    print(\"Wrote to file \\\"\"+file.name+\"\\\"\")\n",
    "\n",
    "# Step 1j: Use NRPy+'s outputC to convert the SymPy expressions for smallb4U, smallb2etk, and PoynSU to C code:\n",
    "#outputC([u0etc.smallb4U[0],u0etc.smallb4U[1],u0etc.smallb4U[2],u0etc.smallb4U[3],u0etc.smallb2etk,\n",
    "outputC([u0etc.smallb4U[0],u0etc.smallb4U[1],u0etc.smallb4U[2],u0etc.smallb4U[3],u0etc.smallb2etk,\n",
    "         u0etc.PoynSU[0],u0etc.PoynSU[1],u0etc.PoynSU[2]],\n",
    "         [gri.gfaccess(\"\",\"smallb4U0\"),gri.gfaccess(\"\",\"smallb4U1\"),gri.gfaccess(\"\",\"smallb4U2\"),gri.gfaccess(\"\",\"smallb4U3\"),\n",
    "          gri.gfaccess(\"\",\"smallb2etk\"),\n",
    "          gri.gfaccess(\"\",\"PoynSU0\"),gri.gfaccess(\"\",\"PoynSU1\"),gri.gfaccess(\"\",\"PoynSU2\")],\n",
    "       filename=\"sbPoynETNRPy/src/smallb4U_smallb2etk_PoynSU.h\",\n",
    "       params=\"outCverbose=False\") # <- Force outCverbose=False for this\n",
    "                                   #    module to avoid gigantic C file filled with the\n",
    "                                   #    non-CSE expressions for the Weyl scalars."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='etk'></a>\n",
    "\n",
    "# Step 2: Build up the needed Einstein Toolkit infrastructure to implement the NRPy+-generated C code kernels \\[Back to [top](#toc)\\]\n",
    "$$\\label{etk}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='etkc'></a> \n",
    "\n",
    "## Step 2.a: Write the C code functions called by the Einstein Toolkit scheduler that incorporate the above \".h\" files \\[Back to [top](#toc)\\]\n",
    "$$\\label{etkc}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:09:00.775209Z",
     "iopub.status.busy": "2021-03-07T17:09:00.774341Z",
     "iopub.status.idle": "2021-03-07T17:09:00.777389Z",
     "shell.execute_reply": "2021-03-07T17:09:00.777830Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting sbPoynETNRPy/src/sbPoynETNRPy.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile sbPoynETNRPy/src/sbPoynETNRPy.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 sbPoynETNRPy_lowlevel(const cGH* restrict const cctkGH,const int *cctk_lsh,\n",
    "                           const CCTK_REAL *gammaDD00GF,const CCTK_REAL *gammaDD01GF,const CCTK_REAL *gammaDD02GF,\n",
    "                           const CCTK_REAL *gammaDD11GF,const CCTK_REAL *gammaDD12GF,const CCTK_REAL *gammaDD22GF,\n",
    "                           const CCTK_REAL *alphaGF,\n",
    "                           const CCTK_REAL *betaU0GF,const CCTK_REAL *betaU1GF,const CCTK_REAL *betaU2GF,\n",
    "                           const CCTK_REAL *vel,const CCTK_REAL *Bvec,\n",
    "\n",
    "                           CCTK_REAL *smallb4U0GF,CCTK_REAL *smallb4U1GF,CCTK_REAL *smallb4U2GF,CCTK_REAL *smallb4U3GF,\n",
    "                           CCTK_REAL *smallb2etkGF,\n",
    "                           CCTK_REAL *PoynSU0GF,CCTK_REAL *PoynSU1GF,CCTK_REAL *PoynSU2GF) {\n",
    "\n",
    "    DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "#pragma omp parallel for\n",
    "    for(int i2=0;i2<cctk_lsh[2];i2++) for(int i1=0;i1<cctk_lsh[1];i1++) for(int i0=0;i0<cctk_lsh[0];i0++) {\n",
    "        const CCTK_REAL gammaDD00 = gammaDD00GF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "        const CCTK_REAL gammaDD01 = gammaDD01GF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "        const CCTK_REAL gammaDD02 = gammaDD02GF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "        const CCTK_REAL gammaDD11 = gammaDD11GF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "        const CCTK_REAL gammaDD12 = gammaDD12GF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "        const CCTK_REAL gammaDD22 = gammaDD22GF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "\n",
    "        const CCTK_REAL alpha = alphaGF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "\n",
    "        const CCTK_REAL betaU0 = betaU0GF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "        const CCTK_REAL betaU1 = betaU1GF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "        const CCTK_REAL betaU2 = betaU2GF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)];\n",
    "\n",
    "        // Valencia 3-velocity may be adjusted due to the velocity ceiling.\n",
    "        CCTK_REAL ValenciavU0 = vel[CCTK_GFINDEX4D(cctkGH, i0,i1,i2, 0)];\n",
    "        CCTK_REAL ValenciavU1 = vel[CCTK_GFINDEX4D(cctkGH, i0,i1,i2, 1)];\n",
    "        CCTK_REAL ValenciavU2 = vel[CCTK_GFINDEX4D(cctkGH, i0,i1,i2, 2)];\n",
    "\n",
    "        const CCTK_REAL BU0 = Bvec[CCTK_GFINDEX4D(cctkGH, i0,i1,i2, 0)];\n",
    "        const CCTK_REAL BU1 = Bvec[CCTK_GFINDEX4D(cctkGH, i0,i1,i2, 1)];\n",
    "        const CCTK_REAL BU2 = Bvec[CCTK_GFINDEX4D(cctkGH, i0,i1,i2, 2)];\n",
    "\n",
    "        CCTK_REAL u0;\n",
    "#include \"u0.h\"\n",
    "#include \"smallb4U_smallb2etk_PoynSU.h\"\n",
    "    }\n",
    "}\n",
    "\n",
    "extern void sbPoynETNRPy(CCTK_ARGUMENTS) {\n",
    "\n",
    "    DECLARE_CCTK_PARAMETERS;\n",
    "    DECLARE_CCTK_ARGUMENTS;\n",
    "\n",
    "    if(sbPoynETNRPy_calc_every<=0 || cctk_iteration%sbPoynETNRPy_calc_every!=0) { return; }\n",
    "\n",
    "    /* Calculate smallb4U[], smallb2etk, and PoynSU[]: */\n",
    "    sbPoynETNRPy_lowlevel(cctkGH,cctk_lsh,\n",
    "                          gxx,gxy,gxz,gyy,gyz,gzz,\n",
    "                          alp,\n",
    "                          betax,betay,betaz,\n",
    "                          vel,Bvec,\n",
    "\n",
    "                          smallb4U0,smallb4U1,smallb4U2,smallb4U3,\n",
    "                          smallb4_sq,\n",
    "                          PoynSU0,PoynSU1,PoynSU2);\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": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:09:00.783226Z",
     "iopub.status.busy": "2021-03-07T17:09:00.782283Z",
     "iopub.status.idle": "2021-03-07T17:09:00.785795Z",
     "shell.execute_reply": "2021-03-07T17:09:00.786293Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting sbPoynETNRPy/interface.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile sbPoynETNRPy/interface.ccl\n",
    "\n",
    "# With \"implements\", we give our thorn its unique name.\n",
    "implements: sbPoynETNRPy\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 HydroBase 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 smallb4U_group type=GF timelevels=3\n",
    "{\n",
    "  smallb4U0,smallb4U1,smallb4U2,smallb4U3\n",
    "} \"smallb4U 4-vector\"\n",
    "\n",
    "public:\n",
    "CCTK_REAL smallb4_sq_group type=GF timelevels=3\n",
    "{\n",
    "  smallb4_sq\n",
    "} \"smallb^{mu} squared == twice the magnetic pressure\"\n",
    "\n",
    "public:\n",
    "CCTK_REAL PoynSU_group type=GF timelevels=3\n",
    "{\n",
    "  PoynSU0,PoynSU1,PoynSU2\n",
    "} \"Poynting 3-vector\""
   ]
  },
  {
   "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 $b^\\mu$, $b^2$, and $S^i$.\n",
    "\n",
    "The third parameter sets the maximum allowed Lorentz factor when computing $u^0$ (i.e., $\\Gamma_{\\rm max}$, as defined in the [previous tutorial notebook](Tutorial-u0_smallb_Poynting-Cartesian.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:09:00.791466Z",
     "iopub.status.busy": "2021-03-07T17:09:00.790727Z",
     "iopub.status.idle": "2021-03-07T17:09:00.793796Z",
     "shell.execute_reply": "2021-03-07T17:09:00.794266Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting sbPoynETNRPy/param.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile sbPoynETNRPy/param.ccl\n",
    "\n",
    "shares: HydroBase\n",
    "USES CCTK_INT timelevels\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 sbPoynETNRPy_calc_every \"Compute these quantities every sbPoynETNRPy_calc_every iterations.\" STEERABLE=ALWAYS\n",
    "{\n",
    "  *:* :: \"\"\n",
    "} 1\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL GAMMA_SPEED_LIMIT \"Maximum Lorentz factor.\"\n",
    "{\n",
    " 1:* :: \"Positive > 1, though you'll likely have troubles in GRMHD far above 10, or far above 2000 in GRFFE.\"\n",
    "} 10.0\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). \n",
    "\n",
    "This file registers the function we wish to call, `sbPoynETNRPy`, with the Einstein Toolkit scheduler."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:09:00.800725Z",
     "iopub.status.busy": "2021-03-07T17:09:00.799663Z",
     "iopub.status.idle": "2021-03-07T17:09:00.804775Z",
     "shell.execute_reply": "2021-03-07T17:09:00.803732Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting sbPoynETNRPy/schedule.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile sbPoynETNRPy/schedule.ccl\n",
    "\n",
    "STORAGE: smallb4U_group[timelevels]\n",
    "STORAGE: smallb4_sq_group[timelevels]\n",
    "STORAGE: PoynSU_group[timelevels]\n",
    "\n",
    "schedule group sbPoynETNRPy_group in MoL_PseudoEvolution after ADMBase_SetADMVars\n",
    "{\n",
    "} \"Schedule sbPoynETNRPy group\"\n",
    "\n",
    "schedule sbPoynETNRPy in sbPoynETNRPy_group\n",
    "{\n",
    "    LANG: C\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",
    "\n",
    "    READS: admbase::alpha(Everywhere)\n",
    "\n",
    "    READS: admbase::betax(Everywhere)\n",
    "    READS: admbase::betay(Everywhere)\n",
    "    READS: admbase::betaz(Everywhere)\n",
    "\n",
    "    READS: HydroBase::vel(Everywhere)\n",
    "    READS: HydroBase::Bvec(Everywhere)\n",
    "\n",
    "    WRITES: sbPoynETNRPy::smallb4U0(Everywhere)\n",
    "    WRITES: sbPoynETNRPy::smallb4U1(Everywhere)\n",
    "    WRITES: sbPoynETNRPy::smallb4U2(Everywhere)\n",
    "    WRITES: sbPoynETNRPy::smallb4U3(Everywhere)\n",
    "\n",
    "    WRITES: sbPoynETNRPy::smallb4_sq(Everywhere)\n",
    "\n",
    "    WRITES: sbPoynETNRPy::PoynSU0(Everywhere)\n",
    "    WRITES: sbPoynETNRPy::PoynSU1(Everywhere)\n",
    "    WRITES: sbPoynETNRPy::PoynSU2(Everywhere)\n",
    "} \"Call sbPoynETNRPy main function, to compute $b^mu$, $b^2$, and $S^i$\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='etksys'></a> \n",
    "\n",
    "## Step 2.c: Inform the Einstein Toolkit build system of the C code \\[Back to [top](#toc)\\]\n",
    "$$\\label{etksys}$$\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": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:09:00.811344Z",
     "iopub.status.busy": "2021-03-07T17:09:00.810335Z",
     "iopub.status.idle": "2021-03-07T17:09:00.814477Z",
     "shell.execute_reply": "2021-03-07T17:09:00.815155Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting sbPoynETNRPy/src/make.code.defn\n"
     ]
    }
   ],
   "source": [
    "%%writefile sbPoynETNRPy/src/make.code.defn\n",
    "\n",
    "SRCS = sbPoynETNRPy.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-u0_smallb_Poynting.pdf](Tutorial-ETK_thorn-u0_smallb_Poynting.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": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:09:00.821342Z",
     "iopub.status.busy": "2021-03-07T17:09:00.820257Z",
     "iopub.status.idle": "2021-03-07T17:09:04.315406Z",
     "shell.execute_reply": "2021-03-07T17:09:04.314578Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ETK_thorn-u0_smallb_Poynting.tex, and compiled LaTeX file\n",
      "    to PDF file Tutorial-ETK_thorn-u0_smallb_Poynting.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-u0_smallb_Poynting\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}