{
 "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",
    "# `GiRaFFE_NRPy`: Source Terms\n",
    "\n",
    "## Author: Patrick Nelson\n",
    "\n",
    "<a id='intro'></a>\n",
    "\n",
    "**Notebook Status:** <font color=green><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This code produces the expected results for generated functions.\n",
    "\n",
    "## This module presents the functionality of [GiRaFFE_NRPy_Source_Terms.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_Source_Terms.py).\n",
    "\n",
    "## Introduction: \n",
    "This writes and documents the C code that `GiRaFFE_NRPy` uses to compute the source terms for the right-hand sides of the evolution equations for the unstaggered prescription.\n",
    "\n",
    "The equations themselves are already coded up in other functions; however, for the $\\tilde{S}_i$ source term, we will need derivatives of the metric. It will be most efficient and accurate to take them using the interpolated metric values that we will have calculated anyway; however, we will need to write our derivatives in a nonstandard way within NRPy+ in order to take advantage of this, writing our own code for memory access."
   ]
  },
  {
   "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](#stilde_source): The $\\tilde{S}_i$ source term\n",
    "1. [Step 2](#code_validation): Code Validation against original C code\n",
    "1. [Step 3](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 0: Add NRPy's directory to the path\n",
    "# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n",
    "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",
    "import cmdline_helper as cmd\n",
    "outdir = os.path.join(\"GiRaFFE_NRPy\",\"GiRaFFE_Ccode_validation\",\"RHSs\")\n",
    "cmd.mkdir(outdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='stilde_source'></a>\n",
    "\n",
    "## Step 1: The $\\tilde{S}_i$ source term \\[Back to [top](#toc)\\]\n",
    "$$\\label{stilde_source}$$\n",
    "\n",
    "We start in the usual way - import the modules we need. We will also import the Levi-Civita symbol from `indexedexp.py` and use it to set the Levi-Civita tensor $\\epsilon^{ijk} = [ijk]/\\sqrt{\\gamma}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 1: The StildeD RHS *source* term\n",
    "from outputC import outputC, outCfunction, add_to_Cfunction_dict # NRPy+: Core C code output module\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import GRHD.equations as GRHD    # NRPy+: Generate general relativistic hydrodynamics equations\n",
    "import GRFFE.equations as GRFFE  # NRPy+: Generate general relativistic force-free electrodynamics equations\n",
    "\n",
    "thismodule = \"GiRaFFE_NRPy_Source_Terms\"\n",
    "\n",
    "def generate_memory_access_code():\n",
    "    # There are several pieces of C code that we will write ourselves because we need to do things\n",
    "    # a little bit outside of what NRPy+ is built for.\n",
    "    # First, we will write general memory access. We will read in values from memory at a given point\n",
    "    # for each quantity we care about.\n",
    "    global general_access\n",
    "    general_access = \"\"\n",
    "    for var in [\"GAMMADD00\", \"GAMMADD01\", \"GAMMADD02\",\n",
    "                \"GAMMADD11\", \"GAMMADD12\", \"GAMMADD22\",\n",
    "                \"BETAU0\", \"BETAU1\", \"BETAU2\",\"ALPHA\",\n",
    "                \"BU0\",\"BU1\",\"BU2\",\n",
    "                \"VALENCIAVU0\",\"VALENCIAVU1\",\"VALENCIAVU2\"]:\n",
    "        lhsvar = var.lower().replace(\"dd\",\"DD\").replace(\"u\",\"U\").replace(\"bU\",\"BU\").replace(\"valencia\",\"Valencia\")\n",
    "        # e.g.,\n",
    "        # const REAL gammaDD00dD0 = auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1,i2)];\n",
    "        general_access += \"const REAL \"+lhsvar+\" = auxevol_gfs[IDX4S(\"+var+\"GF,i0,i1,i2)];\\n\"\n",
    "\n",
    "    # This quick function returns a nearby point for memory access. We need this because derivatives are not local operations.\n",
    "    def idxp1(dirn):\n",
    "        if dirn==0:\n",
    "            return \"i0+1,i1,i2\"\n",
    "        if dirn==1:\n",
    "            return \"i0,i1+1,i2\"\n",
    "        if dirn==2:\n",
    "            return \"i0,i1,i2+1\"\n",
    "\n",
    "    # Next we evaluate needed derivatives of the metric, based on their values at cell faces\n",
    "    global metric_deriv_access\n",
    "    metric_deriv_access = []\n",
    "    for dirn in range(3):\n",
    "        metric_deriv_access.append(\"\")\n",
    "        for var in [\"GAMMA_FACEDDdD00\", \"GAMMA_FACEDDdD01\", \"GAMMA_FACEDDdD02\",\n",
    "                    \"GAMMA_FACEDDdD11\", \"GAMMA_FACEDDdD12\", \"GAMMA_FACEDDdD22\",\n",
    "                    \"BETA_FACEUdD0\", \"BETA_FACEUdD1\", \"BETA_FACEUdD2\",\"ALPHA_FACEdD\"]:\n",
    "            lhsvar = var.lower().replace(\"dddd\",\"DDdD\").replace(\"udd\",\"UdD\").replace(\"dd\",\"dD\").replace(\"u\",\"U\").replace(\"_face\",\"\")\n",
    "            rhsvar = var.replace(\"dD\",\"\")\n",
    "            # e.g.,\n",
    "            # const REAL gammaDDdD000 = (auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1,i2)])/dxx0;\n",
    "            metric_deriv_access[dirn] += \"const REAL \"+lhsvar+str(dirn)+\" = (auxevol_gfs[IDX4S(\"+rhsvar+\"GF,\"+idxp1(dirn)+\")]-auxevol_gfs[IDX4S(\"+rhsvar+\"GF,i0,i1,i2)])/dxx\"+str(dirn)+\";\\n\"\n",
    "        metric_deriv_access[dirn] += \"REAL Stilde_rhsD\"+str(dirn)+\";\\n\"\n",
    "\n",
    "    # This creates the C code that writes to the Stilde_rhs direction specified.\n",
    "    global write_final_quantity\n",
    "    write_final_quantity = []\n",
    "    for dirn in range(3):\n",
    "        write_final_quantity.append(\"\")\n",
    "        write_final_quantity[dirn] += \"rhs_gfs[IDX4S(STILDED\"+str(dirn)+\"GF,i0,i1,i2)] += Stilde_rhsD\"+str(dirn)+\";\"\n",
    "\n",
    "def write_out_functions_for_StildeD_source_term(outdir,outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi):\n",
    "    generate_memory_access_code()\n",
    "    # First, we declare some dummy tensors that we will use for the codegen.\n",
    "    gammaDDdD  = ixp.declarerank3(\"gammaDDdD\",\"sym01\",DIM=3)\n",
    "    betaUdD = ixp.declarerank2(\"betaUdD\",\"nosym\",DIM=3)\n",
    "    alphadD = ixp.declarerank1(\"alphadD\",DIM=3)\n",
    "\n",
    "    # We need to rerun a few of these functions with the reset lists to make sure these functions\n",
    "    # don't cheat by using analytic expressions\n",
    "    GRHD.compute_sqrtgammaDET(gammaDD)\n",
    "    GRHD.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha, betaU, gammaDD, ValenciavU)\n",
    "    GRFFE.compute_smallb4U(gammaDD, betaU, alpha, GRHD.u4U_ito_ValenciavU, BU, sqrt4pi)\n",
    "    GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)\n",
    "    GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4U, GRFFE.smallbsquared,GRHD.u4U_ito_ValenciavU)\n",
    "    GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDDdD,betaUdD,alphadD)\n",
    "    GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD, GRFFE.TEM4UU)\n",
    "    for i in range(3):\n",
    "        desc = \"Adds the source term to StildeD\"+str(i)+\".\"\n",
    "        name = \"calculate_StildeD\"+str(i)+\"_source_term\"\n",
    "        outCfunction(\n",
    "            outfile  = os.path.join(outdir,name+\".h\"), desc=desc, name=name,\n",
    "            params   =\"const paramstruct *params,const REAL *auxevol_gfs, REAL *rhs_gfs\",\n",
    "            body     = general_access \\\n",
    "                      +metric_deriv_access[i]\\\n",
    "                      +outputC(GRHD.S_tilde_source_termD[i],\"Stilde_rhsD\"+str(i),\"returnstring\",params=outCparams)\\\n",
    "                      +write_final_quantity[i],\n",
    "            loopopts =\"InteriorPoints\",\n",
    "            rel_path_to_Cparams=os.path.join(\"../\"))\n",
    "\n",
    "def add_to_Cfunction_dict__functions_for_StildeD_source_term(outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi,\n",
    "                                                             includes=None, rel_path_to_Cparams=os.path.join(\"../\"),\n",
    "                                                             path_from_rootsrcdir_to_this_Cfunc=os.path.join(\"RHSs/\")):\n",
    "    generate_memory_access_code()\n",
    "    # First, we declare some dummy tensors that we will use for the codegen.\n",
    "    gammaDDdD  = ixp.declarerank3(\"gammaDDdD\",\"sym01\",DIM=3)\n",
    "    betaUdD = ixp.declarerank2(\"betaUdD\",\"nosym\",DIM=3)\n",
    "    alphadD = ixp.declarerank1(\"alphadD\",DIM=3)\n",
    "\n",
    "    # We need to rerun a few of these functions with the reset lists to make sure these functions\n",
    "    # don't cheat by using analytic expressions\n",
    "    GRHD.compute_sqrtgammaDET(gammaDD)\n",
    "    GRHD.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha, betaU, gammaDD, ValenciavU)\n",
    "    GRFFE.compute_smallb4U(gammaDD, betaU, alpha, GRHD.u4U_ito_ValenciavU, BU, sqrt4pi)\n",
    "    GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)\n",
    "    GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4U, GRFFE.smallbsquared,GRHD.u4U_ito_ValenciavU)\n",
    "    GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDDdD,betaUdD,alphadD)\n",
    "    GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD, GRFFE.TEM4UU)\n",
    "    for i in range(3):\n",
    "        desc = \"Adds the source term to StildeD\"+str(i)+\".\"\n",
    "        name = \"calculate_StildeD\"+str(i)+\"_source_term\"\n",
    "        params   =\"const paramstruct *params,const REAL *auxevol_gfs, REAL *rhs_gfs\"\n",
    "        body     = general_access \\\n",
    "                  +metric_deriv_access[i]\\\n",
    "                  +outputC(GRHD.S_tilde_source_termD[i],\"Stilde_rhsD\"+str(i),\"returnstring\",params=outCparams)\\\n",
    "                  +write_final_quantity[i]\n",
    "        loopopts =\"InteriorPoints\"\n",
    "        add_to_Cfunction_dict(\n",
    "            includes=includes,\n",
    "            desc=desc,\n",
    "            name=name, params=params,\n",
    "            body=body, loopopts=loopopts,\n",
    "            path_from_rootsrcdir_to_this_Cfunc=path_from_rootsrcdir_to_this_Cfunc,\n",
    "            rel_path_to_Cparams=rel_path_to_Cparams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 2: Code Validation against original C code \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "To validate the code in this tutorial we check for agreement between the files\n",
    "\n",
    "1. that were written in this tutorial and\n",
    "1. those that are stored in `GiRaFFE_NRPy/GiRaFFE_Ccode_library` or generated by `GiRaFFE_NRPy_A2B.py`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function calculate_StildeD0_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_validation\\RHSs\\calculate_StildeD0_source_term.h\n",
      "Output C function calculate_StildeD1_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_validation\\RHSs\\calculate_StildeD1_source_term.h\n",
      "Output C function calculate_StildeD2_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_validation\\RHSs\\calculate_StildeD2_source_term.h\n",
      "Output C function calculate_StildeD0_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_library\\RHSs\\calculate_StildeD0_source_term.h\n",
      "Output C function calculate_StildeD1_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_library\\RHSs\\calculate_StildeD1_source_term.h\n",
      "Output C function calculate_StildeD2_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_library\\RHSs\\calculate_StildeD2_source_term.h\n",
      "Printing difference between original C code and this code...\n",
      "Checking file calculate_StildeD0_source_term.h\n",
      "No difference. TEST PASSED!\n",
      "Checking file calculate_StildeD1_source_term.h\n",
      "No difference. TEST PASSED!\n",
      "Checking file calculate_StildeD2_source_term.h\n",
      "No difference. TEST PASSED!\n"
     ]
    }
   ],
   "source": [
    "# Declare gridfunctions necessary to generate the C code:\n",
    "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gammaDD\",\"sym01\",DIM=3)\n",
    "betaU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"betaU\",DIM=3)\n",
    "alpha = gri.register_gridfunctions(\"AUXEVOL\",\"alpha\",DIM=3)\n",
    "BU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"BU\",DIM=3)\n",
    "ValenciavU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"ValenciavU\",DIM=3)\n",
    "StildeD = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"StildeD\",DIM=3)\n",
    "# Declare this symbol:\n",
    "sqrt4pi = par.Cparameters(\"REAL\",thismodule,\"sqrt4pi\",\"sqrt(4.0*M_PI)\")\n",
    "\n",
    "# First, we generate the file using the functions written in this notebook:\n",
    "outCparams = \"outCverbose=False\"\n",
    "write_out_functions_for_StildeD_source_term(outdir,outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi)\n",
    "\n",
    "# Define the directory that we wish to validate against:\n",
    "valdir = os.path.join(\"GiRaFFE_NRPy\",\"GiRaFFE_Ccode_library\",\"RHSs\")\n",
    "cmd.mkdir(valdir)\n",
    "\n",
    "import GiRaFFE_NRPy.GiRaFFE_NRPy_Source_Terms as source\n",
    "source.write_out_functions_for_StildeD_source_term(valdir,outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi)\n",
    "\n",
    "import difflib\n",
    "import sys\n",
    "\n",
    "print(\"Printing difference between original C code and this code...\")\n",
    "# Open the files to compare\n",
    "files = [\"calculate_StildeD0_source_term.h\",\"calculate_StildeD1_source_term.h\",\"calculate_StildeD2_source_term.h\"]\n",
    "\n",
    "for file in files:\n",
    "    print(\"Checking file \" + file)\n",
    "    with open(os.path.join(valdir,file)) as file1, open(os.path.join(outdir,file)) as file2:\n",
    "        # Read the lines of each file\n",
    "        file1_lines = file1.readlines()\n",
    "        file2_lines = file2.readlines()\n",
    "        num_diffs = 0\n",
    "        for line in difflib.unified_diff(file1_lines, file2_lines, fromfile=os.path.join(valdir+file), tofile=os.path.join(outdir+file)):\n",
    "            sys.stdout.writelines(line)\n",
    "            num_diffs = num_diffs + 1\n",
    "        if num_diffs == 0:\n",
    "            print(\"No difference. TEST PASSED!\")\n",
    "        else:\n",
    "            print(\"ERROR: Disagreement found with .py file. See differences above.\")\n",
    "            sys.exit(1)"
   ]
  },
  {
   "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-GiRaFFE_NRPy_C_code_library-Source_Terms](TTutorial-GiRaFFE_NRPy_C_code_library-Source_Terms.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": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Notebook output to PDF is only supported on Linux systems, with pdflatex installed.\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-GiRaFFE_NRPy-Source_Terms\",location_of_template_file=os.path.join(\"..\"))"
   ]
  }
 ],
 "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": 2
}