Tutorial-ETK_thorn-IDScalarWaveNRPy.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",
    "# IDScalarWaveNRPy: An Einstein Toolkit Initial Data Thorn for the Scalar Wave Equation\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "**Notebook Status and Validation Notes: TODO**\n",
    "\n",
    "### NRPy+ Source Code for this module: [ScalarWave/InitialData.py](../edit/ScalarWave/InitialData.py) [\\[**tutorial**\\]](Tutorial-ScalarWave.ipynb) Contructs the SymPy expressions for plane-wave initial data\n",
    "\n",
    "## Introduction:\n",
    "In this part of the tutorial, we will construct an Einstein Toolkit (ETK) thorn (module) that will set up *initial data* for the scalar wave initial value problem. In a [previous tutorial notebook](Tutorial-ScalarWave.ipynb), we used NRPy+ to contruct the SymPy expressions for plane-wave initial data. This thorn is largely based on and should function similarly to the $\\text{IDScalarWaveC}$ thorn included in the Einstein Toolkit (ETK) $\\text{CactusWave}$ arrangement.\n",
    "\n",
    "We will construct this thorn in two steps.\n",
    "\n",
    "1. Call on NRPy+ to convert the SymPy expressions for the initial data into one C-code kernel.\n",
    "1. Write the C code and linkages to the Einstein Toolkit infrastructure (i.e., the .ccl files) to complete this Einstein Toolkit module."
   ]
  },
  {
   "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): Call on NRPy+ to convert the SymPy expression for the scalar wave initial data into a C-code kernel\n",
    "1. [Step 2](#einstein): Interfacing with the Einstein Toolkit\n",
    "    1. [Step 2.a](#einstein_c): Constructing the Einstein Toolkit C-code calling functions that include the C code kernels\n",
    "    1. [Step 2.b](#einstein_ccl): CCL files - Define how this module interacts and interfaces with the larger Einstein Toolkit infrastructure\n",
    "    1. [Step 2.c](#einstein_list): Add the C code to the 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='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Call on NRPy+ to convert the SymPy expression for the scalar wave initial data into a C-code kernel \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "After importing the core modules, since we are writing an ETK thorn, we'll need to set `\"grid::GridFuncMemAccess\"` to `\"ETK\"`. SymPy expressions for plane wave initial data are written inside [ScalarWave/InitialData.py](../edit/ScalarWave/InitialData.py), and we simply import them for use here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: Call on NRPy+ to convert the SymPy expression\n",
    "#         for the scalar wave initial data into a C-code kernel\n",
    "\n",
    "# Step 1a: Import needed NRPy+ core modules:\n",
    "import NRPy_param_funcs as par  # NRPy+: Parameter interface\n",
    "import indexedexp as ixp        # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import grid as gri              # NRPy+: Functionality for handling numerical grids\n",
    "import finite_difference as fin # NRPy+: Finite difference C code generation module\n",
    "from outputC import lhrh        # NRPy+: Core C code output module\n",
    "import loop                     # NRPy+: Generate C code loops\n",
    "\n",
    "# Step 1b: This is an Einstein Toolkit (ETK) thorn. Here we\n",
    "#          tell NRPy+ that gridfunction memory access will\n",
    "#          therefore be in the \"ETK\" style.\n",
    "par.set_parval_from_str(\"grid::GridFuncMemAccess\",\"ETK\")\n",
    "\n",
    "# Step 1c: Call the InitialData() function from within the\n",
    "#          ScalarWave/InitialData.py module.\n",
    "import ScalarWave.InitialData as swid\n",
    "\n",
    "# Step 1d: Within the ETK, the 3D gridfunctions x, y, and z store the\n",
    "#         Cartesian grid coordinates. Setting the gri.xx[] arrays\n",
    "#         to point to these gridfunctions forces NRPy+ to treat\n",
    "#         the Cartesian coordinate gridfunctions properly --\n",
    "#         reading them from memory as needed.\n",
    "x,y,z = gri.register_gridfunctions(\"AUX\",[\"x\",\"y\",\"z\"])\n",
    "gri.xx[0] = x\n",
    "gri.xx[1] = y\n",
    "gri.xx[2] = z\n",
    "\n",
    "# Step 1e: Set up the plane wave initial data. This sets uu_ID and vv_ID.\n",
    "swid.InitialData(CoordSystem=\"Cartesian\",Type=\"PlaneWave\")\n",
    "\n",
    "# Step 1f: Register uu and vv gridfunctions so they can be written to by NRPy.\n",
    "uu,vv = gri.register_gridfunctions(\"EVOL\",[\"uu\",\"vv\"])\n",
    "\n",
    "# Step 1g: Set the uu and vv gridfunctions to the uu_ID & vv_ID variables\n",
    "#         defined by InitialData().\n",
    "uu = swid.uu_ID\n",
    "vv = swid.vv_ID\n",
    "\n",
    "# Step 1h: Create the C code output kernel.\n",
    "scalar_PWID_to_print = [\\\n",
    "                        lhrh(lhs=gri.gfaccess(\"out_gfs\",\"uu\"),rhs=uu),\\\n",
    "                        lhrh(lhs=gri.gfaccess(\"out_gfs\",\"vv\"),rhs=vv),]\n",
    "scalar_PWID_CcodeKernel = fin.FD_outputC(\"returnstring\",scalar_PWID_to_print)\n",
    "\n",
    "scalar_PWID_looped = loop.loop([\"i2\",\"i1\",\"i0\"],[\"0\",\"0\",\"0\"],[\"cctk_lsh[2]\",\"cctk_lsh[1]\",\"cctk_lsh[0]\"],\\\n",
    "                               [\"1\",\"1\",\"1\"],[\"#pragma omp parallel for\",\"\",\"\"],\"\",\\\n",
    "                               scalar_PWID_CcodeKernel.replace(\"time\",\"cctk_time\"))\n",
    "\n",
    "# Step 1i: Create directories for the thorn if they don't exist.\n",
    "!mkdir IDScalarWaveNRPy     2>/dev/null # 2>/dev/null: Don't throw an error if the directory already exists.\n",
    "!mkdir IDScalarWaveNRPy/src 2>/dev/null # 2>/dev/null: Don't throw an error if the directory already exists.\n",
    "\n",
    "# Step 1j: Write the C code kernel to file.\n",
    "with open(\"IDScalarWaveNRPy/src/ScalarWave_PWID.h\", \"w\") as file:\n",
    "    file.write(str(scalar_PWID_looped))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein'></a>\n",
    "\n",
    "# Step 2: Interfacing with the Einstein Toolkit \\[Back to [top](#toc)\\]\n",
    "$$\\label{einstein}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein_c'></a>\n",
    "\n",
    "## Step 2.a: Constructing the Einstein Toolkit C-code calling functions that include the C code kernels \\[Back to [top](#toc)\\]\n",
    "$$\\label{einstein_c}$$\n",
    "\n",
    "We will write another C file with the functions we need here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing IDScalarWaveNRPy/src/InitialData.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile IDScalarWaveNRPy/src/InitialData.c\n",
    "#include <math.h>\n",
    "#include <stdio.h>\n",
    "\n",
    "#include \"cctk.h\"\n",
    "#include \"cctk_Parameters.h\"\n",
    "#include \"cctk_Arguments.h\"\n",
    "\n",
    "void IDScalarWaveNRPy_param_check(CCTK_ARGUMENTS) {\n",
    "  DECLARE_CCTK_ARGUMENTS;\n",
    "  DECLARE_CCTK_PARAMETERS;\n",
    "\n",
    "  if (kk0 == 0 && kk1 == 0 && kk2 == 0) {\n",
    "     CCTK_WARN(0,\"kk0==kk1==kk2==0: Zero wave vector cannot be normalized. Set one of the kk's to be != 0.\");\n",
    "  }\n",
    "}\n",
    "\n",
    "void IDScalarWaveNRPy_InitialData(CCTK_ARGUMENTS)\n",
    "{\n",
    "  DECLARE_CCTK_ARGUMENTS\n",
    "  DECLARE_CCTK_PARAMETERS\n",
    "\n",
    "  const CCTK_REAL *xGF = x;\n",
    "  const CCTK_REAL *yGF = y;\n",
    "  const CCTK_REAL *zGF = z;\n",
    "#include \"ScalarWave_PWID.h\"\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein_ccl'></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{einstein_ccl}$$\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. Specifically, this file governs the interaction between this thorn and others; more information can be found in the [official Einstein Toolkit documentation](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-178000D2.2). \n",
    "With \"implements\", we give our thorn its unique name. By \"inheriting\" other thorns, we tell the Toolkit that we will rely on variables that exist and are declared \"public\" within those functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing IDScalarWaveNRPy/interface.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile IDScalarWaveNRPy/interface.ccl\n",
    "\n",
    "implements: IDScalarWaveNRPy\n",
    "\n",
    "inherits: WaveToyNRPy grid\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. `param.ccl`: specifies free parameters within the thorn, enabling them to be set at runtime. It is required to provide 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)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing IDScalarWaveNRPy/param.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile IDScalarWaveNRPy/param.ccl\n",
    "shares: grid\n",
    "\n",
    "USES KEYWORD type\n",
    "\n",
    "shares: WaveToyNRPy\n",
    "\n",
    "USES REAL wavespeed\n",
    "\n",
    "restricted:\n",
    "CCTK_KEYWORD initial_data \"Type of initial data\"\n",
    "{\n",
    "  \"plane\"      :: \"Plane wave\"\n",
    "} \"plane\"\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL kk0 \"The wave number in the x-direction\"\n",
    "{\n",
    " *:* :: \"No restriction\"\n",
    "} 4.0\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL kk1 \"The wave number in the y-direction\"\n",
    "{\n",
    " *:* :: \"No restriction\"\n",
    "} 0.0\n",
    "\n",
    "restricted:\n",
    "CCTK_REAL kk2 \"The wave number in the z-direction\"\n",
    "{\n",
    " *:* :: \"No restriction\"\n",
    "} 0.0\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. `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. `schedule.ccl`'s official documentation may be found [here](http://einsteintoolkit.org/usersguide/UsersGuidech12.html#x17-186000D2.4). \n",
    "\n",
    "We specify here the standardized ETK \"scheduling bins\" in which we want each of our thorn's functions to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing IDScalarWaveNRPy/schedule.ccl\n"
     ]
    }
   ],
   "source": [
    "%%writefile IDScalarWaveNRPy/schedule.ccl\n",
    "\n",
    "schedule IDScalarWaveNRPy_param_check at CCTK_PARAMCHECK\n",
    "{\n",
    "  LANG: C\n",
    "  OPTIONS: global\n",
    "} \"Check sanity of parameters\"\n",
    "\n",
    "schedule IDScalarWaveNRPy_InitialData at CCTK_INITIAL as WaveToy_InitialData\n",
    "{\n",
    "  STORAGE:       WaveToyNRPy::scalar_fields[3]\n",
    "  LANG:          C\n",
    "  READS: grid::x(Everywhere)\n",
    "  READS: grid::y(Everywhere)\n",
    "  READS: grid::y(Everywhere)\n",
    "  WRITES: uuGF(Everywhere)\n",
    "  WRITES: vvGF(Everywhere)\n",
    "} \"Initial data for 3D wave equation\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='einstein_list'></a>\n",
    "\n",
    "## Step 2.c: Add the C code to the Einstein Toolkit compilation list \\[Back to [top](#toc)\\]\n",
    "$$\\label{einstein_list}$$\n",
    "\n",
    "We will also need `make.code.defn`, which indicates the list of files that need to be compiled. This thorn only has the one C file to compile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing IDScalarWaveNRPy/src/make.code.defn\n"
     ]
    }
   ],
   "source": [
    "%%writefile IDScalarWaveNRPy/src/make.code.defn\n",
    "\n",
    "SRCS = InitialData.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-IDScalarWaveNRPy.pdf](Tutorial-ETK_thorn-IDScalarWaveNRPy.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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ETK_thorn-IDScalarWaveNRPy.tex, and compiled LaTeX file to\n",
      "    PDF file Tutorial-ETK_thorn-IDScalarWaveNRPy.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-IDScalarWaveNRPy\")"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}