{
 "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",
    "# Custom Functions for NRPy+\n",
    "\n",
    "## Author: Patrick Nelson\n",
    "\n",
    "## Introduction\n",
    "\n",
    "Using SymPy, NRPy+ is able to provide a lot of different functions out of the box, covering many different use cases. However, there are some more obscure functions that we can't access; we may also want to use a different method than the SymPy function uses or use data from multiple different points.\n",
    "\n",
    "The code below first tells SymPy that `nrpyMyFunc` should be treated as a SymPy function. This name is what should be used in the python code. Then, a dictionary is imported from `outputC` so we can add an entry to it. The key for the entry should be identical to the function; the entry (`myFunc_Ccode` in the example below) should be identical to the name of the function in the C code.\n",
    "```python\n",
    "import sympy as sp\n",
    "\n",
    "nrpyMyFunc = sp.Function('nrpyMyFunc')\n",
    "from outputC import custom_functions_for_SymPy_ccode\n",
    "custom_functions_for_SymPy_ccode[\"nrpyMyFunc\"] = \"myFunc_Ccode\"\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above method is not restricted to functions; macros can be used as well. Additionally, SymPy interprets the argument of the function quite generously; the argument of the SymPy function can be more SymPy code, a string, or even nothing (possibly more!). Consider the following examples; the Python code will be on the left and the resulting C code on the right:\n",
    "* `x = nrpyMyFunc(y)` $\\rightarrow$ `x = myFunc_Ccode(y)`\n",
    "    * Here, the Python symbol y can be any SymPy expression; the appropriate corresponding expression will be substituted in the resulting C code, and CSE will be applied as well.\n",
    "* `x = nrpyMyFunc()` $\\rightarrow$ `x = myFunc_Ccode()`\n",
    "* `x = nrpyMyFunc` $\\rightarrow$ `x = myFunc_Ccode`\n",
    "* `x = nrpyMyFunc(\"SOME_GF\")` $\\rightarrow$ `x = myFunc_Ccode(SOME_GF)`\n",
    "\n",
    "As you can see, we have considerable flexibility in how the C code turns out.\n",
    "\n",
    "### Obscure functions\n",
    "\n",
    "Most functions that we may need to use (such as trigonometric or logarithmic functions) are already included in SymPy. However, there are myriad other types of functions that can show up in various problems, such as the dilogarithm function, which is needed for the [Split Monopole initial data in `GiRaFFE_NRPy`](Tutorial-GiRaFFEfood_NRPy-Split_Monopole.ipynb). In this case, we were able to use a preexisting library, but this is not necessary in general; it is just as possible to add a function that you have written yourself. \n",
    "\n",
    "### Changing functionality\n",
    "\n",
    "Sometimes, SymPy may have the function you need, but implemented in a way that is not ideal. This was the case when we wrote [Tutorial-Min_Max_and_Piecewise_Expressions](../Tutorial-Min_Max_and_Piecewise_Expressions.ipynb). The SymPy implementation `sp.Abs` assumed a complex input (even when we specified that the input was real), resulting in a function that was much more computationally expensive than it needed to be and introducing errors into the produced C code. We solved this by adding `nrpyAbs` to the dictionary, pointing to the basic C implementation of the function.\n",
    "\n",
    "### Multiple gridpoints\n",
    "\n",
    "Operations like interpolation are critical to some applications (`GiRaFFE_NRPy` in particular depends on interpolation for its staggered grids). However, the only context in which data is read from multiple gridpoints in NRPy+ is finite-difference derivatives. For the interpolation of metric gridfunctions required in `GiRaFFE_NRPy`, we handwrote an entire function to replicate a macro in the original `GiRaFFE`, requiring several extra gridfunctions of storage compared to the original code. An alternative would be to port the macro to the `GiRaFFE_NRPy` C code and use the above method to allow SymPy-generated kernels to use it. \n",
    "\n",
    "A basic interpolator could be written as follows:\n",
    "```C\n",
    "#define A0 1.0\n",
    "#define A1 1.0\n",
    "#define METRIC   auxevol_gfs[IDX4S(GF_TO_INTERP, i0,i1,i2)]\n",
    "#define METRICp1 auxevol_gfs[IDX4S(GF_TO_INTERP, i0+(flux_dirn==0),i1+(flux_dirn==1),i2+(flux_dirn==2))]\n",
    "#define METRIC_INTERPED(GF_TO_INTERP) (A0*(METRIC) + A1*(METRICp1))\n",
    "```\n",
    "This computes the average of a values of the some gridfunction `GF_TO_INTERP` at the current point and the next one in a given `flux_dirn`, which is a basic interpolation of the value at the cell face. \n",
    "\n",
    "Naturally, one must take care to use the correct gridfunction macros as defined in the C code. A function could be written to automate this because SymPy expressions can be easily cast to strings and NRPy+'s gridfunction `#define`s are straightforwrard to generate. For instance, the following could be used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/*\n",
      " *  Original SymPy expression:\n",
      " *  \"x = nrpyMetricInterp(ALPHAGF)*nrpyMetricInterp(GAMMADD00GF)\"\n",
      " */\n",
      "{\n",
      "  x = METRIC_INTERPED(ALPHAGF)*METRIC_INTERPED(GAMMADD00GF);\n",
      "}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import os,sys\n",
    "nrpy_dir_path = os.path.join(\"..\")\n",
    "if nrpy_dir_path not in sys.path:\n",
    "    sys.path.append(nrpy_dir_path)\n",
    "\n",
    "from outputC import outputC,custom_functions_for_SymPy_ccode # NRPy+: Core C code output module\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\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",
    "\n",
    "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gammaDD\",\"sym01\")\n",
    "alpha = gri.register_gridfunctions(\"AUXEVOL\",\"alpha\")\n",
    "\n",
    "nrpyMetricInterp = sp.Function('nrpyMetricInterp')\n",
    "\n",
    "def gf_id_of(gridfunction):\n",
    "    return (str(gridfunction).upper() + \"GF\")\n",
    "\n",
    "custom_functions_for_SymPy_ccode[\"nrpyMetricInterp\"] = \"METRIC_INTERPED\"\n",
    "\n",
    "x = nrpyMetricInterp(gf_id_of(alpha)) * nrpyMetricInterp(gf_id_of(gammaDD[0][0]))\n",
    "outputC(x,\"x\")"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}