Tutorial-Method_of_Lines-C_Code_Generation_new_way.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",
    "# Generating C Code to implement Method of Lines Timestepping for Explicit Runge Kutta Methods\n",
    "\n",
    "## Authors: Zach Etienne & Brandon Clark\n",
    "\n",
    "## This tutorial notebook generates three blocks of C Code in order to perform Method of Lines timestepping. \n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). All Runge-Kutta (RK) Butcher tables were validated using truncated Taylor series in [a separate module](Tutorial-RK_Butcher_Table_Validation.ipynb). Finally, C-code implementation of RK4 was validated against a trusted version. C-code implementations of other RK methods seem to work as expected in the context of solving the scalar wave equation in Cartesian coordinates.\n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [MoLtimestepping/C_Code_Generation.py](../edit/MoLtimestepping/C_Code_Generation.py)\n",
    "* [MoLtimestepping/RK_Butcher_Table_Dictionary.py](../edit/MoLtimestepping/RK_Butcher_Table_Dictionary.py) ([**Tutorial**](Tutorial-RK_Butcher_Table_Dictionary.ipynb)) Stores the Butcher tables for the explicit Runge Kutta methods\n",
    "\n",
    "## Introduction:\n",
    "\n",
    "When numerically solving a partial differential equation initial-value problem, subject to suitable boundary conditions,  we implement Method of Lines to \"integrate\" the solution forward in time.\n",
    "\n",
    "\n",
    "### The Method of Lines:\n",
    "\n",
    "Once we have the initial data for a PDE, we \"evolve it forward in time\", using the [Method of Lines](https://reference.wolfram.com/language/tutorial/NDSolveMethodOfLines.html). In short, the Method of Lines enables us to handle \n",
    "1. the **spatial derivatives** of an initial value problem PDE using **standard finite difference approaches**, and\n",
    "2. the **temporal derivatives** of an initial value problem PDE using **standard strategies for solving ordinary differential equations (ODEs), like Runge Kutta methods** so long as the initial value problem PDE can be written in the first-order-in-time form\n",
    "$$\\partial_t \\vec{f} = \\mathbf{M}\\ \\vec{f},$$\n",
    "where $\\mathbf{M}$ is an $N\\times N$ matrix containing only *spatial* differential operators that act on the $N$-element column vector $\\vec{f}$. $\\mathbf{M}$ may not contain $t$ or time derivatives explicitly; only *spatial* partial derivatives are allowed to appear inside $\\mathbf{M}$.\n",
    "\n",
    "You may find the next module [Tutorial-ScalarWave](Tutorial-ScalarWave.ipynb) extremely helpful as an example for implementing the Method of Lines for solving the Scalar Wave equation in Cartesian coordinates.\n",
    "\n",
    "### Generating the C code:\n",
    "\n",
    "This module describes how core C functions are generated to implement Method of Lines timestepping for a specified RK method. There are three core functions:\n",
    "\n",
    "1. Allocate memory for gridfunctions.\n",
    "1. Step forward the solution one full timestep.\n",
    "1. Free memory for gridfunctions.\n",
    "\n",
    "The first function is called first, then the second function is repeated within a loop to a fixed \"final\" time (such that the end state of each iteration is the initial state for the next iteration), and the third function is called at the end of the calculation.\n",
    "\n",
    "The generated codes are essential for a number of Start-to-Finish example tutorial notebooks that demonstrate how to numerically solve hyperbolic PDEs."
   ]
  },
  {
   "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): Initialize needed Python/NRPy+ modules\n",
    "1. [Step 2](#diagonal): Checking if Butcher Table is Diagonal\n",
    "1. [Step 3](#ccode): Generating the C Code\n",
    "    1. [Step 3.a](#generategfnames): `generate_gridfunction_names()`: Uniquely and descriptively assign names to sets of gridfunctions\n",
    "    1. [Step 3.b](#alloc): Memory allocation: `MoL_malloc_y_n_gfs()` and `MoL_malloc_non_y_n_gfs()`\n",
    "    1. [Step 3.c](#molstep): Take one Method of Lines time step:  `MoL_step_forward_in_time()`\n",
    "    1. [Step 3.d](#free): Memory deallocation: `MoL_free_memory()`\n",
    "    1. [Step 3.e](#nrpybasicdefines): Define & register `MoL_gridfunctions_struct` in `NRPy_basic_defines.h`: `NRPy_basic_defines_MoL_timestepping_struct()`\n",
    "    1. [Step 3.f](#setupall): Add all MoL C codes to C function dictionary, and add MoL definitions to `NRPy_basic_defines.h`: `register_C_functions_and_NRPy_basic_defines()`\n",
    "1. [Step 4](#code_validation): Code Validation against `MoLtimestepping.RK_Butcher_Table_Generating_C_Code` NRPy+ module \n",
    "1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Initialize needed Python/NRPy+ modules [Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "Let's start by importing all the needed modules from Python/NRPy+:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:00.864214Z",
     "iopub.status.busy": "2021-03-07T17:15:00.862502Z",
     "iopub.status.idle": "2021-03-07T17:15:01.233098Z",
     "shell.execute_reply": "2021-03-07T17:15:01.233622Z"
    }
   },
   "outputs": [],
   "source": [
    "import sympy as sp  # Import SymPy, a computer algebra system written entirely in Python\n",
    "import os, sys      # Standard Python modules for multiplatform OS-level functions\n",
    "from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict\n",
    "from outputC import add_to_Cfunction_dict, indent_Ccode, outC_NRPy_basic_defines_h_dict, superfast_uniq, outputC  # NRPy+: Basic C code output functionality"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='diagonal'></a>\n",
    "\n",
    "# Step 2: Checking if a Butcher table is Diagonal [Back to [top](#toc)\\]\n",
    "$$\\label{diagonal}$$\n",
    "\n",
    "A diagonal Butcher table takes the form \n",
    "\n",
    "$$\\begin{array}{c|cccccc}\n",
    "    0 & \\\\\n",
    "    a_1 & a_1 & \\\\ \n",
    "    a_2 & 0 & a_2 & \\\\\n",
    "    a_3 & 0 & 0 & a_3 & \\\\ \n",
    "    \\vdots & \\vdots & \\ddots & \\ddots & \\ddots \\\\ \n",
    "    a_s & 0 & 0 & 0 & \\cdots & a_s \\\\ \\hline\n",
    "     & b_1 & b_2 & b_3 & \\cdots & b_{s-1} & b_s\n",
    "\\end{array}$$\n",
    "\n",
    "where $s$ is the number of required predictor-corrector steps for a given RK method (see [Butcher, John C. (2008)](https://onlinelibrary.wiley.com/doi/book/10.1002/9780470753767)). One known diagonal RK method is the classic RK4 represented in Butcher table form as:\n",
    "\n",
    "$$\\begin{array}{c|cccc}\n",
    "    0 & \\\\\n",
    "    1/2 & 1/2 & \\\\ \n",
    "    1/2 & 0 & 1/2 & \\\\\n",
    "    1 & 0 & 0 & 1 & \\\\ \\hline\n",
    "     & 1/6 & 1/3 & 1/3 & 1/6\n",
    "\\end{array} $$\n",
    "\n",
    "Diagonal Butcher tables are nice when it comes to saving required memory space. Each new step for a diagonal RK method, when computing the new $k_i$, does not depend on the previous calculation, and so there are ways to save memory. Significantly so in large three-dimensional spatial grid spaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:01.275215Z",
     "iopub.status.busy": "2021-03-07T17:15:01.274055Z",
     "iopub.status.idle": "2021-03-07T17:15:01.278771Z",
     "shell.execute_reply": "2021-03-07T17:15:01.279247Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The MoL method Euler is diagonal!\n",
      "The MoL method RK2 Heun is diagonal!\n",
      "The MoL method RK2 MP is diagonal!\n",
      "The MoL method RK2 Ralston is diagonal!\n",
      "The MoL method RK3 is NOT diagonal!\n",
      "The MoL method RK3 Heun is diagonal!\n",
      "The MoL method RK3 Ralston is diagonal!\n",
      "The MoL method SSPRK3 is NOT diagonal!\n",
      "The MoL method RK4 is diagonal!\n",
      "The MoL method DP5 is NOT diagonal!\n",
      "The MoL method DP5alt is NOT diagonal!\n",
      "The MoL method CK5 is NOT diagonal!\n",
      "The MoL method DP6 is NOT diagonal!\n",
      "The MoL method L6 is NOT diagonal!\n",
      "The MoL method DP8 is NOT diagonal!\n"
     ]
    }
   ],
   "source": [
    "# Check if Butcher Table is diagonal\n",
    "def diagonal(key):\n",
    "    Butcher = Butcher_dict[key][0]\n",
    "    L = len(Butcher)-1  # Establish the number of rows to check for diagonal trait, all bust last row\n",
    "    row_idx = 0  # Initialize the Butcher table row index\n",
    "    for i in range(L):  # Check all the desired rows\n",
    "        for j in range(1,row_idx):  # Check each element before the diagonal element in a row\n",
    "            if Butcher[i][j] != sp.sympify(0):  # If any non-diagonal coeffcient is non-zero,\n",
    "                                                # then the table is not diagonal\n",
    "                return False\n",
    "        row_idx += 1  # Update to check the next row\n",
    "    return True\n",
    "\n",
    "# Loop over all Butcher tables to check whether each is diagonal or not\n",
    "for key, value in Butcher_dict.items():\n",
    "    if diagonal(key) == True:\n",
    "        print(\"The MoL method \"+str(key)+\" is diagonal!\")\n",
    "    else:\n",
    "        print(\"The MoL method \"+str(key)+\" is NOT diagonal!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ccode'></a>\n",
    "\n",
    "# Step 3: Generating the C Code [Back to [top](#toc)\\]\n",
    "$$\\label{ccode}$$\n",
    "\n",
    "The following sections build up the C code for implementing the [Method of Lines timestepping algorithm](http://www.scholarpedia.org/article/Method_of_lines) for solving hyperbolic PDEs.\n",
    "\n",
    "**First an important note on efficiency:**\n",
    "\n",
    "Memory efficiency is incredibly important here, as $\\vec{f}$ is usually the largest object in memory.\n",
    "\n",
    "If we only made use of the Butcher tables without concern for memory efficiency, `generate_gridfunction_names()` and `MoL_step_forward_in_time()` would be very simple functions. \n",
    "\n",
    "It turns out that several of the Runge-Kutta-like methods in MoL can be made more efficient; for example \"RK4\" can be performed using only 4 \"timelevels\" of $\\vec{f}$ in memory (i.e., a total memory usage of `sizeof(f) * 4`). A naive implementation might use 5 or 6 copies. RK-like methods that have diagonal Butcher tables can be made far more efficient than the naive approach.\n",
    "\n",
    "**Exercise to student:** Improve the efficiency of other RK-like methods."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='generategfnames'></a>\n",
    "\n",
    "## Step 3.a: `generate_gridfunction_names()`: Uniquely and descriptively assign names to sets of gridfunctions [Back to [top](#toc)\\]\n",
    "$$\\label{generategfnames}$$\n",
    "\n",
    "`generate_gridfunction_names()` names gridfunctions to be consistent with a given RK substep. For example we might call the set of gridfunctions stored at substep $k_1$ `k1_gfs`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:01.333399Z",
     "iopub.status.busy": "2021-03-07T17:15:01.327431Z",
     "iopub.status.idle": "2021-03-07T17:15:01.357384Z",
     "shell.execute_reply": "2021-03-07T17:15:01.356760Z"
    }
   },
   "outputs": [],
   "source": [
    "# Each MoL method has its own set of names for groups of gridfunctions,\n",
    "#   aiming to be sufficiently descriptive. So for example a set of\n",
    "#   gridfunctions that store \"k_1\" in an RK-like method could be called\n",
    "#   \"k1_gfs\".\n",
    "def generate_gridfunction_names(MoL_method = \"RK4\"):\n",
    "    # Step 3.a: MoL gridfunctions fall into 3 overlapping categories:\n",
    "    #           1) y_n=y_i(t_n) gridfunctions y_n_gfs, which stores data for the vector of gridfunctions y_i at t_n,\n",
    "    #              the start of each MoL timestep.\n",
    "    #           2) non-y_n gridfunctions, needed to compute the data at t_{n+1}. Often labeled with k_i in the name,\n",
    "    #              these gridfunctions are *not* needed at the start of each timestep, so are available for temporary\n",
    "    #              storage when gridfunctions needed for diagnostics are computed at the start of each timestep.\n",
    "    #              These gridfunctions can also be freed during a regrid, to enable storage for the post-regrid\n",
    "    #              destination y_n_gfs.\n",
    "    #           3) Diagnostic output gridfunctions diagnostic_output_gfs, which simply uses the memory from auxiliary\n",
    "    #              gridfunctions at one auxiliary time to compute diagnostics at t_n.\n",
    "\n",
    "    # Here we specify which gridfunctions fall into each category, starting with the obvious: y_n_gridfunctions\n",
    "    y_n_gridfunctions = \"y_n_gfs\"\n",
    "\n",
    "    # Next the less-obvious, which depend on non-y_n_gfs\n",
    "    non_y_n_gridfunctions_list = []\n",
    "\n",
    "    # No matter the method we define gridfunctions \"y_n_gfs\" to store the initial data\n",
    "    if diagonal(MoL_method) and \"RK3\" in MoL_method:\n",
    "        non_y_n_gridfunctions_list.append(\"k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs\")\n",
    "        non_y_n_gridfunctions_list.append(\"k2_or_y_nplus_a32_k2_gfs\")\n",
    "        diagnostic_gridfunctions_point_to = \"k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs\"\n",
    "    else:\n",
    "        if not diagonal(MoL_method):  # Allocate memory for non-diagonal Butcher tables\n",
    "            # Determine the number of k_i steps based on length of Butcher Table\n",
    "            num_k = len(Butcher_dict[MoL_method][0])-1\n",
    "            # For non-diagonal tables an intermediate gridfunction \"next_y_input\" is used for rhs evaluations\n",
    "            non_y_n_gridfunctions_list.append(\"next_y_input_gfs\")\n",
    "            for i in range(num_k): # Need to allocate all k_i steps for a given method\n",
    "                non_y_n_gridfunctions_list.append(\"k\" + str(i + 1) + \"_gfs\")\n",
    "            diagnostic_gridfunctions_point_to = \"k1_gfs\"\n",
    "        else:  # Allocate memory for diagonal Butcher tables, which use a \"y_nplus1_running_total gridfunction\"\n",
    "            non_y_n_gridfunctions_list.append(\"y_nplus1_running_total_gfs\")\n",
    "            if MoL_method != 'Euler':  # Allocate memory for diagonal Butcher tables that aren't Euler\n",
    "                # Need k_odd for k_1,3,5... and k_even for k_2,4,6...\n",
    "                non_y_n_gridfunctions_list.append(\"k_odd_gfs\")\n",
    "                non_y_n_gridfunctions_list.append(\"k_even_gfs\")\n",
    "            diagnostic_gridfunctions_point_to = \"y_nplus1_running_total_gfs\"\n",
    "    non_y_n_gridfunctions_list.append(\"auxevol_gfs\")\n",
    "\n",
    "    return y_n_gridfunctions, non_y_n_gridfunctions_list, diagnostic_gridfunctions_point_to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='alloc'></a>\n",
    "\n",
    "## Step 3.b: Memory allocation: `MoL_malloc_y_n_gfs()` and `MoL_malloc_non_y_n_gfs()`:  [Back to [top](#toc)\\]\n",
    "$$\\label{alloc}$$\n",
    "\n",
    "Generation of C functions `MoL_malloc_y_n_gfs()` and `MoL_malloc_non_y_n_gfs()` read the full list of needed lists of gridfunctions, provided by (Python) function `generate_gridfunction_names()`, to allocate space for all gridfunctions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:01.333399Z",
     "iopub.status.busy": "2021-03-07T17:15:01.327431Z",
     "iopub.status.idle": "2021-03-07T17:15:01.357384Z",
     "shell.execute_reply": "2021-03-07T17:15:01.356760Z"
    }
   },
   "outputs": [],
   "source": [
    "# add_to_Cfunction_dict_MoL_malloc() registers\n",
    "#           MoL_malloc_y_n_gfs() and\n",
    "#           MoL_malloc_non_y_n_gfs(), which allocate memory for\n",
    "#           the indicated sets of gridfunctions\n",
    "def add_to_Cfunction_dict_MoL_malloc(MoL_method, which_gfs):\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    desc  = \"Method of Lines (MoL) for \\\"\" + MoL_method + \"\\\" method: Allocate memory for \\\"\"+which_gfs+\"\\\" gridfunctions\\n\"\n",
    "    desc += \"   * y_n_gfs are used to store data for the vector of gridfunctions y_i at t_n, at the start of each MoL timestep\\n\"\n",
    "    desc += \"   * non_y_n_gfs are needed for intermediate (e.g., k_i) storage in chosen MoL method\\n\"\n",
    "    c_type = \"void\"\n",
    "\n",
    "    y_n_gridfunctions, non_y_n_gridfunctions_list, diagnostic_gridfunctions_point_to = \\\n",
    "        generate_gridfunction_names(MoL_method = MoL_method)\n",
    "\n",
    "    gridfunctions_list = []\n",
    "    if which_gfs == \"y_n_gfs\":\n",
    "        gridfunctions_list = [y_n_gridfunctions]\n",
    "    elif which_gfs == \"non_y_n_gfs\":\n",
    "        gridfunctions_list = non_y_n_gridfunctions_list\n",
    "    else:\n",
    "        print(\"ERROR: which_gfs = \\\"\" + which_gfs + \"\\\" unrecognized.\")\n",
    "        sys.exit(1)\n",
    "    name = \"MoL_malloc_\" + which_gfs\n",
    "    params = \"const paramstruct *restrict params, MoL_gridfunctions_struct *restrict gridfuncs\"\n",
    "    body = \"const int Nxx_plus_2NGHOSTS_tot = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\\n\"\n",
    "    for gridfunctions in gridfunctions_list:\n",
    "        num_gfs = \"NUM_EVOL_GFS\"\n",
    "        if gridfunctions == \"auxevol_gfs\":\n",
    "            num_gfs = \"NUM_AUXEVOL_GFS\"\n",
    "        body += \"gridfuncs->\" + gridfunctions + \" = (REAL *restrict)malloc(sizeof(REAL) * \" + num_gfs + \" * Nxx_plus_2NGHOSTS_tot);\\n\"\n",
    "    body += \"\\ngridfuncs->diagnostic_output_gfs = gridfuncs->\" + diagnostic_gridfunctions_point_to + \";\\n\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=indent_Ccode(body, \"  \"),\n",
    "        rel_path_to_Cparams=os.path.join(\".\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='molstep'></a>\n",
    "\n",
    "## Step 3.c: Take one Method of Lines time step:  `MoL_step_forward_in_time()` [Back to [top](#toc)\\]\n",
    "$$\\label{molstep}$$\n",
    "\n",
    "An MoL step consists in general of a series of Runge-Kutta-like substeps, and the `MoL_step_forward_in_time()` C function pulls together all of these substeps.\n",
    "\n",
    "The basic C code for an MoL substep, set up by the Python function `single_RK_substep()` below, is as follows.\n",
    "\n",
    "1. Evaluate the right-hand side of $\\partial_t \\vec{f}=$ `RHS`, to get the time derivative of the set of gridfunctions $\\vec{f}$ at our current time.\n",
    "1. Perform the Runge-Kutta update, which depends on $\\partial_t \\vec{f}$ on the current and sometimes previous times.\n",
    "1. Call post-right-hand side functions as desired.\n",
    "\n",
    "The `single_RK_substep_input_symbolic()` function generates the C code for performing the above steps, applying substitutions for e.g., `RK_INPUT_GFS` and `RK_OUTPUT_GFS` as appropriate. `single_RK_substep_input_symbolic()` supports SIMD-capable code generation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# single_RK_substep_input_symbolic() performs necessary replacements to\n",
    "#   define C code for a single RK substep\n",
    "#   (e.g., computing k_1 and then updating the outer boundaries)\n",
    "def single_RK_substep_input_symbolic(commentblock, RHS_str, RHS_input_str, RHS_output_str, RK_lhss_list, RK_rhss_list,\n",
    "                                     post_RHS_list, post_RHS_output_list, enable_SIMD=False,\n",
    "                                     enable_griddata=False, gf_aliases=\"\", post_post_RHS_string=\"\"):\n",
    "    return_str  = commentblock + \"\\n\"\n",
    "    if not isinstance(RK_lhss_list, list):\n",
    "        RK_lhss_list = [RK_lhss_list]\n",
    "    if not isinstance(RK_rhss_list, list):\n",
    "        RK_rhss_list = [RK_rhss_list]\n",
    "\n",
    "    if not isinstance(post_RHS_list, list):\n",
    "        post_RHS_list = [post_RHS_list]\n",
    "    if not isinstance(post_RHS_output_list, list):\n",
    "        post_RHS_output_list = [post_RHS_output_list]\n",
    "\n",
    "    indent = \"\"\n",
    "    if enable_griddata:\n",
    "        return_str += \"{\\n\" + indent_Ccode(gf_aliases, \"  \")\n",
    "        indent = \"  \"\n",
    "    # Part 1: RHS evaluation:\n",
    "    return_str += indent_Ccode(str(RHS_str).replace(\"RK_INPUT_GFS\",  str(RHS_input_str).replace(\"gfsL\", \"gfs\")).\n",
    "                               replace(\"RK_OUTPUT_GFS\", str(RHS_output_str).replace(\"gfsL\", \"gfs\"))+\"\\n\", indent=indent)\n",
    "\n",
    "    # Part 2: RK update\n",
    "    if enable_SIMD:\n",
    "        return_str += \"#pragma omp parallel for\\n\"\n",
    "        return_str += indent + \"for(int i=0;i<Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2*NUM_EVOL_GFS;i+=SIMD_width) {\\n\"\n",
    "    else:\n",
    "        return_str += indent + \"LOOP_ALL_GFS_GPS(i) {\\n\"\n",
    "    type = \"REAL\"\n",
    "    if enable_SIMD:\n",
    "        type = \"REAL_SIMD_ARRAY\"\n",
    "    RK_lhss_str_list = []\n",
    "    for i, el in enumerate(RK_lhss_list):\n",
    "        if enable_SIMD:\n",
    "            RK_lhss_str_list.append(indent + \"const REAL_SIMD_ARRAY __RHS_exp_\" + str(i))\n",
    "        else:\n",
    "            RK_lhss_str_list.append(indent + str(el).replace(\"gfsL\", \"gfs[i]\"))\n",
    "    read_list = []\n",
    "    for el in RK_rhss_list:\n",
    "        for read in list(sp.ordered(el.free_symbols)):\n",
    "            read_list.append(read)\n",
    "    read_list_uniq = superfast_uniq(read_list)\n",
    "    for el in read_list_uniq:\n",
    "        if str(el) != \"dt\":\n",
    "            if enable_SIMD:\n",
    "                return_str += indent + \"  const \" + type + \" \" + str(el) + \" = ReadSIMD(&\" + str(el).replace(\"gfsL\", \"gfs[i]\") + \");\\n\"\n",
    "            else:\n",
    "                return_str += indent + \"  const \" + type + \" \" + str(el) + \" = \" + str(el).replace(\"gfsL\", \"gfs[i]\") + \";\\n\"\n",
    "    if enable_SIMD:\n",
    "        return_str += indent + \"  const REAL_SIMD_ARRAY DT = ConstSIMD(dt);\\n\"\n",
    "    preindent = \"1\"\n",
    "    if enable_griddata:\n",
    "        preindent = \"2\"\n",
    "    kernel = outputC(RK_rhss_list, RK_lhss_str_list, filename=\"returnstring\",\n",
    "                     params=\"includebraces=False,preindent=\"+preindent+\",outCverbose=False,enable_SIMD=\"+str(enable_SIMD))\n",
    "    if enable_SIMD:\n",
    "        return_str += kernel.replace(\"dt\", \"DT\")\n",
    "        for i, el in enumerate(RK_lhss_list):\n",
    "            return_str += \"  WriteSIMD(&\" + str(el).replace(\"gfsL\", \"gfs[i]\") + \", __RHS_exp_\" + str(i) + \");\\n\"\n",
    "    else:\n",
    "        return_str += kernel\n",
    "\n",
    "    return_str += indent + \"}\\n\"\n",
    "\n",
    "    # Part 3: Call post-RHS functions\n",
    "    for post_RHS, post_RHS_output in zip(post_RHS_list, post_RHS_output_list):\n",
    "        return_str += indent_Ccode(post_RHS.replace(\"RK_OUTPUT_GFS\", str(post_RHS_output).replace(\"gfsL\", \"gfs\")))\n",
    "\n",
    "    if enable_griddata:\n",
    "        return_str += \"}\\n\"\n",
    "\n",
    "    for post_RHS, post_RHS_output in zip(post_RHS_list, post_RHS_output_list):\n",
    "        return_str += indent_Ccode(post_post_RHS_string.replace(\"RK_OUTPUT_GFS\", str(post_RHS_output).replace(\"gfsL\", \"gfs\")), \"\")\n",
    "\n",
    "    return return_str"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the `add_to_Cfunction_dict_MoL_step_forward_in_time()` Python function below we construct and register the core C function for MoL timestepping: `MoL_step_forward_in_time()`. `MoL_step_forward_in_time()` implements Butcher tables for Runge-Kutta-like methods, leveraging the `single_RK_substep()` helper function above as needed. Again, we aim for maximum memory efficiency so that e.g., RK4 needs to store only 4 levels of $\\vec{f}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:01.333399Z",
     "iopub.status.busy": "2021-03-07T17:15:01.327431Z",
     "iopub.status.idle": "2021-03-07T17:15:01.357384Z",
     "shell.execute_reply": "2021-03-07T17:15:01.356760Z"
    }
   },
   "outputs": [],
   "source": [
    "def add_to_Cfunction_dict_MoL_step_forward_in_time(MoL_method,\n",
    "                                                   RHS_string = \"\", post_RHS_string = \"\", post_post_RHS_string=\"\",\n",
    "                                                   enable_rfm=False, enable_curviBCs=False, enable_SIMD=False,\n",
    "                                                   enable_griddata=False):\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    if enable_SIMD:\n",
    "        includes += [os.path.join(\"SIMD\", \"SIMD_intrinsics.h\")]\n",
    "    desc  = \"Method of Lines (MoL) for \\\"\" + MoL_method + \"\\\" method: Step forward one full timestep.\\n\"\n",
    "    c_type = \"void\"\n",
    "    name = \"MoL_step_forward_in_time\"\n",
    "    if enable_griddata:\n",
    "        params = \"griddata_struct *restrict griddata, const REAL dt\"\n",
    "    else:\n",
    "        params  = \"const paramstruct *restrict params, \"\n",
    "        if enable_rfm:\n",
    "            params += \"const rfm_struct *restrict rfmstruct, \"\n",
    "        else:\n",
    "            params += \"REAL *restrict xx[3], \"\n",
    "        if enable_curviBCs:\n",
    "            params += \"const bc_struct *restrict bcstruct, \"\n",
    "        params += \"MoL_gridfunctions_struct *restrict gridfuncs, const REAL dt\"\n",
    "\n",
    "    indent = \"\"  # We don't bother with an indent here.\n",
    "\n",
    "    body = indent + \"// C code implementation of -={ \" + MoL_method + \" }=- Method of Lines timestepping.\\n\\n\"\n",
    "\n",
    "    y_n_gridfunctions, non_y_n_gridfunctions_list, _throwaway = generate_gridfunction_names(MoL_method)\n",
    "    if enable_griddata:\n",
    "        gf_prefix = \"griddata->gridfuncs.\"\n",
    "    else:\n",
    "        gf_prefix = \"gridfuncs->\"\n",
    "    gf_aliases = \"\"\"// Set gridfunction aliases from gridfuncs struct\n",
    "REAL *restrict \"\"\" + y_n_gridfunctions + \" = \"+gf_prefix + y_n_gridfunctions + \"\"\";  // y_n gridfunctions\n",
    "// Temporary timelevel & AUXEVOL gridfunctions:\\n\"\"\"\n",
    "    for gf in non_y_n_gridfunctions_list:\n",
    "        gf_aliases += \"REAL *restrict \" + gf + \" = \"+gf_prefix + gf + \";\\n\"\n",
    "    if enable_griddata:\n",
    "        gf_aliases += \"paramstruct *restrict params = &griddata->params;\\n\"\n",
    "        gf_aliases += \"const rfm_struct *restrict rfmstruct = &griddata->rfmstruct;\\n\"\n",
    "        gf_aliases += \"const bc_struct *restrict bcstruct = &griddata->bcstruct;\\n\"\n",
    "        for i in [\"0\", \"1\", \"2\"]:\n",
    "            gf_aliases += \"const int Nxx_plus_2NGHOSTS\" + i + \" = griddata->params.Nxx_plus_2NGHOSTS\" + i + \";\\n\"\n",
    "\n",
    "    if not enable_griddata:\n",
    "        body += gf_aliases\n",
    "\n",
    "    # Implement Method of Lines (MoL) Timestepping\n",
    "    Butcher = Butcher_dict[MoL_method][0]  # Get the desired Butcher table from the dictionary\n",
    "    num_steps = len(Butcher)-1  # Specify the number of required steps to update solution\n",
    "\n",
    "    # Diagonal RK3 only!!!\n",
    "    dt = sp.Symbol(\"dt\", real=True)\n",
    "    if diagonal(MoL_method) and \"RK3\" in MoL_method:\n",
    "        #  In a diagonal RK3 method, only 3 gridfunctions need be defined. Below implements this approach.\n",
    "        y_n_gfs = sp.Symbol(\"y_n_gfsL\", real=True)\n",
    "        k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs = sp.Symbol(\"k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfsL\", real=True)\n",
    "        k2_or_y_nplus_a32_k2_gfs = sp.Symbol(\"k2_or_y_nplus_a32_k2_gfsL\", real=True)\n",
    "        # k_1\n",
    "        body += \"\"\"\n",
    "    // In a diagonal RK3 method like this one, only 3 gridfunctions need be defined. Below implements this approach.\n",
    "    // Using y_n_gfs as input, k1 and apply boundary conditions\\n\"\"\"\n",
    "        body += single_RK_substep_input_symbolic(\n",
    "            commentblock=\"\"\"// -={ START k1 substep }=-\n",
    "    // RHS evaluation:\n",
    "    //  1. We will store k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs now as\n",
    "    //     ...  the update for the next rhs evaluation y_n + a21*k1*dt\n",
    "    // Post-RHS evaluation:\n",
    "    //  1. Apply post-RHS to y_n + a21*k1*dt\"\"\",\n",
    "            RHS_str=RHS_string,\n",
    "            RHS_input_str=y_n_gfs,\n",
    "            RHS_output_str=k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs,\n",
    "            RK_lhss_list=[k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs],\n",
    "            RK_rhss_list=[Butcher[1][1]*k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs*dt + y_n_gfs],\n",
    "            post_RHS_list=[post_RHS_string], post_RHS_output_list=[k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs],\n",
    "            enable_SIMD=enable_SIMD, enable_griddata=enable_griddata, gf_aliases=gf_aliases,\n",
    "            post_post_RHS_string=post_post_RHS_string) + \"// -={ END k1 substep }=-\\n\\n\"\n",
    "\n",
    "        # k_2\n",
    "        body += single_RK_substep_input_symbolic(\n",
    "            commentblock=\"\"\"// -={ START k2 substep }=-\n",
    "    // RHS evaluation:\n",
    "    //    1. Reassign k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs to be the running total y_{n+1}; a32*k2*dt to the running total\n",
    "    //    2. Store k2_or_y_nplus_a32_k2_gfs now as y_n + a32*k2*dt\n",
    "    // Post-RHS evaluation:\n",
    "    //    1. Apply post-RHS to both y_n + a32*k2 (stored in k2_or_y_nplus_a32_k2_gfs)\n",
    "    //       ... and the y_{n+1} running total, as they have not been applied yet to k2-related gridfunctions\"\"\",\n",
    "            RHS_str=RHS_string,\n",
    "            RHS_input_str=k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs,\n",
    "            RHS_output_str=k2_or_y_nplus_a32_k2_gfs,\n",
    "            RK_lhss_list=[k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs, k2_or_y_nplus_a32_k2_gfs],\n",
    "            RK_rhss_list=[Butcher[3][1]*(k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs - y_n_gfs)/Butcher[1][1] + y_n_gfs + Butcher[3][2]*k2_or_y_nplus_a32_k2_gfs*dt,\n",
    "                          Butcher[2][2]*k2_or_y_nplus_a32_k2_gfs*dt + y_n_gfs],\n",
    "            post_RHS_list=[post_RHS_string, post_RHS_string],\n",
    "            post_RHS_output_list=[k2_or_y_nplus_a32_k2_gfs,\n",
    "                                  k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs],\n",
    "            enable_SIMD=enable_SIMD, enable_griddata=enable_griddata, gf_aliases=gf_aliases,\n",
    "            post_post_RHS_string=post_post_RHS_string) + \"// -={ END k2 substep }=-\\n\\n\"\n",
    "\n",
    "        # k_3\n",
    "        body += single_RK_substep_input_symbolic(\n",
    "            commentblock=\"\"\"// -={ START k3 substep }=-\n",
    "        // RHS evaluation:\n",
    "        //    1. Add k3 to the running total and save to y_n\n",
    "        // Post-RHS evaluation:\n",
    "        //    1. Apply post-RHS to y_n\"\"\",\n",
    "            RHS_str=RHS_string,\n",
    "            RHS_input_str=k2_or_y_nplus_a32_k2_gfs, RHS_output_str=y_n_gfs,\n",
    "            RK_lhss_list=[y_n_gfs],\n",
    "            RK_rhss_list=[k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs + Butcher[3][3]*y_n_gfs*dt],\n",
    "            post_RHS_list=[post_RHS_string],\n",
    "            post_RHS_output_list=[y_n_gfs],\n",
    "            enable_SIMD=enable_SIMD, enable_griddata=enable_griddata, gf_aliases=gf_aliases,\n",
    "            post_post_RHS_string=post_post_RHS_string) + \"// -={ END k3 substep }=-\\n\\n\"\n",
    "    else:\n",
    "        y_n = sp.Symbol(\"y_n_gfsL\", real=True)\n",
    "        if not diagonal(MoL_method):\n",
    "            for s in range(num_steps):\n",
    "                next_y_input = sp.Symbol(\"next_y_input_gfsL\", real=True)\n",
    "\n",
    "                # If we're on the first step (s=0), we use y_n gridfunction as input.\n",
    "                #      Otherwise next_y_input is input. Output is just the reverse.\n",
    "                if s == 0:  # If on first step:\n",
    "                    RHS_input = y_n\n",
    "                else:  # If on second step or later:\n",
    "                    RHS_input = next_y_input\n",
    "                RHS_output = sp.Symbol(\"k\" + str(s + 1) + \"_gfs\", real=True)\n",
    "                if s == num_steps - 1:  # If on final step:\n",
    "                    RK_lhs = y_n\n",
    "                else:  # If on anything but the final step:\n",
    "                    RK_lhs = next_y_input\n",
    "                RK_rhs = y_n\n",
    "                for m in range(s + 1):\n",
    "                    k_mp1_gfs = sp.Symbol(\"k\" + str(m + 1) + \"_gfsL\")\n",
    "                    if Butcher[s + 1][m + 1] != 0:\n",
    "                        if Butcher[s + 1][m + 1] != 1:\n",
    "                            RK_rhs += dt * k_mp1_gfs*Butcher[s + 1][m + 1]\n",
    "                        else:\n",
    "                            RK_rhs += dt * k_mp1_gfs\n",
    "\n",
    "                post_RHS = post_RHS_string\n",
    "                if s == num_steps - 1:  # If on final step:\n",
    "                    post_RHS_output = y_n\n",
    "                else:  # If on anything but the final step:\n",
    "                    post_RHS_output = next_y_input\n",
    "\n",
    "                body += single_RK_substep_input_symbolic(\n",
    "                    commentblock=\"// -={ START k\" + str(s + 1) + \" substep }=-\",\n",
    "                    RHS_str=RHS_string,\n",
    "                    RHS_input_str=RHS_input, RHS_output_str=RHS_output,\n",
    "                    RK_lhss_list=[RK_lhs], RK_rhss_list=[RK_rhs],\n",
    "                    post_RHS_list=[post_RHS],\n",
    "                    post_RHS_output_list=[post_RHS_output],\n",
    "                    enable_SIMD=enable_SIMD, enable_griddata=enable_griddata, gf_aliases=gf_aliases,\n",
    "                    post_post_RHS_string=post_post_RHS_string) + \"// -={ END k\" + str(s + 1) + \" substep }=-\\n\\n\"\n",
    "        else:\n",
    "            y_n = sp.Symbol(\"y_n_gfsL\", real=True)\n",
    "            y_nplus1_running_total = sp.Symbol(\"y_nplus1_running_total_gfsL\", real=True)\n",
    "            if MoL_method == 'Euler':  # Euler's method doesn't require any k_i, and gets its own unique algorithm\n",
    "                body += single_RK_substep_input_symbolic(\n",
    "                    commentblock=indent + \"// ***Euler timestepping only requires one RHS evaluation***\",\n",
    "                    RHS_str=RHS_string,\n",
    "                    RHS_input_str=y_n, RHS_output_str=y_nplus1_running_total,\n",
    "                    RK_lhss_list=[y_n], RK_rhss_list=[y_n + y_nplus1_running_total*dt],\n",
    "                    post_RHS_list=[post_RHS_string],\n",
    "                    post_RHS_output_list=[y_n],\n",
    "                    enable_SIMD=enable_SIMD, enable_griddata=enable_griddata, gf_aliases=gf_aliases,\n",
    "                    post_post_RHS_string=post_post_RHS_string)\n",
    "            else:\n",
    "                for s in range(num_steps):\n",
    "                    # If we're on the first step (s=0), we use y_n gridfunction as input.\n",
    "                    # and k_odd as output.\n",
    "                    if s == 0:\n",
    "                        RHS_input  = sp.Symbol(\"y_n_gfsL\", real=True)\n",
    "                        RHS_output = sp.Symbol(\"k_odd_gfsL\", real=True)\n",
    "                    # For the remaining steps the inputs and ouputs alternate between k_odd and k_even\n",
    "                    elif s % 2 == 0:\n",
    "                        RHS_input = sp.Symbol(\"k_even_gfsL\", real=True)\n",
    "                        RHS_output = sp.Symbol(\"k_odd_gfsL\", real=True)\n",
    "                    else:\n",
    "                        RHS_input = sp.Symbol(\"k_odd_gfsL\", real=True)\n",
    "                        RHS_output = sp.Symbol(\"k_even_gfsL\", real=True)\n",
    "                    RK_lhs_list = []\n",
    "                    RK_rhs_list = []\n",
    "                    if s != num_steps-1:  # For anything besides the final step\n",
    "                        if s == 0:  # The first RK step\n",
    "                            RK_lhs_list.append(y_nplus1_running_total)\n",
    "                            RK_rhs_list.append(RHS_output*dt*Butcher[num_steps][s+1])\n",
    "\n",
    "                            RK_lhs_list.append(RHS_output)\n",
    "                            RK_rhs_list.append(y_n + RHS_output*dt*Butcher[s+1][s+1])\n",
    "                        else:\n",
    "                            if Butcher[num_steps][s+1] != 0:\n",
    "                                RK_lhs_list.append(y_nplus1_running_total)\n",
    "                                if Butcher[num_steps][s+1] != 1:\n",
    "                                    RK_rhs_list.append(y_nplus1_running_total + RHS_output*dt*Butcher[num_steps][s+1])\n",
    "                                else:\n",
    "                                    RK_rhs_list.append(y_nplus1_running_total + RHS_output*dt)\n",
    "                            if Butcher[s+1][s+1] != 0:\n",
    "                                RK_lhs_list.append(RHS_output)\n",
    "                                if Butcher[s+1][s+1] != 1:\n",
    "                                    RK_rhs_list.append(y_n + RHS_output*dt*Butcher[s+1][s+1])\n",
    "                                else:\n",
    "                                    RK_rhs_list.append(y_n + RHS_output*dt)\n",
    "                        post_RHS_output = RHS_output\n",
    "                    if s == num_steps-1:  # If on the final step\n",
    "                        if Butcher[num_steps][s+1] != 0:\n",
    "                            RK_lhs_list.append(y_n)\n",
    "                            if Butcher[num_steps][s+1] != 1:\n",
    "                                RK_rhs_list.append(y_n + y_nplus1_running_total + RHS_output*dt*Butcher[num_steps][s+1])\n",
    "                            else:\n",
    "                                RK_rhs_list.append(y_n + y_nplus1_running_total + RHS_output*dt)\n",
    "                        post_RHS_output = y_n\n",
    "                    body += single_RK_substep_input_symbolic(\n",
    "                        commentblock=indent + \"// -={ START k\" + str(s + 1) + \" substep }=-\",\n",
    "                        RHS_str=RHS_string,\n",
    "                        RHS_input_str=RHS_input, RHS_output_str=RHS_output,\n",
    "                        RK_lhss_list=RK_lhs_list, RK_rhss_list=RK_rhs_list,\n",
    "                        post_RHS_list=[post_RHS_string],\n",
    "                        post_RHS_output_list=[post_RHS_output],\n",
    "                        enable_SIMD=enable_SIMD, enable_griddata=enable_griddata, gf_aliases=gf_aliases,\n",
    "                        post_post_RHS_string=post_post_RHS_string) + \"// -={ END k\" + str(s + 1) + \" substep }=-\\n\\n\"\n",
    "\n",
    "    enableCparameters=True\n",
    "    if enable_griddata:\n",
    "        enableCparameters=False\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=indent_Ccode(body, \"  \"),\n",
    "        enableCparameters=enableCparameters, rel_path_to_Cparams=os.path.join(\".\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='free'></a>\n",
    "\n",
    "## Step 3.d: Memory deallocation: `MoL_free_memory()` [Back to [top](#toc)\\]\n",
    "$$\\label{free}$$\n",
    "\n",
    "We define the function `MoL_free_memory()` which generates the C code for freeing the memory that was being occupied by the grid functions lists that had been allocated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:01.291671Z",
     "iopub.status.busy": "2021-03-07T17:15:01.290769Z",
     "iopub.status.idle": "2021-03-07T17:15:01.293649Z",
     "shell.execute_reply": "2021-03-07T17:15:01.293107Z"
    }
   },
   "outputs": [],
   "source": [
    "# add_to_Cfunction_dict_MoL_free_memory() registers\n",
    "#           MoL_free_memory_y_n_gfs() and\n",
    "#           MoL_free_memory_non_y_n_gfs(), which free memory for\n",
    "#           the indicated sets of gridfunctions\n",
    "def add_to_Cfunction_dict_MoL_free_memory(MoL_method, which_gfs):\n",
    "    includes = [\"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
    "    desc = \"Method of Lines (MoL) for \\\"\" + MoL_method + \"\\\" method: Free memory for \\\"\" + which_gfs + \"\\\" gridfunctions\\n\"\n",
    "    desc += \"   - y_n_gfs are used to store data for the vector of gridfunctions y_i at t_n, at the start of each MoL timestep\\n\"\n",
    "    desc += \"   - non_y_n_gfs are needed for intermediate (e.g., k_i) storage in chosen MoL method\\n\"\n",
    "    c_type = \"void\"\n",
    "\n",
    "    y_n_gridfunctions, non_y_n_gridfunctions_list, diagnostic_gridfunctions_point_to = \\\n",
    "        generate_gridfunction_names(MoL_method=MoL_method)\n",
    "\n",
    "    gridfunctions_list = []\n",
    "    if which_gfs == \"y_n_gfs\":\n",
    "        gridfunctions_list = [y_n_gridfunctions]\n",
    "    elif which_gfs == \"non_y_n_gfs\":\n",
    "        gridfunctions_list = non_y_n_gridfunctions_list\n",
    "    else:\n",
    "        print(\"ERROR: which_gfs = \\\"\" + which_gfs + \"\\\" unrecognized.\")\n",
    "        sys.exit(1)\n",
    "    name = \"MoL_free_memory_\" + which_gfs\n",
    "    params = \"const paramstruct *restrict params, MoL_gridfunctions_struct *restrict gridfuncs\"\n",
    "    body = \"\"\n",
    "    for gridfunctions in gridfunctions_list:\n",
    "        body += \"    free(gridfuncs->\" + gridfunctions + \");\\n\"\n",
    "    add_to_Cfunction_dict(\n",
    "        includes=includes,\n",
    "        desc=desc,\n",
    "        c_type=c_type, name=name, params=params,\n",
    "        body=indent_Ccode(body, \"  \"),\n",
    "        rel_path_to_Cparams=os.path.join(\".\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='nrpybasicdefines'></a>\n",
    "\n",
    "## Step 3.e: Define & register `MoL_gridfunctions_struct` in `NRPy_basic_defines.h`: `NRPy_basic_defines_MoL_timestepping_struct()` [Back to [top](#toc)\\]\n",
    "$$\\label{nrpybasicdefines}$$\n",
    "\n",
    "`MoL_gridfunctions_struct` stores pointers to all the gridfunctions needed by MoL, and we define this struct within `NRPy_basic_defines.h`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:01.333399Z",
     "iopub.status.busy": "2021-03-07T17:15:01.327431Z",
     "iopub.status.idle": "2021-03-07T17:15:01.357384Z",
     "shell.execute_reply": "2021-03-07T17:15:01.356760Z"
    }
   },
   "outputs": [],
   "source": [
    "# Register MoL_gridfunctions_struct in NRPy_basic_defines\n",
    "def NRPy_basic_defines_MoL_timestepping_struct(MoL_method=\"RK4\"):\n",
    "    y_n_gridfunctions, non_y_n_gridfunctions_list, diagnostic_gridfunctions_point_to = \\\n",
    "        generate_gridfunction_names(MoL_method=MoL_method)\n",
    "    # Step 3.b: Create MoL_timestepping struct:\n",
    "    indent = \"  \"\n",
    "    Nbd = \"typedef struct __MoL_gridfunctions_struct__ {\\n\"\n",
    "    Nbd += indent + \"REAL *restrict \" + y_n_gridfunctions + \";\\n\"\n",
    "    for gfs in non_y_n_gridfunctions_list:\n",
    "        Nbd += indent + \"REAL *restrict \" + gfs + \";\\n\"\n",
    "    Nbd += indent + \"REAL *restrict diagnostic_output_gfs;\\n\"\n",
    "    Nbd += \"} MoL_gridfunctions_struct;\\n\"\n",
    "    Nbd += \"\"\"#define LOOP_ALL_GFS_GPS(ii) _Pragma(\"omp parallel for\") \\\\\n",
    "  for(int (ii)=0;(ii)<Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2*NUM_EVOL_GFS;(ii)++)\\n\"\"\"\n",
    "\n",
    "    outC_NRPy_basic_defines_h_dict[\"MoL\"] = Nbd"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='setupall'></a>\n",
    "\n",
    "## Step 3.f: Add all MoL C codes to C function dictionary, and add MoL definitions to `NRPy_basic_defines.h`: `register_C_functions_and_NRPy_basic_defines()` \\[Back to [top](#toc)\\]\n",
    "$$\\label{setupall}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Finally declare the master registration function\n",
    "def register_C_functions_and_NRPy_basic_defines(MoL_method = \"RK4\",\n",
    "            RHS_string =  \"rhs_eval(Nxx,Nxx_plus_2NGHOSTS,dxx, RK_INPUT_GFS, RK_OUTPUT_GFS);\",\n",
    "            post_RHS_string = \"apply_bcs(Nxx,Nxx_plus_2NGHOSTS, RK_OUTPUT_GFS);\", post_post_RHS_string = \"\",\n",
    "            enable_rfm=False, enable_curviBCs=False, enable_SIMD=False, enable_griddata=False):\n",
    "    for which_gfs in [\"y_n_gfs\", \"non_y_n_gfs\"]:\n",
    "        add_to_Cfunction_dict_MoL_malloc(MoL_method, which_gfs)\n",
    "        add_to_Cfunction_dict_MoL_free_memory(MoL_method, which_gfs)\n",
    "    add_to_Cfunction_dict_MoL_step_forward_in_time(MoL_method, RHS_string, post_RHS_string, post_post_RHS_string,\n",
    "                                                   enable_rfm=enable_rfm, enable_curviBCs=enable_curviBCs,\n",
    "                                                   enable_SIMD=enable_SIMD, enable_griddata=enable_griddata)\n",
    "    NRPy_basic_defines_MoL_timestepping_struct(MoL_method=MoL_method)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 4: Code Validation against `MoLtimestepping.MoL_new_way` NRPy+ module [Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "As a code validation check, we verify agreement in the dictionary of Butcher tables between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [MoLtimestepping.MoL_new_way](../edit/MoLtimestepping/MoL_new_way.py) module.\n",
    "\n",
    "We generate the header files for each RK method and check for agreement with the NRPy+ module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:01.365919Z",
     "iopub.status.busy": "2021-03-07T17:15:01.365081Z",
     "iopub.status.idle": "2021-03-07T17:15:01.504197Z",
     "shell.execute_reply": "2021-03-07T17:15:01.503581Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      " ### BEGIN VALIDATION TESTS ###\n",
      "VALIDATION TEST PASSED on all files from Euler method\n",
      "VALIDATION TEST PASSED on all files from RK2 Heun method\n",
      "VALIDATION TEST PASSED on all files from RK2 MP method\n",
      "VALIDATION TEST PASSED on all files from RK2 Ralston method\n",
      "VALIDATION TEST PASSED on all files from RK3 method\n",
      "VALIDATION TEST PASSED on all files from RK3 Heun method\n",
      "VALIDATION TEST PASSED on all files from RK3 Ralston method\n",
      "VALIDATION TEST PASSED on all files from SSPRK3 method\n",
      "VALIDATION TEST PASSED on all files from RK4 method\n",
      "VALIDATION TEST PASSED on all files from DP5 method\n",
      "VALIDATION TEST PASSED on all files from DP5alt method\n",
      "VALIDATION TEST PASSED on all files from CK5 method\n",
      "VALIDATION TEST PASSED on all files from DP6 method\n",
      "VALIDATION TEST PASSED on all files from L6 method\n",
      "VALIDATION TEST PASSED on all files from DP8 method\n",
      "### END VALIDATION TESTS ###\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import MoLtimestepping.MoL_new_way as MoLC\n",
    "import difflib\n",
    "import pprint\n",
    "# Courtesy https://stackoverflow.com/questions/12956957/print-diff-of-python-dictionaries ,\n",
    "#   which itself is an adaptation of some Cpython core code\n",
    "def compare_dicts(d1, d2):\n",
    "    return ('\\n' + '\\n'.join(difflib.ndiff(\n",
    "                   pprint.pformat(d1).splitlines(),\n",
    "                   pprint.pformat(d2).splitlines())))\n",
    "\n",
    "\n",
    "print(\"\\n\\n ### BEGIN VALIDATION TESTS ###\")\n",
    "import filecmp\n",
    "for key, value in Butcher_dict.items():\n",
    "    register_C_functions_and_NRPy_basic_defines(key,\n",
    "                            \"rhs_eval(Nxx,Nxx_plus_2NGHOSTS,dxx, RK_INPUT_GFS, RK_OUTPUT_GFS);\",\n",
    "                           \"apply_bcs(Nxx,Nxx_plus_2NGHOSTS, RK_OUTPUT_GFS);\")\n",
    "    from outputC import outC_function_dict, outC_function_master_list\n",
    "    notebook_dict = outC_function_dict.copy()\n",
    "    notebook_master_list = list(outC_function_master_list)\n",
    "    outC_function_dict.clear()\n",
    "    del outC_function_master_list[:]\n",
    "    from outputC import outC_function_dict\n",
    "    if outC_function_dict != {}:\n",
    "        print(\"Error in clearing outC_function_dict.\")\n",
    "        sys.exit(1)\n",
    "    MoLC.register_C_functions_and_NRPy_basic_defines(key,\n",
    "                            \"rhs_eval(Nxx,Nxx_plus_2NGHOSTS,dxx, RK_INPUT_GFS, RK_OUTPUT_GFS);\",\n",
    "                           \"apply_bcs(Nxx,Nxx_plus_2NGHOSTS, RK_OUTPUT_GFS);\")\n",
    "    from outputC import outC_function_dict\n",
    "    python_module_dict = outC_function_dict\n",
    "\n",
    "    if notebook_dict != python_module_dict:\n",
    "        print(\"VALIDATION TEST FAILED.\\n\")\n",
    "        print(compare_dicts(notebook_dict, python_module_dict))\n",
    "        sys.exit(1)\n",
    "    print(\"VALIDATION TEST PASSED on all files from \"+str(key)+\" method\")\n",
    "print(\"### END VALIDATION TESTS ###\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: Output this notebook to $\\LaTeX$-formatted PDF \\[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-RK_Butcher_Table_Generating_C_Code.pdf](Tutorial-RK_Butcher_Table_Generating_C_Code.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": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:15:01.508961Z",
     "iopub.status.busy": "2021-03-07T17:15:01.508216Z",
     "iopub.status.idle": "2021-03-07T17:15:04.996994Z",
     "shell.execute_reply": "2021-03-07T17:15:04.996025Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Method_of_Lines-C_Code_Generation_new_way.tex, and\n",
      "    compiled LaTeX file to PDF file Tutorial-Method_of_Lines-\n",
      "    C_Code_Generation_new_way.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Method_of_Lines-C_Code_Generation_new_way\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}