Tutorial-Numerical_Grids.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",
    "# Numerical Grids in NRPy+\n",
    "\n",
    "## Author: Zach Etienne\n",
    "\n",
    "### NRPy+ Source Code for this module: [grid.py](../edit/grid.py)\n",
    "\n",
    "Solving partial differential equations on the computer with finite difference techniques requires that we sample our solutions to these equations on numerical grids. We call all sampled functions that are stored on our numerical grid *gridfunctions*. NRPy+'s grid module adds the capability of registering gridfunctions in NRPy, setting basic parameters of a numerical grid, and providing functions to other modules regarding reading and writing of gridfunctions to memory in the C code.\n",
    "\n",
    "Parameters in this module include:\n",
    "* **grid::DIM** -- the dimension of the grid (e.g., a 3D numerical grid will have grid::DIM=3)\n",
    "* **grid::Nx\\[DIM\\]** -- an integer array yielding the size of the grid in each direction. E.g., in Cartesian coordinates Nx\\[0\\] will be set to the number of grid points in the x direction. *This is a parameter that is set at C runtime, not in NRPy+; NRPy+ simply generates the declaration of this parameter in the C code.*\n",
    "* **grid::MemAllocStyle** -- how the gridfunction is allocated in memory in the C code. This is used when generating the C code to ensure that gridfunctions are read as sequentially in memory as possible, to avoid [cache misses](https://en.wikipedia.org/wiki/CPU_cache#CACHE-MISS). There are currently two MemAllocStyles supported:\n",
    "    * If the following loop accesses the grid function in a contiguous manner in memory, we set \"grid::MemAllocStyle=012\":\n",
    "        * for(int i0=0;i0<Nx\\[0\\];i0++) for(int i1=0;i1<Nx\\[1\\];i1++) for(int i2=0;i2<Nx\\[2\\];i2++) {\n",
    "    * Alternatively, the \"grid::MemAllocStyle=210\" is the reverse:\n",
    "        * for(int i2=0;i2<Nx\\[2\\];i2++) for(int i1=0;i1<Nx\\[1\\];i1++) for(int i0=0;i0<Nx\\[0\\];i0++) {\n",
    "* **grid::GridFuncMemAccess** -- specifies how gridfunction data is accessed from memory. For example,\n",
    "    * In the Einstein Toolkit (grid::GridFuncMemAccess=\"ETK\"), the datum from gridfunction dummy at point (i0,j0,k0) is accessed via \"dummyGF\\[CCTK_GFINDEX3D(cctkGH,i0,j0,k0)\\]\".\n",
    "    * In SENR (grid::GridFuncMemAccess=\"SENRlike\"), the datum from gridfunction dummy in gridfunction array gfarry at point (i0,j0,k0) is accessed via \"gfarray\\[IDX4D(DUMMYGF,i0,j0,k0)\\]\".\n",
    "    * *Special note*: NRPy+ is code agnostic; additional types can be easily be added by modifying the function gfaccess() in the NRPy+ grid module. It should only require about 5 lines of code.\n",
    "\n",
    "Functions in this module include:\n",
    "* **gfaccess(gfarrayname = \"\",varname = \"\",ijklstring = \"\")**: Given a gridfunction array name, a variable name, and a string indicating the coordinates of the point in memory, return the string to access the gridfunction in memory at this data point. See grid::GridFuncMemAccess parameter above for a complete description.\n",
    "* **register_gridfunctions(gf_type,gf_names)**: returns either a single gridfunction of list of gridfunctions as SymPy variables\n",
    "    * gf_type can be set to either \"EVOL\" or \"AUX\".\n",
    "        * Setting to \"EVOL\" denotes a grid function that will need to be stepped forward in time within the C code's timestepping algorithm.\n",
    "        * Setting to \"AUX\" denotes any other grid function.\n",
    "    * gf_names can be a single gridfunction or a set of gridfunctions, all with the same gf_type.\n",
    "* **variable_type(var)**: First searches the list of registered gridfunctions and parameters for gridfunctions or C parameters; then outputs \"gridfunction\" or \"Cparameter\" if one or the other is found, respectively. Otherwise it will output \"other\"\n",
    "* **output__gridfunction_defines_h__return_gf_lists(outdir)**: Given output directory \"outdir\" as input, output a C code file called \"outdir/gridfunction_defines.h\", which will `#define` all the gridfunction aliases. \n",
    "    + The function returns three lists, corresponding to the names (strings) of the evolved, auxiliary, and \"auxevol\" gridfunction names respectively. Each is interpreted by NRPy+ as follows:\n",
    "        + \"evolved\" gridfunctions in NRPy+ will be automatically handled with the [MoLtimestepping](Tutorial-Method_of_Lines-C_Code_Generation.ipynb) module, including memory allocation (`RK_Allocate_Memory.h`), MoL updates (`RK_MoL.h`), and memory deallocation (`RK_Free_Memory.h`).\n",
    "        + \"auxevol\" gridfunctions in NRPy+ must be manually allocated and deallocated by the user. These gridfunctions provide additional storage for gridfunctions needed within the MoL step.\n",
    "        + \"auxiliary\" gridfunctions in NRPy+, soon to be ***deprecated*** and replaced by the \"diagnostic\" gridfunction type, are used for diagnostic purposes only, at the end of each timestep. As only a single MoL timelevel is needed to launch the next timestep, \"auxiliary\" gridfunctions use memory allocated for \"evolved\" gridfunctions at a different timelevel and can be reached via the `*diagnostic_gfs` pointer. \"auxiliary\"/\"diagnostic\" gridfunctions therefore cannot be used with \"Euler\" timestepping, and one cannot use more \"auxiliary\"/\"diagnostic\" gridfunctions than the number of \"evolved\" gridfunctions.\n",
    "    + The aliases `#define`'d in \"outdir/gridfunction_defines.h\" are meant to be human-friendly, so that accessing each gridfunction in C code can be done by its human-friendly alias (e.g., test_gfs\\[IDX4(VVGF,i0,i1,i2\\] instead of the less-friendly test_gfs\\[IDX4(1,i0,i1,i2\\]). \n",
    "    * Example: if we register with NRPy+ only two gridfunctions uu and vv, which are evolved quantities (i.e., represent the solution of the PDEs we are solving, and are registered with gftype == \"EVOL\"), then the first returned list (all gridfunctions  registered as EVOL) will be \\[\"uu\",\"vv\"\\], and the second (all gridfunctions registered as AUX) will be the empty list: \\[\\]. Also, this function will create a file with the following content:\n",
    "\n",
    "```C\n",
    "/* This file is automatically generated by NRPy+. Do not edit. */\n",
    "/* EVOLVED VARIABLES: */\n",
    "#define NUM_EVOL_GFS 2\n",
    "#define UUGF 0\n",
    "#define VVGF 1\n",
    "\n",
    "/* AUXILIARY VARIABLES: */\n",
    "#define NUM_AUX_GFS 0\n",
    "\n",
    "/* AUXEVOL VARIABLES: */\n",
    "#define NUM_AUXEVOL_GFS 0\n",
    "```\n",
    "\n",
    "Let's now register a gridfunction called \"phi\" in the \"in_gfs\" gridfunction array, then print C code that specifies how to access the gridfunction from the single point in memory (i0,i1,i2). Next we will demonstrate that phi is a normal SymPy variable, and finally we will confirm that its type is \"gridfunction\" using variable_type():"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Here's how to access the gridfunction phi, in grid array in_gfs, at point (i0,i1,i2):\n",
      "SENR-like memory access: aux_gfs[IDX4(PHIGF, i0,i1,i2)]\n",
      "ETK memory access: phiGF[CCTK_GFINDEX3D(cctkGH, i0,i1,i2)]\n",
      "\n",
      "\n",
      "To demonstrate that phi is just a regular SymPy variable, we will now print its square:\n",
      "phi**2\n",
      "\n",
      "\n",
      "\"phi\" is of type \"gridfunction\"\n",
      "\n",
      "\n",
      "\n",
      "Here is the list of registered evolved variables:  []\n",
      "... and here is the list of registered auxiliary variables:  ['phi']\n",
      "... and here is the list of registered auxevol variables:  []\n",
      "\n",
      "\n",
      " ... and here is the output gridfunction_defines.h file:\n",
      "/* This file is automatically generated by NRPy+. Do not edit. */\n",
      "\n",
      "/* EVOLVED VARIABLES: */\n",
      "#define NUM_EVOL_GFS 0\n",
      "\n",
      "\n",
      " /* AUXILIARY VARIABLES: */\n",
      "#define NUM_AUX_GFS 1\n",
      "#define PHIGF\t0\n",
      "\n",
      "\n",
      " /* AUXEVOL VARIABLES: */\n",
      "#define NUM_AUXEVOL_GFS 0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import grid as gri\n",
    "import NRPy_param_funcs as par\n",
    "import cmdline_helper as cmd\n",
    "\n",
    "# Register gridfunction phi, as type \"AUX\".\n",
    "# WARNING: register_gridfunctions can only be run once on a given gridfunction;\n",
    "#          you'll need to reset the Jupyter kernel before trying again:\n",
    "phi = gri.register_gridfunctions(\"AUX\",\"phi\")\n",
    "\n",
    "print(\"Here's how to access the gridfunction phi, in grid array in_gfs, at point (i0,i1,i2):\")\n",
    "print(\"SENR-like memory access: \"+gri.gfaccess(\"in_gfs\",\"phi\",\"i0,i1,i2\"))\n",
    "\n",
    "par.set_paramsvals_value(\"grid::GridFuncMemAccess = ETK\")\n",
    "print(\"ETK memory access: \"+gri.gfaccess(\"in_gfs\",\"phi\",\"i0,i1,i2\"))\n",
    "\n",
    "# Note that phi can now be used as a usual SymPy variable:\n",
    "print(\"\\n\\nTo demonstrate that phi is just a regular SymPy variable, we will now print its square:\")\n",
    "print(phi**2)\n",
    "\n",
    "print(\"\\n\\n\\\"phi\\\" is of type \\\"\"+ gri.variable_type(phi)+\"\\\"\")\n",
    "\n",
    "cmd.mkdir(\"NumGrids-tmp\")\n",
    "evolved_variables_list,auxiliary_variables_list, auxevol_variables_list = \\\n",
    "        gri.output__gridfunction_defines_h__return_gf_lists(\"NumGrids-tmp\")\n",
    "print(\"\\n\\n\")\n",
    "print(\"Here is the list of registered evolved variables: \",evolved_variables_list)\n",
    "print(\"... and here is the list of registered auxiliary variables: \",auxiliary_variables_list)\n",
    "print(\"... and here is the list of registered auxevol variables: \",auxevol_variables_list)\n",
    "print(\"\\n\\n ... and here is the output gridfunction_defines.h file:\")\n",
    "f = open(\"NumGrids-tmp/gridfunction_defines.h\",\"r\")\n",
    "print(f.read())\n",
    "f.close()"
   ]
  }
 ],
 "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
}