Deprecated/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 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 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",
    "This module describes how three C code blocks are written to implement Method of Lines timestepping for a specified RK method. The first block is dedicated to allocating memory for the appropriate number of grid function lists needed for the given RK method. The second block will implement the Runge Kutta numerical scheme based on the corresponding Butcher table. The third block will free up the previously allocated memory after the Method of Lines run is complete. These blocks of code are stored within the following three header files respectively\n",
    "\n",
    "1. `MoLtimestepping/RK_Allocate_Memory.h`\n",
    "1. `MoLtimestepping/RK_MoL.h`\n",
    "1. `MoLtimestepping/RK_Free_Memory.h`\n",
    "\n",
    "The generated code is then included in future Start-to-Finish example tutorial notebooks when solving PDEs numerically."
   ]
  },
  {
   "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](#allocate): Allocating Memory, `MoLtimestepping/RK_Allocate_Memory.h`\n",
    "    1. [Step 3.b](#rkmol): Implementing the Runge Kutta Scheme for Method of Lines Timestepping,  `MoLtimestepping/RK_MoL.h`\n",
    "    1. [Step 3.c](#free): Freeing Allocated Memory, `MoLtimestepping/RK_Free_Memory.h`\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": "2022-08-05T18:54:11.342313Z",
     "iopub.status.busy": "2022-08-05T18:54:11.341797Z",
     "iopub.status.idle": "2022-08-05T18:54:11.874029Z",
     "shell.execute_reply": "2022-08-05T18:54:11.873236Z"
    }
   },
   "outputs": [],
   "source": [
    "import sympy as sp  # Import SymPy, a computer algebra system written entirely in Python\n",
    "import os, shutil   # Standard Python modules for multiplatform OS-level functions\n",
    "from MoLtimestepping.RK_Butcher_Table_Dictionary import Butcher_dict"
   ]
  },
  {
   "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": "2022-08-05T18:54:11.878844Z",
     "iopub.status.busy": "2022-08-05T18:54:11.878616Z",
     "iopub.status.idle": "2022-08-05T18:54:11.886330Z",
     "shell.execute_reply": "2022-08-05T18:54:11.885112Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The RK method Euler is diagonal!\n",
      "The RK method RK2 Heun is diagonal!\n",
      "The RK method RK2 MP is diagonal!\n",
      "The RK method RK2 Ralston is diagonal!\n",
      "The RK method RK3 is NOT diagonal!\n",
      "The RK method RK3 Heun is diagonal!\n",
      "The RK method RK3 Ralston is diagonal!\n",
      "The RK method SSPRK3 is NOT diagonal!\n",
      "The RK method RK4 is diagonal!\n",
      "The RK method DP5 is NOT diagonal!\n",
      "The RK method DP5alt is NOT diagonal!\n",
      "The RK method CK5 is NOT diagonal!\n",
      "The RK method DP6 is NOT diagonal!\n",
      "The RK method L6 is NOT diagonal!\n",
      "The RK method DP8 is NOT diagonal!\n"
     ]
    }
   ],
   "source": [
    "def diagonal(key):\n",
    "    diagonal = True #  Start with the Butcher table is diagonal\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 element is non-zero, then the table is not diagonal\n",
    "                diagonal = False\n",
    "                break\n",
    "        row_idx += 1 # Update to check the next row\n",
    "    return diagonal\n",
    "\n",
    "# State whether each Butcher table is diagonal or not\n",
    "for key, value in Butcher_dict.items():\n",
    "    if diagonal(key) == True:\n",
    "        print(\"The RK method \"+str(key)+\" is diagonal!\")\n",
    "    else:\n",
    "        print(\"The RK 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 for solving PDEs. To see what the C code looks like for a particular method, simply change the `RK_method` below, otherwise it will default to `\"RK4\"`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-08-05T18:54:11.938672Z",
     "iopub.status.busy": "2022-08-05T18:54:11.937470Z",
     "iopub.status.idle": "2022-08-05T18:54:11.941842Z",
     "shell.execute_reply": "2022-08-05T18:54:11.941020Z"
    }
   },
   "outputs": [],
   "source": [
    "# Choose a method to see the C code print out for\n",
    "RK_method = \"RK3 Ralston\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='free'></a>\n",
    "\n",
    "## Step 3.a: Freeing Allocated Memory, `MoLtimestepping/RK_Free_Memory.h` [Back to [top](#toc)\\]\n",
    "$$\\label{free}$$\n",
    "\n",
    "We define the function `RK_free()` which generates the C code for freeing the memory that was being occupied by the grid functions lists that had been allocated. The function writes the C code to the header file `MoLtimestepping/RK_Free_Memory.h`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-08-05T18:54:11.946251Z",
     "iopub.status.busy": "2022-08-05T18:54:11.945977Z",
     "iopub.status.idle": "2022-08-05T18:54:11.951576Z",
     "shell.execute_reply": "2022-08-05T18:54:11.950621Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.a: When allocating memory, we populate a list malloced_gridfunctions,\n",
    "#           which is used here to determine which gridfunctions need memory freed,\n",
    "#           via the free() command. Free the mallocs!\n",
    "def free_allocated_memory(outdir,RK_method,malloced_gridfunctions):\n",
    "    # This step is made extremely easy, as we had to\n",
    "    with open(os.path.join(outdir, \"RK_Free_Memory.h\"), \"w\") as file:\n",
    "        file.write(\"// Code snippet freeing gridfunction memory for \\\"\" + RK_method + \"\\\" method:\\n\")\n",
    "\n",
    "        for gridfunction in malloced_gridfunctions:\n",
    "            file.write(\"free(\" + gridfunction + \");\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rkmol'></a>\n",
    "\n",
    "## Step 3.b: Implementing the Runge Kutta Scheme for Method of Lines Timestepping,  `MoLtimestepping/RK_MoL.h` [Back to [top](#toc)\\]\n",
    "$$\\label{rkmol}$$\n",
    "\n",
    "We define the function `RK_MoL()` which generates the C code for implementing Method of Lines using a specified Runge Kutta scheme. The function writes the C code to the header file `MoLtimestepping/RK_MoL.h`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-08-05T18:54:11.956004Z",
     "iopub.status.busy": "2022-08-05T18:54:11.955712Z",
     "iopub.status.idle": "2022-08-05T18:54:12.037408Z",
     "shell.execute_reply": "2022-08-05T18:54:12.036625Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This is the MoL timestepping RK scheme C code for the RK3 Ralston method: \n",
      "\n",
      "// C code implementation of RK3 Ralston Method of Lines timestepping.\n",
      "\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",
      "// ***k1 substep:***\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_eval(Nxx,Nxx_plus_2NGHOSTS,dxx, y_n_gfs, k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs);\n",
      "LOOP_ALL_GFS_GPS(i) {\n",
      "  k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs[i] = (1.0/2.0)*k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs[i]*dt + y_n_gfs[i];\n",
      "}\n",
      "apply_bcs(Nxx,Nxx_plus_2NGHOSTS, k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs);\n",
      "\n",
      "\n",
      "// ***k2 substep:***\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",
      "\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_eval(Nxx,Nxx_plus_2NGHOSTS,dxx, k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs, k2_or_y_nplus_a32_k2_gfs);\n",
      "LOOP_ALL_GFS_GPS(i) {\n",
      "  k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs[i] = (2.0/9.0)*(k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs[i] - y_n_gfs[i])/(1.0/2.0) + y_n_gfs[i] + (1.0/3.0)*k2_or_y_nplus_a32_k2_gfs[i]*dt;\n",
      "  k2_or_y_nplus_a32_k2_gfs[i] = (3.0/4.0)*k2_or_y_nplus_a32_k2_gfs[i]*dt + y_n_gfs[i];\n",
      "}\n",
      "apply_bcs(Nxx,Nxx_plus_2NGHOSTS, k2_or_y_nplus_a32_k2_gfs);\n",
      "apply_bcs(Nxx,Nxx_plus_2NGHOSTS, k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs);\n",
      "\n",
      "\n",
      "// ***k3 substep:***\n",
      "//    1. Add k3 to the running total and save to y_n\n",
      "\n",
      "// Post-RHS evaluation:\n",
      "//    1. Apply post-RHS to y_n\n",
      "rhs_eval(Nxx,Nxx_plus_2NGHOSTS,dxx, k2_or_y_nplus_a32_k2_gfs, y_n_gfs);\n",
      "LOOP_ALL_GFS_GPS(i) {\n",
      "  y_n_gfs[i] = k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs[i] + (4.0/9.0)*y_n_gfs[i]*dt;\n",
      "}\n",
      "apply_bcs(Nxx,Nxx_plus_2NGHOSTS, y_n_gfs);\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Step 3.b: Main driver function for outputting all the MoL C Code\n",
    "def MoL_C_Code_Generation(RK_method = \"RK4\", RHS_string = \"\", post_RHS_string = \"\",outdir=\"MoLtimestepping/\",\n",
    "                          MemAllocOnly=False):\n",
    "\n",
    "\n",
    "####### Step 3.b.i: Allocating Memory\n",
    "    malloc_str = \"// Code snippet allocating gridfunction memory for \\\"\" + RK_method + \"\\\" method:\\n\"\n",
    "\n",
    "    # Loop over grids\n",
    "    malloced_gridfunctions = []\n",
    "    # Set gridfunction type\n",
    "    type_str = \"REAL *restrict \"\n",
    "    # Define a couple useful functions for outputting the needed C code for allocating memory\n",
    "    def malloc_gfs_str(varname):\n",
    "        malloced_gridfunctions.append(varname)\n",
    "        memory_alloc_str = \" = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS_tot\"+\")\"\n",
    "        return type_str + varname + memory_alloc_str + \";\\n\"\n",
    "    def diagnostic_output_gfs_equal_to(gfs):\n",
    "        return type_str + \"diagnostic_output_gfs\"+\" = \"+gfs + \";\\n\"\n",
    "\n",
    "    # No matter the method we define gridfunctions \"y_n_gfs\" to store the initial data\n",
    "    malloc_str += malloc_gfs_str(\"y_n_gfs\")\n",
    "    if diagonal(RK_method) == True and \"RK3\" in RK_method:\n",
    "        malloc_str += malloc_gfs_str(\"k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs\")\n",
    "        malloc_str += malloc_gfs_str(\"k2_or_y_nplus_a32_k2_gfs\")\n",
    "        malloc_str += diagnostic_output_gfs_equal_to(\"k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs\")\n",
    "    else:\n",
    "        if diagonal(RK_method) == False: # 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[RK_method][0])-1\n",
    "            # For non-diagonal tables an intermediate gridfunction \"next_y_input\" is used for rhs evaluations\n",
    "            malloc_str += malloc_gfs_str(\"next_y_input_gfs\")\n",
    "            for i in range(num_k): # Need to allocate all k_i steps for a given method\n",
    "                malloc_str += malloc_gfs_str(\"k\"+str(i+1)+\"_gfs\")\n",
    "            malloc_str += diagnostic_output_gfs_equal_to(\"k1_gfs\")\n",
    "        else: # Allocate memory for diagonal Butcher tables, which use a \"y_nplus1_running_total gridfunction\"\n",
    "            malloc_str += malloc_gfs_str(\"y_nplus1_running_total_gfs\")\n",
    "            if RK_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",
    "                malloc_str += malloc_gfs_str(\"k_odd_gfs\")\n",
    "                malloc_str += malloc_gfs_str(\"k_even_gfs\")\n",
    "            malloc_str += diagnostic_output_gfs_equal_to(\"y_nplus1_running_total_gfs\")\n",
    "    with open(os.path.join(outdir,\"RK_Allocate_Memory.h\"), \"w\") as file:\n",
    "        file.write(malloc_str)\n",
    "\n",
    "    if MemAllocOnly:\n",
    "        free_allocated_memory(outdir,RK_method,malloced_gridfunctions)\n",
    "        return\n",
    "########################################################################################################################\n",
    "# EXAMPLE\n",
    "# ODE: y' = f(t,y), y(t_0) = y_0\n",
    "# Starting at time t_n with solution having value y_n and trying to update to y_nplus1 with timestep dt\n",
    "\n",
    "# Example of scheme for RK4 with k_1, k_2, k_3, k_4 (Using non-diagonal algorithm) Notice this requires storage of\n",
    "# y_n, y_nplus1, k_1 through k_4\n",
    "\n",
    "# k_1      = dt*f(t_n, y_n)\n",
    "# k_2      = dt*f(t_n + 1/2*dt, y_n + 1/2*k_1)\n",
    "# k_3      = dt*f(t_n + 1/2*dt, y_n + 1/2*k_2)\n",
    "# k_4      = dt*f(t_n + dt, y_n + k_3)\n",
    "# y_nplus1 = y_n + 1/3k_1 + 1/6k_2 + 1/6k_3 + 1/3k_4\n",
    "\n",
    "# Example of scheme RK4 using only k_odd and k_even (Diagonal algorithm) Notice that this only requires storage\n",
    "\n",
    "# k_odd     = dt*f(t_n, y_n)\n",
    "# y_nplus1  = 1/3*k_odd\n",
    "# k_even    = dt*f(t_n + 1/2*dt, y_n + 1/2*k_odd)\n",
    "# y_nplus1 += 1/6*k_even\n",
    "# k_odd     = dt*f(t_n + 1/2*dt, y_n + 1/2*k_even)\n",
    "# y_nplus1 += 1/6*k_odd\n",
    "# k_even    = dt*f(t_n + dt, y_n + k_odd)\n",
    "# y_nplus1 += 1/3*k_even\n",
    "########################################################################################################################\n",
    "\n",
    "####### Step 3.b.ii: Implementing the Runge Kutta Scheme for Method of Lines Timestepping\n",
    "    Butcher = Butcher_dict[RK_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",
    "    # Diagonal RK3 only!!!\n",
    "\n",
    "    def single_RK_substep(commentblock, RHS_str, RHS_input_str, RHS_output_str, RK_lhss_list, RK_rhss_list,\n",
    "                          post_RHS_list, post_RHS_output_list, indent = \"  \"):\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",
    "        # Part 1: RHS evaluation:\n",
    "        return_str += RHS_str.replace(\"RK_INPUT_GFS\", RHS_input_str).\\\n",
    "                              replace(\"RK_OUTPUT_GFS\",RHS_output_str)+\"\\n\"\n",
    "\n",
    "        # Part 2: RK update\n",
    "        return_str += \"LOOP_ALL_GFS_GPS\"+\"(i) {\\n\"\n",
    "        for lhs,rhs in zip(RK_lhss_list,RK_rhss_list):\n",
    "            return_str += indent + lhs + \"[i] = \" + rhs.replace(\"_gfs\",\"_gfs\") + \";\\n\"\n",
    "        return_str += \"}\\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 += post_RHS.replace(\"RK_OUTPUT_GFS\",post_RHS_output)+\"\\n\"\n",
    "\n",
    "        return return_str+\"\\n\"\n",
    "\n",
    "    RK_str = \"// C code implementation of \" + RK_method + \" Method of Lines timestepping.\\n\"\n",
    "\n",
    "    if diagonal(RK_method) == True and \"RK3\" in RK_method:\n",
    "        #  In a diagonal RK3 method, only 3 gridfunctions need be defined. Below implements this approach.\n",
    "\n",
    "        # k_1\n",
    "        RK_str += \"\"\"\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",
    "\n",
    "        RK_str += single_RK_substep(\n",
    "            commentblock = \"\"\"\n",
    "// ***k1 substep:***\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\", 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 = [\"(\"+sp.ccode(Butcher[1][1]).replace(\"L\",\"\")+\")*k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs[i]*dt + y_n_gfs[i]\"],\n",
    "            post_RHS_list = [post_RHS_string], post_RHS_output_list = [\"k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs\"])\n",
    "\n",
    "        # k_2\n",
    "        RK_str += single_RK_substep(\n",
    "            commentblock=\"\"\"\n",
    "// ***k2 substep:***\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",
    "\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\", 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=[\"(\"+sp.ccode(Butcher[3][1]).replace(\"L\",\"\")+\")*(k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs[i] - y_n_gfs[i])/(\"+sp.ccode(Butcher[1][1]).replace(\"L\",\"\")+\") + y_n_gfs[i] + (\"+sp.ccode(Butcher[3][2]).replace(\"L\",\"\")+\")*k2_or_y_nplus_a32_k2_gfs[i]*dt\",\n",
    "                          \"(\"+sp.ccode(Butcher[2][2]).replace(\"L\",\"\")+\")*k2_or_y_nplus_a32_k2_gfs[i]*dt + y_n_gfs[i]\"],\n",
    "            post_RHS_list=[post_RHS_string,post_RHS_string],\n",
    "            post_RHS_output_list=[\"k2_or_y_nplus_a32_k2_gfs\",\"k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs\"])\n",
    "\n",
    "        # k_3\n",
    "        RK_str += single_RK_substep(\n",
    "            commentblock=\"\"\"\n",
    "// ***k3 substep:***\n",
    "//    1. Add k3 to the running total and save to y_n\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\",\"k2_or_y_nplus_a32_k2_gfs\"],\n",
    "            RK_rhss_list=[\"k1_or_y_nplus_a21_k1_or_y_nplus1_running_total_gfs[i] + (\"+sp.ccode(Butcher[3][3]).replace(\"L\",\"\")+\")*y_n_gfs[i]*dt\"],\n",
    "            post_RHS_list=[post_RHS_string],\n",
    "            post_RHS_output_list=[\"y_n_gfs\"])\n",
    "    else:\n",
    "        y_n           = \"y_n_gfs\"\n",
    "        if diagonal(RK_method) == False:\n",
    "            for s in range(num_steps):\n",
    "                next_y_input  = \"next_y_input_gfs\"\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 = \"k\" + str(s + 1) + \"_gfs\"\n",
    "                if s == num_steps-1: # If on final step:\n",
    "                    RK_lhs = y_n\n",
    "                    RK_rhs = y_n + \"[i] + dt*(\"\n",
    "                else:                # If on anything but the final step:\n",
    "                    RK_lhs = next_y_input\n",
    "                    RK_rhs = y_n + \"[i] + dt*(\"\n",
    "                for m in range(s+1):\n",
    "                    if Butcher[s+1][m+1] != 0:\n",
    "                        if Butcher[s+1][m+1] != 1:\n",
    "                            RK_rhs += \" + k\"+str(m+1)+\"_gfs[i]*(\"+sp.ccode(Butcher[s+1][m+1]).replace(\"L\",\"\")+\")\"\n",
    "                        else:\n",
    "                            RK_rhs += \" + k\"+str(m+1)+\"_gfs[i]\"\n",
    "                RK_rhs += \" )\"\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",
    "                RK_str += single_RK_substep(\n",
    "                    commentblock=\"// ***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",
    "        else:\n",
    "            y_nplus1_running_total = \"y_nplus1_running_total_gfs\"\n",
    "            if RK_method == 'Euler': # Euler's method doesn't require any k_i, and gets its own unique algorithm\n",
    "                RK_str += single_RK_substep(\n",
    "                    commentblock=\"// ***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+\"[i] + \"+y_nplus1_running_total+\"[i]*dt\"],\n",
    "                    post_RHS_list=[post_RHS_string],\n",
    "                    post_RHS_output_list=[y_n])\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  = \"y_n_gfs\"\n",
    "                        RHS_output = \"k_odd_gfs\"\n",
    "                    # For the remaining steps the inputs and outputs alternate between k_odd and k_even\n",
    "                    elif s%2 == 0:\n",
    "                        RHS_input = \"k_even_gfs\"\n",
    "                        RHS_output = \"k_odd_gfs\"\n",
    "                    else:\n",
    "                        RHS_input = \"k_odd_gfs\"\n",
    "                        RHS_output = \"k_even_gfs\"\n",
    "\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+\"[i]*dt*(\"+sp.ccode(Butcher[num_steps][s+1]).replace(\"L\",\"\")+\")\")\n",
    "\n",
    "                            RK_lhs_list.append(RHS_output)\n",
    "                            RK_rhs_list.append(y_n+\"[i] + \"+RHS_output+\"[i]*dt*(\"+sp.ccode(Butcher[s+1][s+1]).replace(\"L\",\"\")+\")\")\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+\"[i] + \"+RHS_output+\"[i]*dt*(\"+sp.ccode(Butcher[num_steps][s+1]).replace(\"L\",\"\")+\")\")\n",
    "                                else:\n",
    "                                    RK_rhs_list.append(y_nplus1_running_total+\"[i] + \"+RHS_output+\"[i]*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+\"[i] + \"+RHS_output+\"[i]*dt*(\"+sp.ccode(Butcher[s+1][s+1]).replace(\"L\",\"\")+\")\")\n",
    "                                else:\n",
    "                                    RK_rhs_list.append(y_n+\"[i] + \"+RHS_output+\"[i]*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+\"[i] + \"+y_nplus1_running_total+\"[i] + \"+RHS_output+\"[i]*dt*(\"+sp.ccode(Butcher[num_steps][s+1]).replace(\"L\",\"\")+\")\")\n",
    "                            else:\n",
    "                                RK_rhs_list.append(y_n+\"[i] + \"+y_nplus1_running_total+\"[i] + \"+RHS_output+\"[i]*dt)\")\n",
    "                        post_RHS_output = y_n\n",
    "\n",
    "                    RK_str += single_RK_substep(\n",
    "                        commentblock=\"// ***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",
    "\n",
    "    with open(os.path.join(outdir,\"RK_MoL.h\"), \"w\") as file:\n",
    "        file.write(RK_str)\n",
    "\n",
    "####### Step 3.b.iii: Freeing Allocated Memory\n",
    "    free_allocated_memory(outdir,RK_method,malloced_gridfunctions)\n",
    "\n",
    "MoL_C_Code_Generation(RK_method,\"rhs_eval(Nxx,Nxx_plus_2NGHOSTS,dxx, RK_INPUT_GFS, RK_OUTPUT_GFS);\",\n",
    "                                \"apply_bcs(Nxx,Nxx_plus_2NGHOSTS, RK_OUTPUT_GFS);\")\n",
    "print(\"This is the MoL timestepping RK scheme C code for the \"+str(RK_method)+\" method: \\n\")\n",
    "with open(os.path.join(\"MoLtimestepping/\",\"RK_MoL.h\"), \"r\") as file:\n",
    "    print(file.read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 4: Code Validation against `MoLtimestepping.RK_Butcher_Table_Generating_C_Code` 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.RK_Butcher_Table_Generating_C_Code](../edit/MoLtimestepping/RK_Butcher_Table_Generating_C_Code.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": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-08-05T18:54:12.040184Z",
     "iopub.status.busy": "2022-08-05T18:54:12.040001Z",
     "iopub.status.idle": "2022-08-05T18:54:12.154115Z",
     "shell.execute_reply": "2022-08-05T18:54:12.153215Z"
    }
   },
   "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.C_Code_Generation as MoLC\n",
    "\n",
    "print(\"\\n\\n ### BEGIN VALIDATION TESTS ###\")\n",
    "import filecmp\n",
    "for key, value in Butcher_dict.items():\n",
    "    MoL_C_Code_Generation(key,\"rhs_eval(Nxx,Nxx_plus_2NGHOSTS,dxx, RK_INPUT_GFS, RK_OUTPUT_GFS);\",\n",
    "                              \"apply_bcs(Nxx,Nxx_plus_2NGHOSTS, RK_OUTPUT_GFS);\")\n",
    "    for filename in [\"RK_Allocate_Memory.h\",\"RK_MoL.h\",\"RK_Free_Memory.h\"]:\n",
    "        shutil.copy(os.path.join(\"MoLtimestepping/\",filename), os.path.join(\"MoLtimestepping/\",filename+key+\".h\"))\n",
    "\n",
    "    MoLC.MoL_C_Code_Generation(key, \"rhs_eval(Nxx,Nxx_plus_2NGHOSTS,dxx, RK_INPUT_GFS, RK_OUTPUT_GFS);\",\n",
    "                                    \"apply_bcs(Nxx,Nxx_plus_2NGHOSTS, RK_OUTPUT_GFS);\")\n",
    "    for filename in [\"RK_Allocate_Memory.h\",\"RK_MoL.h\",\"RK_Free_Memory.h\"]:\n",
    "        if filecmp.cmp(os.path.join(\"MoLtimestepping/\",filename), os.path.join(\"MoLtimestepping/\",filename+key+\".h\")) == False:\n",
    "            print(\"VALIDATION TEST FAILED ON files: \"+os.path.join(\"MoLtimestepping/\",filename)+\" and \"+\n",
    "                  os.path.join(\"MoLtimestepping/\",filename+key+\".h\"))\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": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-08-05T18:54:12.158185Z",
     "iopub.status.busy": "2022-08-05T18:54:12.157905Z",
     "iopub.status.idle": "2022-08-05T18:54:18.204029Z",
     "shell.execute_reply": "2022-08-05T18:54:18.202686Z"
    }
   },
   "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",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}