Tutorial-Method_of_Lines-C_Code_Generation.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 notebook details the generation of C Code for executing the Method of Lines (MoL) timestepping using Explicit Runge Kutta methods to solve partial differential equations (PDEs). The approach encompasses creating three core C functions for memory management and solution evolution of gridfunctions. This contributes to efficient MoL timestepping and finds use in various numerical applications, notably hyperbolic PDEs.\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, the 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 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. and 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 below.\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": [],
   "source": [
    "# Check if Butcher Table is diagonal\n",
    "def diagonal(key):\n",
    "    # Get the Butcher table corresponding to the provided key\n",
    "    Butcher = Butcher_dict[key][0]\n",
    "\n",
    "    # Establish the number of rows to check for diagonal trait, all but the last row\n",
    "    L = len(Butcher) - 1\n",
    "\n",
    "    # Initialize the Butcher table row index\n",
    "    row_idx = 0\n",
    "\n",
    "    # Check all the desired rows\n",
    "    for i in range(L):\n",
    "        # Check each element before the diagonal element in a row\n",
    "        for j in range(1, row_idx):\n",
    "            # If any non-diagonal coefficient is non-zero, then the table is not diagonal\n",
    "            if Butcher[i][j] != sp.sympify(0):\n",
    "                return False\n",
    "\n",
    "        # Update to check the next row\n",
    "        row_idx += 1\n",
    "\n",
    "    # If all non-diagonal elements are zero, the table is diagonal\n",
    "    return True"
   ]
  },
  {
   "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": [
    "def generate_gridfunction_names(MoL_method=\"RK4\"):\n",
    "    \"\"\"\n",
    "    Generate gridfunction names for the specified Method of Lines (MoL) method.\n",
    "\n",
    "    :param MoL_method: The MoL method to generate gridfunction names for.\n",
    "    :return: A tuple containing y_n_gridfunctions, non_y_n_gridfunctions_list,\n",
    "             diagnostic_gridfunctions_point_to, and diagnostic_gridfunctions2_point_to.\n",
    "    \"\"\"\n",
    "    # y_n_gridfunctions store data for the vector of gridfunctions y_i at t_n,\n",
    "    # the start of each MoL timestep.\n",
    "    y_n_gridfunctions = \"y_n_gfs\"\n",
    "\n",
    "    # non_y_n_gridfunctions are needed to compute the data at t_{n+1}.\n",
    "    # Often labeled with k_i in the name, these gridfunctions are *not*\n",
    "    # needed at the start of each timestep.\n",
    "    non_y_n_gridfunctions_list = []\n",
    "\n",
    "    # Diagnostic output gridfunctions diagnostic_output_gfs & diagnostic_output_gfs2.\n",
    "    diagnostic_gridfunctions2_point_to = \"\"\n",
    "\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",
    "        diagnostic_gridfunctions2_point_to = \"k2_or_y_nplus_a32_k2_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{}_gfs\".format(i + 1))\n",
    "            diagnostic_gridfunctions_point_to = \"k1_gfs\"\n",
    "            if \"k2_gfs\" in non_y_n_gridfunctions_list:\n",
    "                diagnostic_gridfunctions2_point_to = \"k2_gfs\"\n",
    "            else:\n",
    "                print(\"MoL WARNING: No gridfunction group available for diagnostic_output_gfs2\")\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",
    "            diagnostic_gridfunctions2_point_to = \"k_odd_gfs\"\n",
    "    non_y_n_gridfunctions_list.append(\"auxevol_gfs\")\n",
    "\n",
    "    return y_n_gridfunctions, non_y_n_gridfunctions_list, \\\n",
    "        diagnostic_gridfunctions_point_to, diagnostic_gridfunctions2_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",
    "    # Create a description for the function\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, diagnostic_gridfunctions2_point_to = \\\n",
    "        generate_gridfunction_names(MoL_method=MoL_method)\n",
    "\n",
    "    # Determine which gridfunctions to allocate memory for\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",
    "\n",
    "    name = \"MoL_malloc_\" + which_gfs\n",
    "    params = \"const paramstruct *restrict params, MoL_gridfunctions_struct *restrict gridfuncs\"\n",
    "\n",
    "    # Generate the body of the function\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",
    "    body += \"\\ngridfuncs->diagnostic_output_gfs2 = gridfuncs->\" + diagnostic_gridfunctions2_point_to + \";\\n\"\n",
    "\n",
    "    # Add the function to the C function dictionary\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(comment_block, substep_time_offset_dt, RHS_str, RHS_input_str, RHS_output_str, RK_lhs_list, RK_rhs_list,\n",
    "                                     post_RHS_list, post_RHS_output_list, enable_SIMD=False,\n",
    "                                     gf_aliases=\"\", post_post_RHS_string=\"\"):\n",
    "    return_str = comment_block + \"\\n\"\n",
    "    substep_time_offset_str = \"{:.17e}\".format(float(substep_time_offset_dt))\n",
    "    return_str += \"griddata->params.time = time_start + \" + substep_time_offset_str + \" * griddata->params.dt;\\n\"\n",
    "\n",
    "    # Ensure all input lists are lists\n",
    "    RK_lhs_list = [RK_lhs_list] if not isinstance(RK_lhs_list, list) else RK_lhs_list\n",
    "    RK_rhs_list = [RK_rhs_list] if not isinstance(RK_rhs_list, list) else RK_rhs_list\n",
    "    post_RHS_list = [post_RHS_list] if not isinstance(post_RHS_list, list) else post_RHS_list\n",
    "    post_RHS_output_list = [post_RHS_output_list] if not isinstance(post_RHS_output_list, list) else post_RHS_output_list\n",
    "\n",
    "    return_str += \"{\\n\" + indent_Ccode(gf_aliases, \"  \")\n",
    "    indent = \"  \"\n",
    "\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",
    "\n",
    "    var_type = \"REAL_SIMD_ARRAY\" if enable_SIMD else \"REAL\"\n",
    "    RK_lhs_str_list = [indent + \"const REAL_SIMD_ARRAY __RHS_exp_\" + str(i) if enable_SIMD else indent + str(el).replace(\"gfsL\", \"gfs[i]\") for i, el in enumerate(RK_lhs_list)]\n",
    "\n",
    "    read_list = [read for el in RK_rhs_list for read in list(sp.ordered(el.free_symbols))]\n",
    "    read_list_unique = superfast_uniq(read_list)\n",
    "\n",
    "    for el in read_list_unique:\n",
    "        if str(el) != \"params->dt\":\n",
    "            if enable_SIMD:\n",
    "                simd_el = str(el).replace(\"gfsL\", \"gfs[i]\")\n",
    "                return_str += \"{}  const {} {} = ReadSIMD(&{});\\n\".format(indent, var_type, str(el), simd_el)\n",
    "            else:\n",
    "                return_str += \"{}  const {} {} = {};\\n\".format(indent, var_type, str(el),\n",
    "                                                               str(el).replace(\"gfsL\", \"gfs[i]\"))\n",
    "\n",
    "    if enable_SIMD:\n",
    "        return_str += \"{}  const REAL_SIMD_ARRAY DT = ConstSIMD(params->dt);\\n\".format(indent)\n",
    "\n",
    "    pre_indent = \"2\"\n",
    "    kernel = outputC(RK_rhs_list, RK_lhs_str_list, filename=\"returnstring\",\n",
    "                     params=\"includebraces=False,preindent=\"+pre_indent+\",outCverbose=False,enable_SIMD=\"+str(enable_SIMD))\n",
    "    if enable_SIMD:\n",
    "        return_str += kernel.replace(\"params->dt\", \"DT\")\n",
    "        for i, el in enumerate(RK_lhs_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",
    "    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",
    "    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",
    "    params = \"griddata_struct *restrict griddata\"\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",
    "    body += \"// First set the initial time:\\n\"\n",
    "    body += \"const REAL time_start = griddata->params.time;\\n\"\n",
    "\n",
    "    y_n_gridfunctions, non_y_n_gridfunctions_list, _throwaway, _throwaway2 = generate_gridfunction_names(MoL_method)\n",
    "\n",
    "    gf_prefix = \"griddata->gridfuncs.\"\n",
    "\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",
    "\n",
    "    gf_aliases += \"paramstruct *restrict params = &griddata->params;\\n\"\n",
    "    if enable_rfm:\n",
    "        gf_aliases += \"const rfm_struct *restrict rfmstruct = &griddata->rfmstruct;\\n\"\n",
    "    else:\n",
    "        gf_aliases += \"REAL *restrict xx[3]; for(int ww=0;ww<3;ww++) xx[ww] = griddata->xx[ww];\\n\"\n",
    "    if enable_curviBCs:\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",
    "    # 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",
    "    dt = sp.Symbol(\"params->dt\", real=True)\n",
    "\n",
    "    if diagonal(MoL_method) and \"RK3\" in MoL_method:\n",
    "        # Diagonal RK3 only!!!\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",
    "            comment_block=\"\"\"// -={ 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",
    "            substep_time_offset_dt=Butcher[0][0],\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_lhs_list=[k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs],\n",
    "            RK_rhs_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, 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",
    "            comment_block=\"\"\"// -={ 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",
    "            substep_time_offset_dt=Butcher[1][0],\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_lhs_list=[k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs, k2_or_y_nplus_a32_k2_gfs],\n",
    "            RK_rhs_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, 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",
    "            comment_block=\"\"\"// -={ 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",
    "            substep_time_offset_dt=Butcher[2][0],\n",
    "            RHS_str=RHS_string,\n",
    "            RHS_input_str=k2_or_y_nplus_a32_k2_gfs, RHS_output_str=y_n_gfs,\n",
    "            RK_lhs_list=[y_n_gfs],\n",
    "            RK_rhs_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, 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",
    "                    comment_block=\"// -={ START k\" + str(s + 1) + \" substep }=-\",\n",
    "                    substep_time_offset_dt=Butcher[s][0],\n",
    "                    RHS_str=RHS_string,\n",
    "                    RHS_input_str=RHS_input, RHS_output_str=RHS_output,\n",
    "                    RK_lhs_list=[RK_lhs], RK_rhs_list=[RK_rhs],\n",
    "                    post_RHS_list=[post_RHS],\n",
    "                    post_RHS_output_list=[post_RHS_output],\n",
    "                    enable_SIMD=enable_SIMD, 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",
    "                    comment_block=indent + \"// ***Euler timestepping only requires one RHS evaluation***\",\n",
    "                    substep_time_offset_dt=Butcher[0][0],\n",
    "                    RHS_str=RHS_string,\n",
    "                    RHS_input_str=y_n, RHS_output_str=y_nplus1_running_total,\n",
    "                    RK_lhs_list=[y_n], RK_rhs_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, 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",
    "                        comment_block=indent + \"// -={ START k\" + str(s + 1) + \" substep }=-\",\n",
    "                        substep_time_offset_dt=Butcher[s][0],\n",
    "                        RHS_str=RHS_string,\n",
    "                        RHS_input_str=RHS_input, RHS_output_str=RHS_output,\n",
    "                        RK_lhs_list=RK_lhs_list, RK_rhs_list=RK_rhs_list,\n",
    "                        post_RHS_list=[post_RHS_string],\n",
    "                        post_RHS_output_list=[post_RHS_output],\n",
    "                        enable_SIMD=enable_SIMD, gf_aliases=gf_aliases,\n",
    "                        post_post_RHS_string=post_post_RHS_string) + \"// -={ END k\" + str(s + 1) + \" substep }=-\\n\\n\"\n",
    "\n",
    "    body += \"\"\"\n",
    "// To minimize roundoff error (from adding dt to params.time lots of times),\n",
    "//   here we set time based on the iteration number.\n",
    "griddata->params.time = (REAL)(griddata->params.n + 1) * griddata->params.dt;\n",
    "\n",
    "// Finally, increment the timestep n:\n",
    "griddata->params.n++;\n",
    "\"\"\"\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",
    "        _diagnostic_gridfunctions2_point_to = generate_gridfunction_names(MoL_method=MoL_method)\n",
    "\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",
    "        _diagnostic_gridfunctions2_point_to = 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 += indent + \"REAL *restrict diagnostic_output_gfs2;\\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):\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)\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` 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](../edit/MoLtimestepping/MoL.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 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",
    "for key, value in Butcher_dict.items():\n",
    "    if key not in {\"AHE\", \"ABS\", \"ARKF\", \"ACK\", \"ADP5\", \"ADP8\", \"AB\"}:\n",
    "        # This validation does not work on anything other than standard RK methods, so they are excluded.\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.tex, and compiled LaTeX\n",
      "    file to PDF file Tutorial-Method_of_Lines-C_Code_Generation.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\")"
   ]
  }
 ],
 "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}