Tutorial-Solving_the_Scalar_Wave_Equation_with_NumPy.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",
    "# NRPy+ Tutorial: Solving the Scalar Wave Equation with `NumPy`\n",
    "\n",
    "## Authors: Zach Etienne, Patrick Nelson, Terrence Pierre Jacques, Thiago Assumpção, Leo Werneck, Brandon Clark\n",
    "\n",
    "(*Note on Authors*: Zach wrote the NRPy+ infrastructure, as well as this notebook and its Python code; Patrick wrote the first version of the NRPy+-based Einstein Toolkit thorns for solving the scalar wave equation in 2018; Terrence rewrote these thorns to latest NRPy+ standards in 2020, along the lines of the [`BaikalETK` thorns](Tutorial-BaikalETK.ipynb); Thiago extended the scalar wave initial data infrastructure and contributed fixes to the original NRPy+ scalar wave notebooks; Leo created the boundary condition animation below; and Brandon established NRPy+ notebook formatting standards.)\n",
    "\n",
    "This notebook was first written as a tutorial to introduce NRPy+ during the 2020 Einstein Toolkit Workshop.\n",
    "\n",
    "## In this tutorial we will construct and validate a code that solves the scalar wave equation $\\partial_t^2 u = c^2 \\nabla^2 u$ using `NumPy`, to both review the basic components of a numerical PDE solver and motivate the use of NRPy+\n",
    "\n",
    "This notebook was written to explicitly outline the basic algorithms for numerically solving [hyperbolic PDEs](https://en.wikipedia.org/wiki/Hyperbolic_partial_differential_equation), **including Einstein's equations of general relativity in e.g., the BSSN formalism**.\n",
    "\n",
    "While the codes here are written by hand, the objective of the notebook is motivate the use of NRPy+, which can be used to generate much of this code either *automatically* or from validated templates, greatly reducing the possibility of human error. \n",
    "\n",
    "**Notebook Status:** <font color='green'><b>Validated</b></font>\n",
    "\n",
    "**Validation Notes:** The code developed in this tutorial notebook has been confirmed to converge to the exact solution at the expected rate, as documented below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "1. [Step 0](#prereqs): Prerequisites & Additional Reading\n",
    "1. [Step 1](#intro): Introduction: The scalar wave equation\n",
    "    1. [Step 1.a](#problem_statement): Mathematical problem statement\n",
    "    1. [Step 1.b](#solution): Chosen solution to the scalar wave equation\n",
    "    1. [Step 1.c](#initial_condition): Initial condition\n",
    "1. [Step 2](#mol): The Method of Lines (MoL)\n",
    "1. [Step 3](#basicalg): `NumPy` Implementation: Basic Algorithm\n",
    "    1. [Step 3.a](#numgrid_freeparams): Set up the numerical grid and free parameters\n",
    "    1. [Step 3.b](#numpy_id): Set up the initial data\n",
    "    1. [Step 3.c](#numpy_gfs): Allocate memory for the gridfunctions storing $u$ and $v$, and define the indexing macro function\n",
    "    1. [Step 3.d](#numpy_rhss): Define the right-hand sides of the PDEs\n",
    "    1. [Step 3.e](#numpy_bcs): Boundary Conditions\n",
    "    1. [Step 3.f](#numpy_mol): The Method of Lines\n",
    "    1. [Step 3.g](#numpy_driver): The main driver function\n",
    "1. [Step 4](#too_slow): Argh, the code is SLOW! Why use NRPy+ instead?\n",
    "1. [Step 5](#error_analysis): Error analysis & code validation: Confirming numerical errors converge to zero at the expected rate\n",
    "1. [Step 6](#student_exercises): Exercises for students\n",
    "1. [Step 7](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='prereqs'></a>\n",
    "\n",
    "# Step 0: Prerequisites & Additional Reading \\[Back to [top](#toc)\\]\n",
    "$$\\label{prereqs}$$\n",
    "\n",
    "\n",
    "This tutorial assumes basic familiarity with computer programming, undergraduate mathematics, and computational physics or numerical analysis.\n",
    "\n",
    "For additional reading, please consult the following links:\n",
    "\n",
    "* [Online Resources for Numerical Analysis & Basic Mathematics](https://etienneresearch.com/MATH521-f2018/notes__additional_reading.html)\n",
    "* [Numerical Recipes in C](http://www.numerical.recipes/) Feel free to use the *free* \"obsolete\" (but not really!) versions. Note that Numerical Recipes was written by numerical relativists!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='intro'></a>\n",
    "\n",
    "# Step 1: Introduction: The scalar wave equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{intro}$$\n",
    "\n",
    "We will write a Python code (based in [NumPy](https://numpy.org/)) to numerically solve the scalar wave equation\n",
    "\n",
    "$$\\partial_t^2 u(t,x,y,z) = c^2 \\nabla^2 u(t,x,y,z),$$\n",
    "\n",
    "given initial conditions $u(t=t_0,x,y,z)$.\n",
    "\n",
    "<a id='problem_statement'></a>\n",
    "\n",
    "## Step 1.a: Mathematical problem statement \\[Back to [top](#toc)\\]\n",
    "$$\\label{problem_statement}$$\n",
    "\n",
    "We will numerically solve the scalar wave equation as an [initial value problem](https://en.wikipedia.org/wiki/Initial_value_problem) in Cartesian coordinates:\n",
    "$$\\partial_t^2 u = c^2 \\nabla^2 u,$$\n",
    "where $u$ (the amplitude of the wave) is a function of time and space: $u = u(t,x,y,...)$ (spatial dimension as-yet unspecified) and $c$ is the wave speed, subject to some initial condition\n",
    "\n",
    "$$u(0,x,y,...) = f(x,y,...)$$\n",
    "\n",
    "and suitable spatial boundary conditions (we'll stick with simple extrapolation boundary conditions at first).\n",
    "\n",
    "As described in the next section, we will find it quite useful to define\n",
    "$$v(t,x,y,...) = \\partial_t u(t,x,y,...).$$\n",
    "\n",
    "In this way, the second-order PDE is reduced to a set of two coupled first-order PDEs in time\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_t u &= v \\\\\n",
    "\\partial_t v &= c^2 \\nabla^2 u.\n",
    "\\end{align}\n",
    "\n",
    "We will use NRPy+ to generate efficient C codes capable of generating both initial data $u(0,x,y,...) = f(x,y,...)$; $v(0,x,y,...)=g(x,y,...)$, as well as finite-difference expressions for the right-hand sides of the above expressions. These expressions are needed within the *Method of Lines* to \"integrate\" the solution forward in time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='solution'></a>\n",
    "\n",
    "## Step 1.b: Chosen solution to the scalar wave equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{solution}$$\n",
    "\n",
    "Here we will implement the spherical Gaussian solution to the scalar wave equation, which consists of ingoing and outgoing wave fronts:\n",
    "\\begin{align}\n",
    "u(r,t) &= u_{\\rm out}(r,t) + u_{\\rm in}(r,t) + u_{\\rm offset},\\ \\ \\text{where}\\\\\n",
    "u_{\\rm out}(r,t) &=\\frac{r-ct}{r} \\exp\\left[\\frac{-(r-ct)^2}{2 \\sigma^2}\\right] \\\\\n",
    "u_{\\rm in}(r,t) &=\\frac{r+ct}{r} \\exp\\left[\\frac{-(r+ct)^2}{2 \\sigma^2}\\right] \\\\\n",
    "\\end{align}\n",
    "where $c$ is the wavespeed, $u_{\\rm offset}$ is a freely specifiable constant offset, $\\sigma$ is the width of the Gaussian (i.e., the \"standard deviation\"). Next we'll demonstrate using SymPy (a computer algebra system written in Python, on which NRPy+ depends) that the above indeed is a solution to the scalar wave equation.\n",
    "\n",
    "In Cartesian coordinates we have\n",
    "$$\n",
    "\\partial_t^2 u(t,x,y,z) = c^2 \\nabla^2 u(t,x,y,z),\n",
    "$$\n",
    "and we know that\n",
    "$$\n",
    "r = \\sqrt{x^2+y^2+z^2},\n",
    "$$\n",
    "which we implement below using [SymPy](https://www.sympy.org/), the Python computer algebra package that NRPy+ uses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:32.820582Z",
     "iopub.status.busy": "2021-03-07T17:24:32.819429Z",
     "iopub.status.idle": "2021-03-07T17:24:33.200196Z",
     "shell.execute_reply": "2021-03-07T17:24:33.199489Z"
    }
   },
   "outputs": [],
   "source": [
    "import sympy as sp  # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "\n",
    "# Declare independent variables t,x,y,z as (real-valued) SymPy symbols\n",
    "t,x,y,z = sp.symbols('t x y z', real=True)\n",
    "\n",
    "# Declare the parameters c, sigma, and u_offset as (real-valued) SymPy symbols as well.\n",
    "#   In NRPy+ we'd declare these as NRPy+ *parameters*\n",
    "c, sigma, u_offset = sp.symbols('c sigma u_offset', real=True)\n",
    "\n",
    "# Then define r:\n",
    "r = sp.sqrt(x**2 + y**2 + z**2)\n",
    "\n",
    "# Next set up the solution u(t,x,y,z):\n",
    "# First the outgoing wave\n",
    "u_out = (r - c*t)/r * sp.exp(-(r - c*t)**2 / (2*sigma**2))\n",
    "# ... and then the ingoing wave\n",
    "u_in  = (r + c*t)/r * sp.exp(-(r + c*t)**2 / (2*sigma**2))\n",
    "\n",
    "u_exact = u_out + u_in + u_offset\n",
    "\n",
    "# ... and then v, which is the time derivative of u:\n",
    "v_exact = sp.diff(u_exact, t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a visualization of this solution over time $u_{\\rm exact}(t)$ in the $x-z$ plane, generated using [matplotlib](https://matplotlib.org/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:33.209972Z",
     "iopub.status.busy": "2021-03-07T17:24:33.209303Z",
     "iopub.status.idle": "2021-03-07T17:24:33.214118Z",
     "shell.execute_reply": "2021-03-07T17:24:33.214587Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"560\"\n",
       "            height=\"315\"\n",
       "            src=\"https://www.youtube.com/embed/TJOo2JIW53g?rel=0&amp;controls=0&amp;showinfo=0\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7f4ca830b9a0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import IFrame\n",
    "# Youtube\n",
    "IFrame('https://www.youtube.com/embed/TJOo2JIW53g?rel=0&amp;controls=0&amp;showinfo=0',560,315)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's confirm the solution solves the PDE:\n",
    "$$\n",
    "\\partial_t^2 u(t,x,y,z) = c^2 \\nabla^2 u(t,x,y,z),\n",
    "$$\n",
    "by confirming that\n",
    "$$\n",
    "\\partial_t^2 u(t,x,y,z) - c^2 \\nabla^2 u(t,x,y,z) = 0,\n",
    "$$\n",
    "using SymPy to compute the second derivatives of $u$ symbolically. To make it easier for SymPy to simplify the resulting expression, we will split up the above equation into:\n",
    "\n",
    "$$\n",
    "\\partial_t^2 (u_{\\rm in} + u_{\\rm out} + u_{\\rm offset}) - c^2 \\nabla^2 (u_{\\rm in} + u_{\\rm out} + u_{\\rm offset}) = 0,\n",
    "$$\n",
    "and confirm that\n",
    "\\begin{align}\n",
    "\\partial_t^2 u_{\\rm in} - c^2 \\nabla^2 u_{\\rm in} &= 0 \\\\\n",
    "\\partial_t^2 u_{\\rm out} - c^2 \\nabla^2 u_{\\rm out} &= 0 \\\\\n",
    "\\partial_t^2 u_{\\rm offset} - c^2 \\nabla^2 u_{\\rm offset} &= 0,\n",
    "\\end{align}\n",
    "which must be the case since the scalar wave equation is a [linear PDE](https://en.wikipedia.org/wiki/Partial_differential_equation), in which each of the waves (ingoing, outgoing, and constant) must satisfy the wave equation separately:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:33.315662Z",
     "iopub.status.busy": "2021-03-07T17:24:33.279928Z",
     "iopub.status.idle": "2021-03-07T17:24:45.982006Z",
     "shell.execute_reply": "2021-03-07T17:24:45.982578Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(rhs - lhs) = 0\n"
     ]
    }
   ],
   "source": [
    "# Finally confirm that the solution indeed solves the PDE,\n",
    "#   by subtracting the left-hand-side from the right-hand-side\n",
    "#   of the equation and simplifying; we should get zero\n",
    "scalarwave_lhs_in  = sp.diff(u_in,  t, 2)\n",
    "scalarwave_lhs_out = sp.diff(u_out, t, 2)\n",
    "scalarwave_lhs_ost = sp.diff(u_offset, t, 2)\n",
    "scalarwave_rhs_in  = c**2 * (sp.diff(u_in,  x, 2) + sp.diff(u_in,  y, 2) + sp.diff(u_in,  z, 2))\n",
    "scalarwave_rhs_out = c**2 * (sp.diff(u_out, x, 2) + sp.diff(u_out, y, 2) + sp.diff(u_out, z, 2))\n",
    "scalarwave_rhs_ost = c**2 * (sp.diff(u_offset, x, 2) + sp.diff(u_offset, y, 2) + sp.diff(u_offset, z, 2))\n",
    "\n",
    "scalarwave_lhs_minus_rhs_in  = sp.simplify(scalarwave_lhs_in  - scalarwave_rhs_in)\n",
    "scalarwave_lhs_minus_rhs_out = sp.simplify(scalarwave_lhs_out - scalarwave_rhs_out)\n",
    "scalarwave_lhs_minus_rhs_ost = sp.simplify(scalarwave_lhs_ost - scalarwave_rhs_ost)\n",
    "\n",
    "print(\"(rhs - lhs) = %s\" % (scalarwave_lhs_minus_rhs_in+scalarwave_lhs_minus_rhs_out+scalarwave_lhs_minus_rhs_ost))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initial_condition'></a>\n",
    "\n",
    "## Step 1.c: Initial Condition \\[Back to [top](#toc)\\]\n",
    "$$\\label{initial_condition}$$\n",
    "\n",
    "We will choose the above solution at $t=0$, $u(t=0,x,y,z)$, as our initial data and adopt the Method of Lines (described next) to advance the solution forward in time (i.e., solve the initial value problem)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mol'></a>\n",
    "\n",
    "# Step 2: The Method of Lines (MoL) \\[Back to [top](#toc)\\]\n",
    "$$\\label{mol}$$\n",
    "\n",
    "Once we have set our initial conditions (usually referred to as our \"initial data\"), 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)**, so long as the initial value problem PDE can be written in the form\n",
    "$$\\partial_t \\vec{f} = \\mathbf{M}\\ \\vec{f},$$\n",
    "where $\\mathbf{M}$ is an $N\\times N$ matrix filled with 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}$. The scalar wave equation as written in the [previous module](Tutorial-ScalarWave.ipynb)\n",
    "\\begin{equation}\n",
    "\\partial_t \n",
    "\\begin{bmatrix}\n",
    "u \\\\\n",
    "v \n",
    "\\end{bmatrix}=\n",
    "\\begin{bmatrix}\n",
    "0 & 1 \\\\\n",
    "c^2 \\nabla^2  & 0 \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "u \\\\\n",
    "v \n",
    "\\end{bmatrix}\n",
    "\\end{equation}\n",
    "satisfies this requirement. \n",
    "\n",
    "Thus we can treat the spatial derivatives $\\nabla^2 u$ of the scalar wave equation  using **standard finite-difference approaches**, and the temporal derivatives $\\partial_t u$ and $\\partial_t v$ using **standard approaches for solving ODEs**. \n",
    "\n",
    "Here we will apply the highly robust [explicit Runge-Kutta fourth-order scheme](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) (RK4), used widely for numerically solving ODEs, to \"march\" (integrate) the solution vector $\\vec{f}$ forward in time from its initial value (\"initial data\")."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how MoL works.\n",
    "\n",
    "The RK4 method is usually presented for solving the ODE\n",
    "$$\n",
    "y'(t) = f(y,t)\n",
    "$$\n",
    "as follows. Given initial data $y(t_0)=y_0$, one can construct the solution at any later time via the algorithm:\n",
    "\\begin{align}\n",
    "k_1 &= f(y_n, t_n), \\\\\n",
    "k_2 &= f(y_n + \\frac{1}{2}\\Delta tk_1, t_n + \\frac{\\Delta t}{2}), \\\\\n",
    "k_3 &= f(y_n + \\frac{1}{2}\\Delta tk_2, t_n + \\frac{\\Delta t}{2}), \\\\\n",
    "k_4 &= f(y_n + \\Delta tk_3, t_n + \\Delta t), \\\\\n",
    "y_{n+1} &= y_n + \\frac{1}{6}\\Delta t(k_1 + 2k_2 + 2k_3 + k_4) + \\mathcal{O}\\big((\\Delta t)^5\\big).\n",
    "\\end{align}\n",
    "\n",
    "Our PDE involves two variables $u$ and $v$, and the algorithm generalizes in exactly the same manner as it would if we were solving a system of coupled ODEs with two variables. Further our PDE does not contain explicit time dependence, which simplifies the algorithm a bit:\n",
    "\\begin{align}\n",
    "k_{1,u} &= f_u(u_n,v_n) = f_u(v_n) = v_n, \\\\\n",
    "k_{1,v} &= f_v(u_n,v_n) = f_v(u_n) = c^2\\nabla^2 u_n, \\\\\n",
    "k_{2,u} &= f_u\\left(v_n + \\frac{1}{2}\\Delta tk_{1,v}\\right) = v_n + \\frac{1}{2}\\Delta tk_{1,v}\\\\\n",
    "k_{2,v} &= f_v\\left(u_n + \\frac{1}{2}\\Delta tk_{1,u}\\right) = c^2\\nabla^2 \\left(u_n + \\frac{1}{2}\\Delta tk_{1,u}\\right), \\\\\n",
    "k_{3,u} &= f_u\\left(v_n + \\frac{1}{2}\\Delta tk_{2,v}\\right) = v_n + \\frac{1}{2}\\Delta tk_{2,v}\\\\\n",
    "k_{3,v} &= f_v\\left(u_n + \\frac{1}{2}\\Delta tk_{2,u}\\right) = c^2\\nabla^2 \\left(u_n + \\frac{1}{2}\\Delta tk_{2,u}\\right), \\\\\n",
    "k_{4,u} &= f_u(v_n + \\Delta tk_{3,v}) = v_n + \\Delta tk_{3,v}\\\\\n",
    "k_{4,v} &= f_v(u_n + \\Delta tk_{3,u}) = c^2\\nabla^2 \\left(u_n + \\Delta tk_{3,u}\\right), \\\\\n",
    "u_{n+1} &= u_n + \\frac{1}{6}\\Delta t(k_{1,u} + 2k_{2,u} + 2k_{3,u} + k_{4,u}) + \\mathcal{O}\\big((\\Delta t)^5\\big)\\\\\n",
    "v_{n+1} &= v_n + \\frac{1}{6}\\Delta t(k_{1,v} + 2k_{2,v} + 2k_{3,v} + k_{4,v}) + \\mathcal{O}\\big((\\Delta t)^5\\big).\n",
    "\\end{align}\n",
    "\n",
    "Thus, given initial data $u_0$ and $v_0$, we can use the above algorithm to advance the solution forward in time by one timestep, to $u_1$ and $v_1$. Recall the $\\nabla^2 u$ terms in the above expressions are computed using finite-difference derivatives. Since finite-difference derivatives require neighboring points be evaluated, we only evaluate the $k_i$'s in the interior of the grid; at each step we apply boundary conditions to fill in the outermost neighboring points (called ghost zones)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='basicalg'></a>\n",
    "\n",
    "# Step 3: `NumPy` Implementation: Basic Algorithm \\[Back to [top](#toc)\\]\n",
    "$$\\label{basicalg}$$\n",
    "\n",
    "We will store the numerical solution $u$ and its time derivative $v$, *at a given instant in time* on a three-dimensional numerical grid. Since these variables are defined at each point on the numerical grid, we call them **gridfunctions**.\n",
    "\n",
    "We refer to the right-hand side of the equation $\\partial_t \\vec{f} = \\mathbf{M}\\ \\vec{f}$ as the RHS. In this case, we refer to the $\\mathbf{M}\\ \\vec{f}$ as the **scalar wave RHSs**.\n",
    "\n",
    "Armed with these definitions, the basic algorithm for solving the scalar wave equation [initial value problem](https://en.wikipedia.org/wiki/Initial_value_problem), based on the Method of Lines (see section above) is outlined below.\n",
    "\n",
    "1. Set up the numerical grid and free parameters\n",
    "1. Allocate memory for gridfunctions, including temporary storage needed for the RK4 time integration.\n",
    "1. Set gridfunction values to initial data.\n",
    "1. Evolve the system forward in time using RK4 time integration. At each RK4 substep, do the following:\n",
    "    1. Evaluate scalar wave RHS expressions.\n",
    "    1. Apply boundary conditions.\n",
    "\n",
    "In the following sections we will implement this algorithm to solve the scalar wave equation in 3D *by hand* using [NumPy](https://numpy.org/), and to motivate the use of NRPy+."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='numgrid_freeparams'></a>\n",
    "\n",
    "## Step 3.a: `NumPy` Implementation: Set up the numerical grid and free parameters \\[Back to [top](#toc)\\]\n",
    "$$\\label{numgrid_freeparams}$$\n",
    "\n",
    "We will solve the scalar wave equation on a uniform Cartesian grid with `Nx` by `Ny` by `Nz` coordinate points in the $x$, $y$, and $z$ directions respectively. Since the grid is uniform, we can describe the $x$ coordinate of any gridpoint with a single integer $i$, and the same holds true for $y$ and $z$, for which we will use integers $j$ and $k$. Thus we will label each gridpoint $(x_i,y_j,z_k)$.\n",
    "\n",
    "Let's choose a \"cell-centered\" grid, which will store the solution at points\n",
    "$$\n",
    "x_i \\in \\{..., -\\frac{3}{2} \\Delta x, -\\frac{1}{2} \\Delta x, +\\frac{1}{2} \\Delta x, +\\frac{3}{2} \\Delta x  ...\\};\n",
    "$$\n",
    "and we will allow for two additional ghost zones on the outer boundary to account for the fourth-order finite differencing we will implement to numerically compute $\\nabla^2 u$. Thus the expression for computing $x_i$ will be\n",
    "\n",
    "$$\n",
    "x_i = x_{\\rm min} + \\left( (i-\\text{NGHOSTS}) + \\frac{1}{2} \\right) \\Delta x,\n",
    "$$\n",
    "where $\\Delta x$ is the spacing between gridpoints, and $x_{\\rm min}$ denotes the minimum grid value in $x$. We will solve this equation on a cube centered on the origin with the `domain_size=10`, where $x_{\\rm min}=$ `-domain_size`, $y_{\\rm min}=$ `-domain_size`, $z_{\\rm min}=$ `-domain_size`, $x_{\\rm max}=$ `+domain_size`, and so forth. We'll also choose `Nx=Ny=Nz`, so that\n",
    "\n",
    "$$\n",
    "\\Delta x = \\Delta y = \\Delta z = \\frac{x_{\\rm max}-x_{\\rm min}}{\\text{Nx}}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:45.990737Z",
     "iopub.status.busy": "2021-03-07T17:24:45.989971Z",
     "iopub.status.idle": "2021-03-07T17:24:46.128129Z",
     "shell.execute_reply": "2021-03-07T17:24:46.127285Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np  # NumPy: A numerical methods Python module\n",
    "domain_size = 1.0\n",
    "\n",
    "Nx = Ny = Nz = 24\n",
    "# We add two ghostzones for the outer boundary; needed because our\n",
    "#   finite-differencing stencils are two gridpoints wide.\n",
    "NGHOSTS = 2\n",
    "\n",
    "xx = np.zeros(Nx+2*NGHOSTS)\n",
    "yy = np.zeros(Ny+2*NGHOSTS)\n",
    "zz = np.zeros(Nz+2*NGHOSTS)\n",
    "\n",
    "xmin = ymin = zmin = -domain_size\n",
    "xmax = ymax = zmax = +domain_size\n",
    "\n",
    "dx = (xmax - xmin) / Nx\n",
    "for i in range(Nx + 2*NGHOSTS):\n",
    "    xx[i] = xmin + (i - NGHOSTS + 0.5) * dx\n",
    "dy = (ymax - ymin) / Ny\n",
    "for j in range(Ny + 2*NGHOSTS):\n",
    "    yy[j] = ymin + (j - NGHOSTS + 0.5) * dy\n",
    "dz = (zmax - zmin) / Nz\n",
    "for k in range(Nz + 2*NGHOSTS):\n",
    "    zz[k] = zmin + (k - NGHOSTS + 0.5) * dz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we set the free parameters for the scalar wave solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.133210Z",
     "iopub.status.busy": "2021-03-07T17:24:46.132204Z",
     "iopub.status.idle": "2021-03-07T17:24:46.134676Z",
     "shell.execute_reply": "2021-03-07T17:24:46.135257Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set free parameters\n",
    "freeparam_c     = 1.0  # wave speed\n",
    "freeparam_sigma = 3.0  # width of Gaussian\n",
    "freeparam_u_offset=1.0  # offset of solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we set the timestep, which is governed by the [CFL condition](https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition), and the final time `t_final`, relative to the chosen start time $t_0$ (usually $t_0=0$), so that the points closest to origin aren't affected by the approximate boundary condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.140996Z",
     "iopub.status.busy": "2021-03-07T17:24:46.139780Z",
     "iopub.status.idle": "2021-03-07T17:24:46.142295Z",
     "shell.execute_reply": "2021-03-07T17:24:46.142883Z"
    }
   },
   "outputs": [],
   "source": [
    "dt = 0.5*min(dx,min(dy,dz))/freeparam_c\n",
    "t_final = domain_size*0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='numpy_id'></a>\n",
    "\n",
    "## Step 3.b: `NumPy` Implementation: Set up the initial data \\[Back to [top](#toc)\\]\n",
    "$$\\label{numpy_id}$$\n",
    "\n",
    "Now we'll set up `exact_solution_all_points(time, u, v)`, which numerically evaluates the solution for $u$ and $v$ at all gridpoints at a given numerical time `time`.\n",
    "\n",
    "Recall the exact solution is given by\n",
    "\\begin{align}\n",
    "u(r,t) &= u_{\\rm out}(r,t) + u_{\\rm in}(r,t) + 1,\\ \\ \\text{where}\\\\\n",
    "u_{\\rm out}(r,t) &=\\frac{r-ct}{r} \\exp\\left[\\frac{-(r-ct)^2}{2 \\sigma^2}\\right] \\\\\n",
    "u_{\\rm in}(r,t) &=\\frac{r+ct}{r} \\exp\\left[\\frac{-(r+ct)^2}{2 \\sigma^2}\\right].\n",
    "\\end{align}\n",
    "\n",
    "*Exercise for students: Prove that at $t=0$, $v=\\partial_t u \\equiv 0$.*\n",
    "\n",
    "The problem is, SymPy expressions need to be converted to NumPy expressions; otherwise using functions like `sp.N()` will be *incredibly slow*. So we attempt to fix this by some simple string manipulations, some for $v$ were done by hand using the below Python function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.160337Z",
     "iopub.status.busy": "2021-03-07T17:24:46.159714Z",
     "iopub.status.idle": "2021-03-07T17:24:46.162615Z",
     "shell.execute_reply": "2021-03-07T17:24:46.163093Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "u_offset + (-freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))*np.exp(-(-freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k) + (freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))*np.exp(-(freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k)\n",
      "c*np.exp(-(freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k) - c*np.exp(-(-freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k) + c*(-freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2*np.exp(-(-freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/(freeparam_sigma**2*np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k)) - c*(freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2*np.exp(-(freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/(freeparam_sigma**2*np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))\n"
     ]
    }
   ],
   "source": [
    "def opt_string_replace(input):\n",
    "    return input.replace(\"sqrt\",\"np.sqrt\").replace(\"exp\",\"np.exp\").\\\n",
    "      replace(\"x**2\",\"x_i*x_i\").replace(\"y**2\",\"y_j*y_j\").replace(\"z**2\",\"z_k*z_k\").\\\n",
    "      replace(\"c*t\", \"freeparam_c*time\").replace(\"sigma\", \"freeparam_sigma\")\n",
    "\n",
    "print(opt_string_replace(str(u_exact)))\n",
    "print(opt_string_replace(str(v_exact)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.178217Z",
     "iopub.status.busy": "2021-03-07T17:24:46.177520Z",
     "iopub.status.idle": "2021-03-07T17:24:46.179643Z",
     "shell.execute_reply": "2021-03-07T17:24:46.180148Z"
    }
   },
   "outputs": [],
   "source": [
    "def exact_solution_single_pt_u(time, x_i,y_j,z_k):\n",
    "                # Kludge: The following expressions were pasted from above:\n",
    "                return (-freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))*np.exp(-(-freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k) + (freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))*np.exp(-(freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k) + freeparam_u_offset\n",
    "def exact_solution_single_pt_v(time, x_i,y_j,z_k):\n",
    "                # Kludge: The following expressions were pasted from above, and edited slightly by hand\n",
    "                #         to convert the symbol c to the numerical value for c, freeparam_c\n",
    "                return freeparam_c*np.exp(-(freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k) - freeparam_c*np.exp(-(-freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k) + freeparam_c*(-freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2*np.exp(-(-freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/(freeparam_sigma**2*np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k)) - freeparam_c*(freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2*np.exp(-(freeparam_c*time + np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))**2/(2*freeparam_sigma**2))/(freeparam_sigma**2*np.sqrt(x_i*x_i + y_j*y_j + z_k*z_k))\n",
    "\n",
    "def exact_solution_all_points(time, u, v):\n",
    "    for k in range(0, Nz+2*NGHOSTS):\n",
    "        z_k = zz[k]\n",
    "        for j in range(0, Ny+2*NGHOSTS):\n",
    "            y_j = yy[j]\n",
    "            for i in range(0, Nx+2*NGHOSTS):\n",
    "                x_i = xx[i]\n",
    "                u[IDX3D(i,j,k)] = exact_solution_single_pt_u(time, x_i,y_j,z_k)\n",
    "                v[IDX3D(i,j,k)] = exact_solution_single_pt_v(time, x_i,y_j,z_k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To store the solution $u$ and $v$ at all gridpoints on our numerical grid cube requires \n",
    "\n",
    "`2*Nx*Ny*Nz*double`\n",
    "\n",
    "bytes of memory, where `double` is the amount of memory storage (in bytes) needed to store one [double-precision number](https://en.wikipedia.org/wiki/Double-precision_floating-point_format) (this is 8, by the way)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='numpy_gfs'></a>\n",
    "\n",
    "## Step 3.c: `NumPy` Implementation: Allocate memory for the gridfunctions storing $u$ and $v$, and define the indexing macro function \\[Back to [top](#toc)\\]\n",
    "$$\\label{numpy_gfs}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.185227Z",
     "iopub.status.busy": "2021-03-07T17:24:46.184536Z",
     "iopub.status.idle": "2021-03-07T17:24:46.186941Z",
     "shell.execute_reply": "2021-03-07T17:24:46.187508Z"
    }
   },
   "outputs": [],
   "source": [
    "# Allocate memory for gridfunctions. We need ghostzones\n",
    "u = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))\n",
    "v = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As is done in the Einstein Toolkit and native NRPy+ codes, instead of declaring multi-dimensional arrays (e.g., a 3D array), we will instead declare $u$ and $v$ as *one-dimensional* arrays `u[ijk]` and `v[ijk]`, each with `(Nx+2*NGHOSTS)*(Ny+2*NGHOSTS)*(Nz+2*NGHOSTS)` gridpoints. To access data an arbitrary point $(x_i,y_j,z_k)$, we need only call a simple function to find the correct index `ijk` given the grid indices `i`, `j`, and `k`, which label the point $(x_i,y_j,z_k)$:\n",
    "\n",
    "$$\n",
    "\\verb|\n",
    "(i,j,k) = i + (Nx+2*NGHOSTS)*j + (Nx+2*NGHOSTS)*(Ny+2*NGHOSTS)*k = i + (Nx+2*NGHOSTS)*(j + (Ny+2*NGHOSTS)*k)|\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.192468Z",
     "iopub.status.busy": "2021-03-07T17:24:46.191706Z",
     "iopub.status.idle": "2021-03-07T17:24:46.193901Z",
     "shell.execute_reply": "2021-03-07T17:24:46.194416Z"
    }
   },
   "outputs": [],
   "source": [
    "# Define the indexing macro function\n",
    "def IDX3D(i,j,k):\n",
    "    return i + (Nx+2*NGHOSTS)*(j + (Ny+2*NGHOSTS)*k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='numpy_rhss'></a>\n",
    "\n",
    "## Step 3.d: `NumPy` Implementation: Define the right-hand sides of the PDEs \\[Back to [top](#toc)\\]\n",
    "$$\\label{numpy_rhss}$$\n",
    "\n",
    "Next we define the right-hand sides of the $u$ and $v$ equations:\n",
    "\\begin{align}\n",
    "\\partial_t u &= v \\\\\n",
    "\\partial_t v &= c^2 \\nabla^2 u.\n",
    "\\end{align}\n",
    "\n",
    "Again we'll approximate the $\\nabla^2 u$ using fourth-order [finite-difference derivatives](https://en.wikipedia.org/wiki/Finite_difference) (also see [the NRPy+ tutorial on how to compute these expressions automatically or by hand using simple matrix methods](Tutorial-Finite_Difference_Derivatives.ipynb)).\n",
    "\n",
    "Here we'll just use the [Wikipedia article on finite-difference coefficients](https://en.wikipedia.org/wiki/Finite_difference_coefficient) to construct the expressions for\n",
    "\n",
    "$$\n",
    "(\\nabla u)_{i,j,k} = (\\partial_x^2 u)_{i,j,k} + (\\partial_y^2 u)_{i,j,k} + (\\partial_z^2 u)_{i,j,k}\n",
    "$$\n",
    "by hand:\n",
    "\n",
    "The fourth-order finite difference stencil for $(\\partial_x^2 u)_{i,j,k}$ is written\n",
    "\\begin{align}\n",
    "(\\partial_x^2 u)_{i,j,k} &= \\left[-\\frac{1}{12} u_{i-2,j,k} + \\frac{4}{3} u_{i-1,j,k} - \\frac{5}{2} u_{i,j,k} + \\frac{4}{3} u_{i+1,j,k} - \\frac{1}{12} u_{i+2,j,k}\\right]\\frac{1}{(\\Delta x)^2} \\\\\n",
    "&= \\left[-\\frac{1}{12} \\left(u_{i-2,j,k} + u_{i+2,j,k}\\right) + \\frac{4}{3} \\left(u_{i-1,j,k}+u_{i+1,j,k}\\right) - \\frac{5}{2} u_{i,j,k}\\right]\\frac{1}{(\\Delta x)^2},\n",
    "\\end{align}\n",
    "and the expressions can be written for $(\\partial_y^2 u)_{i,j,k}$ and $(\\partial_z^2 u)_{i,j,k}$ can be immediately written based on this pattern:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.205759Z",
     "iopub.status.busy": "2021-03-07T17:24:46.205090Z",
     "iopub.status.idle": "2021-03-07T17:24:46.207213Z",
     "shell.execute_reply": "2021-03-07T17:24:46.207716Z"
    }
   },
   "outputs": [],
   "source": [
    "def eval_rhs_all_interior_points(u, v, u_rhs, v_rhs):\n",
    "    # Notice that if we looped from e.g., k=0, then u[IDX3D(i,j,k-2)] would be OUT OF BOUNDS.\n",
    "    for k in range(NGHOSTS, Nz+NGHOSTS):  # Recall the valid range of k is 0 to Nz+2*NGHOSTS, ...\n",
    "        for j in range(NGHOSTS, Ny+NGHOSTS):  # ... similarly for j and i\n",
    "            for i in range(NGHOSTS, Nx+NGHOSTS):\n",
    "                u_rhs[IDX3D(i,j,k)] = v[IDX3D(i,j,k)]\n",
    "                # First the x-component of nabla\n",
    "                v_rhs[IDX3D(i,j,k)] = freeparam_c**2 * (-1./12. * (u[IDX3D(i-2,j,k)] + u[IDX3D(i+2,j,k)])\n",
    "                                                        +4./3.  * (u[IDX3D(i-1,j,k)] + u[IDX3D(i+1,j,k)])\n",
    "                                                        -5./2.  * u[IDX3D(i,j,k)])/(dx*dx)\n",
    "                # Then the y-component of nabla\n",
    "                v_rhs[IDX3D(i,j,k)]+= freeparam_c**2 * (-1./12. * (u[IDX3D(i,j-2,k)] + u[IDX3D(i,j+2,k)])\n",
    "                                                        +4./3.  * (u[IDX3D(i,j-1,k)] + u[IDX3D(i,j+1,k)])\n",
    "                                                        -5./2.  * u[IDX3D(i,j,k)])/(dy*dy)\n",
    "                # and finally the y-component of nabla\n",
    "                v_rhs[IDX3D(i,j,k)]+= freeparam_c**2 * (-1./12. * (u[IDX3D(i,j,k-2)] + u[IDX3D(i,j,k+2)])\n",
    "                                                        +4./3.  * (u[IDX3D(i,j,k-1)] + u[IDX3D(i,j,k+1)])\n",
    "                                                        -5./2.  * u[IDX3D(i,j,k)])/(dz*dz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='numpy_bcs'></a>\n",
    "\n",
    "## Step 3.e: `NumPy` Implementation: Boundary Conditions \\[Back to [top](#toc)\\]\n",
    "$$\\label{numpy_bcs}$$\n",
    "\n",
    "Notice that the above code does not fill the input gridfunctions $u$ and $v$ in the ghostzones, which will be updated at each Runge-Kutta substep (as outlined next). We will need to apply our spatial boundary conditions to fill in these points. For simplicity let's choose quadratic extrapolation boundary conditions.\n",
    "\n",
    "For example, suppose we are on the lower boundary point in $x$: $u_{1,j,k}$. Then this boundary condition will be written as the quadratic [polynomial extrapolation](https://en.wikipedia.org/wiki/Polynomial_interpolation) taking data from the interior:\n",
    "$$\n",
    "u_{1,j,k} = 3 u_{2,j,k} - 3 u_{3,j,k} + u_{4,j,k}.\n",
    "$$\n",
    "\n",
    "Similarly for the upper boundary point in $x$, the condition becomes:\n",
    "$$\n",
    "u_{\\text{Nx}-2,j,k} = 3 u_{\\text{Nx}-3,j,k} - 3 u_{\\text{Nx}-4,j,k} + u_{\\text{Nx}-5,j,k}.\n",
    "$$\n",
    "\n",
    "\n",
    "We'll apply this algorithm from the innermost boundary point outward, using the approach of filling in the (green-colored) ghost zones as illustrated here in 2 dimensions (*courtesy Leo Werneck*). Extension to 3 dimensions is straightforward.\n",
    "\n",
    "<img src=\"CurviBoundaryConditions/bdrycond_general_algorithm.gif\" width=\"450\" align=\"center\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.228736Z",
     "iopub.status.busy": "2021-03-07T17:24:46.227282Z",
     "iopub.status.idle": "2021-03-07T17:24:46.230267Z",
     "shell.execute_reply": "2021-03-07T17:24:46.230821Z"
    }
   },
   "outputs": [],
   "source": [
    "def bc_face_update(gf, i0min,i0max, i1min,i1max, i2min,i2max, FACEX0,FACEX1,FACEX2):\n",
    "    for i2 in range(i2min,i2max):\n",
    "        for i1 in range(i1min,i1max):\n",
    "            for i0 in range(i0min,i0max):\n",
    "                gf[IDX3D(i0,i1,i2)] = (+3.0*gf[IDX3D(i0+1*FACEX0,i1+1*FACEX1,i2+1*FACEX2)]\n",
    "                                       -3.0*gf[IDX3D(i0+2*FACEX0,i1+2*FACEX1,i2+2*FACEX2)]\n",
    "                                       +1.0*gf[IDX3D(i0+3*FACEX0,i1+3*FACEX1,i2+3*FACEX2)])\n",
    "\n",
    "MAXFACE = -1  # Interp stencil reaches in the negative direction on upper (max) face\n",
    "NUL     = +0\n",
    "MINFACE = +1  # Interp stencil reaches in the positive direction on lower (min) face\n",
    "def apply_extrapolation_bcs(u, v):\n",
    "    for gf in [u,v]:\n",
    "        imin = [NGHOSTS, NGHOSTS, NGHOSTS]\n",
    "        imax = [Nx+NGHOSTS, Ny+NGHOSTS, Nz+NGHOSTS]\n",
    "        for which_gz in range(NGHOSTS):\n",
    "            # After updating each face, adjust imin[] and imax[]\n",
    "            #   to reflect the newly-updated face extents.\n",
    "            bc_face_update(gf, imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2], MINFACE,NUL,NUL); imin[0]-=1\n",
    "            bc_face_update(gf, imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2], MAXFACE,NUL,NUL); imax[0]+=1\n",
    "            bc_face_update(gf, imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2], NUL,MINFACE,NUL); imin[1]-=1\n",
    "            bc_face_update(gf, imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2], NUL,MAXFACE,NUL); imax[1]+=1\n",
    "            bc_face_update(gf, imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2], NUL,NUL,MINFACE); imin[2]-=1\n",
    "            bc_face_update(gf, imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1, NUL,NUL,MAXFACE); imax[2]+=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='numpy_mol'></a>\n",
    "\n",
    "## Step 3.f: `NumPy` Implementation: The Method of Lines \\[Back to [top](#toc)\\]\n",
    "$$\\label{numpy_mol}$$\n",
    "\n",
    "Next we'll set up the Method of Lines (MoL) routine for Runge-Kutta fourth order (RK4), which takes the solution at a given iteration in time $n$, and enables us to advance the solution forward to iteration $n+1$, as outlined above:\n",
    "\n",
    "\\begin{align}\n",
    "k_{1,u} &= f_u(u_n,v_n) = f_u(v_n) = v_n, \\\\\n",
    "k_{1,v} &= f_v(u_n,v_n) = f_v(u_n) = c^2\\nabla^2 u_n, \\\\\n",
    "k_{2,u} &= f_u\\left(v_n + \\frac{1}{2}\\Delta tk_{1,v}\\right) = v_n + \\frac{1}{2}\\Delta tk_{1,v}\\\\\n",
    "k_{2,v} &= f_v\\left(u_n + \\frac{1}{2}\\Delta tk_{1,u}\\right) = c^2\\nabla^2 \\left(u_n + \\frac{1}{2}\\Delta tk_{1,u}\\right), \\\\\n",
    "k_{3,u} &= f_u\\left(v_n + \\frac{1}{2}\\Delta tk_{2,v}\\right) = v_n + \\frac{1}{2}\\Delta tk_{2,v}\\\\\n",
    "k_{3,v} &= f_v\\left(u_n + \\frac{1}{2}\\Delta tk_{2,u}\\right) = c^2\\nabla^2 \\left(u_n + \\frac{1}{2}\\Delta tk_{2,u}\\right), \\\\\n",
    "k_{4,u} &= f_u(v_n + \\Delta tk_{3,v}) = v_n + \\Delta tk_{3,v}\\\\\n",
    "k_{4,v} &= f_v(u_n + \\Delta tk_{3,u}) = c^2\\nabla^2 \\left(u_n + \\Delta tk_{3,u}\\right), \\\\\n",
    "u_{n+1} &= u_n + \\frac{1}{6}\\Delta t(k_{1,u} + 2k_{2,u} + 2k_{3,u} + k_{4,u}) + \\mathcal{O}\\big((\\Delta t)^5\\big)\\\\\n",
    "v_{n+1} &= v_n + \\frac{1}{6}\\Delta t(k_{1,v} + 2k_{2,v} + 2k_{3,v} + k_{4,v}) + \\mathcal{O}\\big((\\Delta t)^5\\big).\n",
    "\\end{align}\n",
    "\n",
    "We will store $k_1$ through $k_4$ as additional gridfunctions, one each for $u$ and $v$, and another gridfunction for $u$ and $v$ (`u_tmp` and `v_tmp`, respectively) for the input into $f_u()$ and $f_v()$ functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.239931Z",
     "iopub.status.busy": "2021-03-07T17:24:46.239242Z",
     "iopub.status.idle": "2021-03-07T17:24:46.242883Z",
     "shell.execute_reply": "2021-03-07T17:24:46.243376Z"
    }
   },
   "outputs": [],
   "source": [
    "u_k1 = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))\n",
    "v_k1 = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))\n",
    "u_k2 = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))\n",
    "v_k2 = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))\n",
    "u_k3 = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))\n",
    "v_k3 = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))\n",
    "u_k4 = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))\n",
    "v_k4 = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))\n",
    "u_tmp = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))\n",
    "v_tmp = np.zeros((Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... then implement a single timestep by calling the `eval_rhs_all_interior_points()` function above with appropriate inputs. Recall that the RK4 algorithm is given by\n",
    "\\begin{align}\n",
    "k_1 &= f(y_n, t_n), \\\\\n",
    "k_2 &= f(y_n + \\frac{1}{2}\\Delta tk_1, t_n + \\frac{\\Delta t}{2}), \\\\\n",
    "k_3 &= f(y_n + \\frac{1}{2}\\Delta tk_2, t_n + \\frac{\\Delta t}{2}), \\\\\n",
    "k_4 &= f(y_n + \\Delta tk_3, t_n + \\Delta t), \\\\\n",
    "y_{n+1} &= y_n + \\frac{1}{6}\\Delta t(k_1 + 2k_2 + 2k_3 + k_4) + \\mathcal{O}\\big((\\Delta t)^5\\big).\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.255245Z",
     "iopub.status.busy": "2021-03-07T17:24:46.254483Z",
     "iopub.status.idle": "2021-03-07T17:24:46.256628Z",
     "shell.execute_reply": "2021-03-07T17:24:46.257115Z"
    }
   },
   "outputs": [],
   "source": [
    "def one_RK_step():\n",
    "    # Compute k_1\n",
    "    eval_rhs_all_interior_points(u, v, u_k1, v_k1)\n",
    "    # Compute inputs into k_2\n",
    "    for idx in range(0, (Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS)):\n",
    "        u_tmp[idx] = u[idx] + 0.5*dt*u_k1[idx]\n",
    "        v_tmp[idx] = v[idx] + 0.5*dt*v_k1[idx]\n",
    "    # Apply BCs to u_tmp and v_tmp:\n",
    "    apply_extrapolation_bcs(u_tmp, v_tmp)\n",
    "\n",
    "    # Compute k_2\n",
    "    eval_rhs_all_interior_points(u_tmp, v_tmp, u_k2, v_k2)\n",
    "    # Compute inputs into k_3\n",
    "    for idx in range(0, (Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS)):\n",
    "        u_tmp[idx] = u[idx] + 0.5*dt*u_k2[idx]\n",
    "        v_tmp[idx] = v[idx] + 0.5*dt*v_k2[idx]\n",
    "    # Apply BCs to u_tmp and v_tmp:\n",
    "    apply_extrapolation_bcs(u_tmp, v_tmp)\n",
    "\n",
    "    # Compute k_3\n",
    "    eval_rhs_all_interior_points(u_tmp, v_tmp, u_k3, v_k3)\n",
    "    # Compute inputs into k_4\n",
    "    for idx in range(0, (Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS)):\n",
    "        u_tmp[idx] = u[idx] + dt*u_k3[idx]\n",
    "        v_tmp[idx] = v[idx] + dt*v_k3[idx]\n",
    "    # Apply BCs to u_tmp and v_tmp:\n",
    "    apply_extrapolation_bcs(u_tmp, v_tmp)\n",
    "\n",
    "    # Compute k_4\n",
    "    eval_rhs_all_interior_points(u_tmp, v_tmp, u_k4, v_k4)\n",
    "    # Finally compute y_{n+1}\n",
    "    for idx in range(0, (Nx+2*NGHOSTS) * (Ny+2*NGHOSTS) * (Nz+2*NGHOSTS)):\n",
    "        u[idx] = u[idx] + (1.0/6.0)*dt*(u_k1[idx] + 2*u_k2[idx] + 2*u_k3[idx] + u_k4[idx])\n",
    "        v[idx] = v[idx] + (1.0/6.0)*dt*(v_k1[idx] + 2*v_k2[idx] + 2*v_k3[idx] + v_k4[idx])\n",
    "    # ... and apply BCs to the updated u and v:\n",
    "    apply_extrapolation_bcs(u, v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='numpy_driver'></a>\n",
    "\n",
    "## Step 3.g: `NumPy` Implementation: The main driver function \\[Back to [top](#toc)\\]\n",
    "$$\\label{numpy_driver}$$\n",
    "\n",
    "Finally we'll write the main driver function, which as a diagnostic outputs the relative error between numerical and exact solutions at the closest point to the center of the numerical grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:24:46.331285Z",
     "iopub.status.busy": "2021-03-07T17:24:46.295351Z",
     "iopub.status.idle": "2021-03-07T17:25:00.619384Z",
     "shell.execute_reply": "2021-03-07T17:25:00.619865Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Outputting data at (x,y,z) = (0.04,0.04,0.04)\n",
      "# Integrating forward in time, to time 0.042 . ETA: N/A seconds\n",
      "# Integrating forward in time, to time 0.083 . ETA: 0m7.51s seconds\n",
      "# Integrating forward in time, to time 0.125 . ETA: 0m7.29s seconds\n",
      "# Integrating forward in time, to time 0.167 . ETA: 0m6.27s seconds\n",
      "# Integrating forward in time, to time 0.208 . ETA: 0m5.38s seconds\n",
      "# Integrating forward in time, to time 0.250 . ETA: 0m4.58s seconds\n",
      "# Integrating forward in time, to time 0.292 . ETA: 0m3.93s seconds\n",
      "# Integrating forward in time, to time 0.333 . ETA: 0m3.25s seconds\n",
      "# Integrating forward in time, to time 0.375 . ETA: 0m2.59s seconds\n",
      "# Integrating forward in time, to time 0.417 . ETA: 0m1.94s seconds\n",
      "# Integrating forward in time, to time 0.458 . ETA: 0m1.30s seconds\n",
      "# Integrating forward in time, to time 0.500 . ETA: 0m0.65s seconds\n",
      "# Results, output to file output_experiment_resolution_24_cubed.txt\n",
      "0.00 -16.00 2.999421380013 2.999421380013\n",
      "0.04 -10.70 2.998843001874 2.998843001814\n",
      "0.08 -10.06 2.997108425138 2.997108424880\n",
      "0.12 -9.70 2.994219322036 2.994219321440\n",
      "0.17 -9.45 2.990178477110 2.990178476038\n",
      "0.21 -9.25 2.984989783456 2.984989781772\n",
      "0.25 -9.09 2.978658237475 2.978658235046\n",
      "0.29 -8.95 2.971189932138 2.971189928832\n",
      "0.33 -8.84 2.962592048775 2.962592044465\n",
      "0.38 -8.73 2.952872847425 2.952872841973\n",
      "0.42 -8.64 2.942041655765 2.942041648971\n",
      "0.46 -8.54 2.930108856521 2.930108848127\n",
      "0.50 -8.48 2.917085872939 2.917085863230\n",
      "\n",
      "CPU times: user 8.52 s, sys: 271 µs, total: 8.52 s\n",
      "Wall time: 8.52 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "initial_time = 0.0\n",
    "\n",
    "# First set up the initial data:\n",
    "exact_solution_all_points(initial_time, u, v)\n",
    "\n",
    "# Store the indices at the point closest to the origin\n",
    "i_o = int((Nx+2*NGHOSTS)/2)\n",
    "j_o = int((Ny+2*NGHOSTS)/2)\n",
    "k_o = int((Nz+2*NGHOSTS)/2)\n",
    "print(\"# Outputting data at (x,y,z) = (%.2f,%.2f,%.2f)\" % (xx[i_o],yy[j_o],zz[k_o]))\n",
    "\n",
    "def diagnostics(n):\n",
    "    # Print the time and the value of the solution closest to the origin\n",
    "    curr_time = initial_time + n*dt\n",
    "    num   = u[IDX3D(i_o, j_o, k_o)]\n",
    "    exact = exact_solution_single_pt_u(curr_time, xx[i_o],yy[j_o],zz[k_o])\n",
    "    log10relerror = np.log10(max(1e-16, np.abs((num-exact)/exact)))\n",
    "    return \"%.2f %.2f %.12f %.12f\\n\" % (curr_time, log10relerror, num, exact)\n",
    "\n",
    "# Output diagnostics at the initial time.\n",
    "out_str = diagnostics(0)\n",
    "\n",
    "# Then integrate forward in time:\n",
    "n_final = int(t_final/dt + 0.5)  # add 0.5 to correct for rounding.\n",
    "n_out_every = int(Nx/24.)  # Output every timestep for Nx=24; every other timestep for Nx=48; etc\n",
    "import time  # for live benchmarking & estimates\n",
    "start = time.time()\n",
    "for n in range(0,n_final):\n",
    "    ETA = \"N/A\"\n",
    "    if n > 0:\n",
    "        time_elapsed_in_seconds = time.time() - start\n",
    "        seconds_per_n = time_elapsed_in_seconds/n\n",
    "        time_remaining_m_field = int((n_final - n)*seconds_per_n/60)\n",
    "        time_remaining_s_field = (n_final - n)*seconds_per_n - time_remaining_m_field*60\n",
    "        ETA = str(time_remaining_m_field)+\"m\"+ '%.2f' % time_remaining_s_field + \"s\"\n",
    "    print(\"# Integrating forward in time, to time %.3f . ETA: %s seconds\" % ((n+1)*dt, ETA))\n",
    "    one_RK_step()\n",
    "    # After the RK step we are at iteration n+1\n",
    "    if((n+1) % n_out_every == 0):\n",
    "        out_str += diagnostics(n+1)\n",
    "\n",
    "experiment_filename = \"output_experiment_resolution_\"+str(Nx)+\"_cubed.txt\"\n",
    "print(\"# Results, output to file \" + experiment_filename)\n",
    "print(out_str)\n",
    "with open(experiment_filename, \"w\") as file:\n",
    "    file.write(out_str)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='too_slow'></a>\n",
    "\n",
    "# Step 4: Argh, the code is SLOW! Why use NRPy+ instead? \\[Back to [top](#toc)\\]\n",
    "$$\\label{too_slow}$$\n",
    "\n",
    "By default the above code outputs data on the `Nx=Ny=Nz=24` = $24^3$ numerical grid, and it takes around 16 seconds to complete (on mybinder). \n",
    "\n",
    "If we were to double the resolution to $48^3$ (keeping the domain size fixed), the number of gridpoints we need to update increases by a factor of 8, and the timestep reduces by a factor of 2; hence the total cost is about 16x higher. Thus it will take roughly 16 seconds, times 16, or roughly 4 minutes to complete a $48^3$ resolution simulation on the same CPU. Similarly, it should take a bit over an hour to complete a simulation with $96^3$ resolution!\n",
    "\n",
    "One reason we wish to use `NRPy+` to convert human-friendly Python expressions (written in `SymPy`) to highly optimized C code is the speed. As we'll see, a C code implementing exactly the same algorithm for solving the scalar wave equation can generate results roughly 10,000x faster!\n",
    "\n",
    "Given how slowly the above Python code solves the scalar wave equation, solving Einstein's equations of general relativity in 3 dimensions with a Python code would be a futile effort. However, speed of execution isn't the only reason to use NRPy+. Here are some more reasons:\n",
    "1. NRPy+ contains a rigid syntax for SymPy symbols that \n",
    "    1. enables you to specify derivatives (e.g., $\\partial_j u$= `u_dD[j]`) and output C code at arbitrary finite differencing (FD) order \n",
    "        1. [**tutorial on FD derivatives**](Tutorial-Finite_Difference_Derivatives.ipynb);\n",
    "        1. [**tutorial on computing FD coefficients**](Tutorial-How_NRPy_Computes_Finite_Difference_Coeffs.ipynb); \n",
    "        1. [**sample C code tutorial**](Tutorial-Start_to_Finish-Finite_Difference_Playground.ipynb)\n",
    "    1. allows for tensorial expressions to be input unambiguously in Einstein-like notation (e.g., $\\gamma_{ij}=$ `gammaDD[i][j]`) \n",
    "        1. [**tutorial on indexed (e.g., tensorial) expressions**](Tutorial-Indexed_Expressions.ipynb)\n",
    "1. NRPy+ automatically implements 15 different RK-like timestepping methods for MoL \n",
    "    1. [**tutorial on NRPy+ dictionary of RK methods**](Tutorial-RK_Butcher_Table_Dictionary.ipynb); \n",
    "    1. [**tutorial validating the NRPy+ RK dictionary**](Tutorial-RK_Butcher_Table_Validation.ipynb); \n",
    "    1. [**tutorial on C code MoL implementation**](Tutorial-Method_of_Lines-C_Code_Generation.ipynb)\n",
    "1. NRPy+ supports boundary conditions in Cartesian and singular curvilinear coordinates\n",
    "    1. [**tutorial on general boundary condition driver**](Tutorial-Start_to_Finish-Curvilinear_BCs.ipynb); \n",
    "1. NRPy+ provides support for solving the scalar wave equation stably in curvilinear coordinates, including Cartesian, spherical-like, cylindrical-like, and prolate-spheroidal-like coordinates \n",
    "    1. [**tutorial on Cartesian scalar wave equation in NRPy+**](Tutorial-ScalarWave.ipynb); \n",
    "    1. [**tutorial on C code scalar wave implementation**](Tutorial-Start_to_Finish-ScalarWave.ipynb); \n",
    "    1. [**tutorial on NRPy+ curvilinear coordinates (reference metric) support**](Tutorial-Reference_Metric.ipynb); \n",
    "    1. [**tutorial on scalar wave equation in curvilinear coordinates**](Tutorial-ScalarWaveCurvilinear.ipynb)\n",
    "    1. Einstein Toolkit thorns\n",
    "        1. [**tutorial on `WaveToyNRPy`, for solving the scalar wave equation**](Tutorial-ETK_thorn-WaveToyNRPy.ipynb)\n",
    "        1. [**tutorial on `IDScalarWaveNRPy`, for setting up scalar wave initial data**](Tutorial-ETK_thorn-IDScalarWaveNRPy.ipynb)\n",
    "1. NRPy+ implements a covariant BSSN formulation that supports Cartesian, spherical-like, and cylindrical-like coordinates, and its boundary condition driver automatically sets up correct boundary conditions for any tensor in any orthogonal coordinate system! \n",
    "    1. [**BSSN overview tutorial; contains links to several other tutorials**](Tutorial-BSSN_formulation.ipynb)\n",
    "    1. Einstein Toolkit thorns\n",
    "        1. [**tutorial on `Baikal` & `BaikalVacuum`, for solving Einstein's equations in Cartesian coordinates**](Tutorial-BaikalETK.ipynb)\n",
    "1. NRPy+ contains multiple initial data sets for Einstein's equations of general relativity, and a means to quickly validate they satisfy the Einstein constraint equations on numerical grids\n",
    "    1. [**BSSN initial data validation**](Tutorial-Start_to_Finish-BSSNCurvilinear-Exact_Initial_Data.ipynb); \n",
    "    1. [**static trumpet initial data**](Tutorial-ADM_Initial_Data-StaticTrumpet.ipynb); \n",
    "    1. [**\"UIUC\" spinning black hole initial data**](Tutorial-ADM_Initial_Data-UIUC_BlackHole.ipynb); \n",
    "    1. [**shifted Kerr-Schild initial data**](Tutorial-ADM_Initial_Data-ShiftedKerrSchild.ipynb); \n",
    "    1. [**Brill-Lindquist initial data**](Tutorial-ADM_Initial_Data-Brill-Lindquist.ipynb);\n",
    "    1. [**Fishbone-Moncrief black hole accretion disk initial data**](Tutorial-FishboneMoncriefID.ipynb); \n",
    "    1. [**piecewise-polytrope TOV initial data**](Tutorial-ADM_Initial_Data-TOV.ipynb)\n",
    "1. NRPy+ contains multiple diagnostics for spacetime evolutions\n",
    "    1. [**computing $\\psi_4$ in Cartesian coordinates**](Tutorial-WeylScalarsInvariants-Cartesian.ipynb)\n",
    "    1. [**computing $\\psi_4$ in curvilinear coordinates**](Tutorial-Psi4.ipynb)\n",
    "    1. [**$\\psi_4$ tetrads in curvilinear coordinates**](Tutorial-Psi4_tetrads.ipynb)\n",
    "    1. Einstein Toolkit thorn\n",
    "        1. [**`WeylScal4NRPy`, a `WeylScal4` clone written in NRPy+**](Tutorial-ETK_thorn-Weyl_Scalars_and_Spacetime_Invariants.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='error_analysis'></a>\n",
    "\n",
    "# Step 5: Error analysis & code validation: Confirming numerical errors converge to zero at the expected rate \\[Back to [top](#toc)\\]\n",
    "$$\\label{error_analysis}$$\n",
    "\n",
    "The above code is *really slow*, so to bypass the long wait, results at $24^3$, $48^3$, and $96^3$ are *precomputed* and output to files in the next cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:00.626577Z",
     "iopub.status.busy": "2021-03-07T17:25:00.625669Z",
     "iopub.status.idle": "2021-03-07T17:25:00.628633Z",
     "shell.execute_reply": "2021-03-07T17:25:00.629234Z"
    }
   },
   "outputs": [],
   "source": [
    "# Pasted results, assuming u_offset=1 and Nx=Ny=Nz=24\n",
    "with open(\"output_resolution_24cubed.txt\", \"w\") as file:\n",
    "    file.write(\"\"\"0.00 -16.00 2.999421380013 2.999421380013\n",
    "0.04 -10.70 2.998843001874 2.998843001814\n",
    "0.08 -10.06 2.997108425138 2.997108424880\n",
    "0.12 -9.70 2.994219322036 2.994219321440\n",
    "0.17 -9.45 2.990178477110 2.990178476038\n",
    "0.21 -9.25 2.984989783456 2.984989781772\n",
    "0.25 -9.09 2.978658237475 2.978658235046\n",
    "0.29 -8.95 2.971189932138 2.971189928832\n",
    "0.33 -8.84 2.962592048775 2.962592044465\n",
    "0.38 -8.73 2.952872847425 2.952872841973\n",
    "0.42 -8.64 2.942041655765 2.942041648971\n",
    "0.46 -8.54 2.930108856521 2.930108848127\n",
    "0.50 -8.48 2.917085872939 2.917085863230\n",
    "\"\"\")\n",
    "\n",
    "# Pasted results, assuming u_offset=1 and Nx=Ny=Nz=48 <- required 2 minutes on fast computer\n",
    "with open(\"output_resolution_48cubed.txt\", \"w\") as file:\n",
    "    file.write(\"\"\"0.00 -16.00 2.999855329307 2.999855329307\n",
    "0.04 -11.87 2.999276741878 2.999276741874\n",
    "0.08 -11.25 2.997541537534 2.997541537518\n",
    "0.12 -10.89 2.994651389354 2.994651389316\n",
    "0.17 -10.64 2.990609083291 2.990609083222\n",
    "0.21 -10.45 2.985418514411 2.985418514305\n",
    "0.25 -10.29 2.979084681648 2.979084681494\n",
    "0.29 -10.15 2.971613681053 2.971613680844\n",
    "0.33 -10.04 2.963012697588 2.963012697316\n",
    "0.38 -9.93 2.953289995451 2.953289995108\n",
    "0.42 -9.84 2.942454906957 2.942454906535\n",
    "0.46 -9.76 2.930517820003 2.930517819496\n",
    "0.50 -9.69 2.917490164124 2.917490163524\n",
    "\"\"\")\n",
    "\n",
    "# Pasted results, assuming u_offset=1 and Nx=Ny=Nz=96 <- required\n",
    "with open(\"output_resolution_96cubed.txt\", \"w\") as file:\n",
    "    file.write(\"\"\"0.00 -16.00 2.999963831346 2.999963831346\n",
    "0.04 -13.06 2.999385191594 2.999385191594\n",
    "0.08 -12.45 2.997649830354 2.997649830353\n",
    "0.12 -12.09 2.994759420914 2.994759420911\n",
    "0.17 -11.84 2.990716749579 2.990716749575\n",
    "0.21 -11.65 2.985525711912 2.985525711906\n",
    "0.25 -11.49 2.979191307475 2.979191307466\n",
    "0.29 -11.36 2.971719633092 2.971719633079\n",
    "0.33 -11.24 2.963117874631 2.963117874614\n",
    "0.38 -11.14 2.953394297333 2.953394297311\n",
    "0.42 -11.05 2.942558234689 2.942558234663\n",
    "0.46 -10.97 2.930620075904 2.930620075872\n",
    "0.50 -10.89 2.917591251949 2.917591251911\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We say that our scheme is fourth-order-accurate in the truncation error, so the numerical solution at a given point $(t,x,y,z)$, $u_{\\rm num}$, should satisfy the equation\n",
    "\n",
    "$$u_{\\rm num} = u_{\\rm exact} + \\mathcal{O}(\\Delta x^4) + \\mathcal{O}(\\Delta t^4),$$\n",
    "\n",
    "where $u_{\\rm exact}$ is the exact solution, and $\\mathcal{O}(\\Delta x^4)$ and $\\mathcal{O}(\\Delta t^4)$ are terms proportional to $\\Delta x^4$ and $\\Delta t^4$, respectively. However note that the [CFL condition](https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition) for this PDE requires that $\\Delta x \\propto \\Delta t$, so we can simplify the above expression to\n",
    "\n",
    "$$u_{\\rm num} = u_{\\rm exact} + \\mathcal{O}(\\Delta x^4).$$\n",
    "\n",
    "Therefore, the [relative error](https://en.wikipedia.org/wiki/Approximation_error) between the numerical and the exact solution should be given to good approximation by\n",
    "\n",
    "\\begin{align}\n",
    "E_{\\rm Rel} &= \\left| \\frac{u_{\\rm num} - u_{\\rm exact}}{u_{\\rm exact}}\\right| \\\\\n",
    "&\\propto \\Delta x^4 \\\\\n",
    "\\implies E_{\\rm Rel} &= k \\Delta x^4,\n",
    "\\end{align}\n",
    "where $k$ is the proportionality constant, divided by $u_{\\rm exact}$.\n",
    "\n",
    "Therefore, taking the logarithm of both sides of the equation, we get:\n",
    "\n",
    "\\begin{align}\n",
    "\\log_{10} E_{\\rm Rel} &= \\log_{10} (k [\\Delta x]^4) \\\\\n",
    "&= \\log_{10} ([\\Delta x]^4) + \\log_{10} (k) \\\\\n",
    "&= 4 \\log_{10} (\\Delta x) + \\log_{10} (k)\n",
    "\\end{align}\n",
    "\n",
    "$\\Delta x$ is proportional to `1/Nx`, so if we perform the simulation at twice the resolution (i.e., $\\Delta x\\to \\Delta x/2$), $\\log_{10} E_{\\rm Rel}$ should drop by \n",
    "\n",
    "\\begin{align}\n",
    "4 \\log_{10} (\\Delta x) - 4 \\log_{10} (\\Delta x/2) &= 4 \\log_{10} \\frac{\\Delta x}{\\Delta x/2} \\\\\n",
    "&= 4 \\log_{10} 2 \\approx 1.20.\n",
    "\\end{align}\n",
    "\n",
    "In the below plot we show that when the logarithmic relative error $\\log_{10} E_{\\rm Rel}$ versus time in the $48^3$ case is shifted upward by 1.2 (or for the $96^3$ case by 2.4), we observe perfect overlap with $\\log_{10} E_{\\rm Rel}$ in the $24^3$ case (except at $t=0$ when the numerical solution is set to the exact solution). This is a common way in numerical relativity to present convergence of numerical errors to zero, demonstrating that our code is working as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:00.644094Z",
     "iopub.status.busy": "2021-03-07T17:25:00.643369Z",
     "iopub.status.idle": "2021-03-07T17:25:00.990083Z",
     "shell.execute_reply": "2021-03-07T17:25:00.990528Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Moving from 24^3 grid to 48^3 grid, logarithmic drop in numerical error should be 1.20\n",
      "Moving from 24^3 grid to 96^3 grid, logarithmic drop in numerical error should be 2.41\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f4c6aba0610>]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAArq0lEQVR4nO3deVzV953v8dcHBBTcQdwQcd+NiajRLG3aLG0yWUzqNG22tplx5rYzvb3dZnrTyaO9097OtPd2OpOmt7X7cnvTNgmYrcnENEsNaFyiARV3jfuCCAgIAp/7x++oKBiOwDk/OLyfj8d5cM739zvn9/nJ8vG7m7sjIiLSUlLYAYiISPej5CAiIq0oOYiISCtKDiIi0oqSg4iItNIn7AC6SlZWlufl5YUdhohIj7Ju3brj7j7s4vKESQ55eXmsXbs27DBERHoUM9vbVrmalUREpBUlBxERaUXJQUREWlFyEBGRVpQcRESkFSUHERFpRclBRERaSZh5DiIiiai5GSor4eRJOHj0FPvf3cmhAzs4dnw3FVX7qDh9kG99+XvkTRzdpddVchARibH6eqioCP7An/16sryJ6v3HOHFkB5UVO6k+tZdTDfupbT5MbdIxalNOcKpfJaf61VKVcZoTGWeoS23xoSlAZvD0/nX3KTmIiHQXDQ1w6BAcOACHdtVRXnaEE3t3c/zIDk6d2kVNw37qkw/T2Pc4DRknOZ1RRXVGHZUZDZRnNFGeDj4YGHzh5/ZpgqzaJDJPp5LV0I9J1SMZWj2IzJQsstKzGT54JKOycxk7Zhy54yaQPnFKl9+bkoOISBtqqpo4vKmcQ5sOsG/nNg4d3EX5yb1UnT5IdfMRavqUU9u3kpqMGqoy6jme0czRDKifCExs/XmD6pPJbEgjqymDCTaSRUlDGZ6UxaghIxiZlUN2Zi7Zw8eTPWoSgzNHk5SUHPd7bknJQUR6j8ZGmg4eYV/JTnaWbWHfgZ0cOfEux2sOc7LxOJVJFVSlnqKyXx0VGWc4kgEn+wF9gfEXflRqk5F5Oo3MxnSybDiT0oYyqn82ozJHMiJrLNnZ48geOZFhmblkZ2STmpzaVkTdlpKDiPR8dXWc3n2QXSVb2bWjjP0Hd3Kk8l2O1x+iguOc7FPJiX41HO9/hkP94VRa5H2DIg/AHAbV9WHI6X4MbhzG6ObBzG7IIjttBDnDc8gbm0fu2Ilkj5hA9oARDEgdgJmFdccxp+QgIt2TO1RWUrltH9tKNrNnzzYOHN3N8ap9lJ85QkVSORWpVZSn13G0fxOH+8OZsy0x2ZEHkF6fxNDafgw+ncnwukymNgwnMyOHEUPzyBkzmYlTJjNjci4jBmXSJ0l/Es/Sv4SIxFdDA5V7drNj8xb27NjBwUO7OVaxn/LThznZdJyq5EpOptZwMr2eoxnOsQxwI/hrNSryAAbV9mFobTqDG0YxsTqL+fUjyOqfy4hh4xk7bgqTp01lRt5oBqf3D/Fmey4lBxHptMamMxw/tIs9ZZvZtWMHBw7u4eiJ/ZTXHaGi8Tgnk4I/+Cf6BR23tS2b3wdGHkBGfRKDa9MYWN+fgWdyuKIyi8zTo8gePJaRw8eTN3Ea06dOZErOSNL69Kw2/J5GyUFE2uTuHK8+wvata9hWtoX9B/Zw5MQBjtUc4kRjORVWde4PfkW/5uB/92clA8OCIZmZNckMqk1jYP0AJtSOYc6JLAanjWDogByys8eRM3Yi46dOZdqE0Qwd2C+s25WLKDmI9GKnKo6w6e1VvFP6Ntv3bWZf9R4O+kEOpVVwYEAdtal+4RsyYEhS5A9+TRo5VQOY2jiG/p7FwJSRDBk4hmFZeYwaM5lx06czceoYhmYmkaSFenocJQeRROXO6QP72bj6LTZv2cCuQ2Xsq9nNITvM4b7BH//yjBZ//PvDgBQYXZHCyPKBzNyfS6aPIbvfJEYMm8zInCmMmTyFETNyGDEmhbS0S19aer5ulxzM7Argh0B/YA9wn7tXhRqUSHfkTv3ufbzz6irKtr3NnqNbOVC7hyPJhznUr4IDg05zcAA0n/1fezakNMHIk6kMrx7I/JN5DLMxjOg/mdxRs5g6Yx5T5s1g5JgUksOdfyXdQLdLDsBPgC+6++tm9ingS8A/hRyTSGgaKyooeWMlb68tYsf+EvbX7eBwymEODKhi91CnLgVIB/KCsfrDqlPJPjWQSSfGc23lWEYNnEzemNlMv2IeV101jSFDuuOvvXQ33fGnZDLwRuT5y8BLKDlIgjtTX0fZ2jdZu+pNtu7ewLvVOziUdJAD/SvZO7iJhj4Ev615kNpojDqZwfC6cdx8bBxjMqczcfwcZs+Zz7wZk+jfT+090nndMTlsAu4ECoElwJhLnWhmS4GlALm5ufGITaTD6s+cZseOdawtWsnmbW+zu2IbB9jPgfST7B90hqazzT/DIH0QjK5IZ3hVDrOq8hg9eDqTJs4j/9rrmDd7PH2S1cMrsWXu3v5ZXX1RsxXAiDYOPQJsBf6DYDHaZ4DPuntme5+Zn5/va9eu7dI4RTriZF0F69f9iaLiP1Gyfz27G/dwoO8JDg1ouGC458DTMOZEX7KrsxjWlMuI/tMYNy6fKxdcT/6iaWRkJO7SDNJ9mNk6d8+/uDyUmoO739jOKTcDmNlk4LbYRyRy+erO1FFa9iYr31jBxp1r2H56G7vSj3J4QMO5c/pnwMRjKcw8lsn7z+SQnTaV3NFXMeOq67nihivIHtmHBF6eR3qwbtesZGbZ7n7UzJKArxKMXBIJTWNzI9v2vs2br/8n67asYlv1FnanHeTdgXXnRgKlZsCkumTm7M1iRP14cgbMZdqUG7jyjg8wce5gUlLCvQeRy9XtkgPwMTP7TOT508DPwwxGeg93591j21n155dYvWElZSdK2ZW8n92DqoIOYSCpL4yvNcYfHsw1O6cyuu9sJo1/H1fecDPTrh1FupqCJEF0u+Tg7v8O/HvYcUhiO32mjjXrX2bFn56j5OB6dvhudg88yam05uCEPjAmDcYeHcDsvdMYlTyTiTmLuGLRLcy8aQpDMtUhLImt2yUHka52pukMJWVv8MorhazeWcwm386OwVU0RiZ6ZabDxCP9uOnAWEb6dPKGLeCKeTcz65NXMWJMivoEpFdScpCE0tTcRNnut1j5RiF/3vQGJQ3b2DroBPWRn/Qh/WDagQzu2TOLvL4LmDP7dubd/0HGzcjQ+j8iLSg5SI/l7uzc/w6rip/mzxtfZ331ZjYPOH5usbj+qTC9vC+37ZpKTlI+0ybdxoI7b2XmwoHqIBZph5KD9Ajuzr5j21lT/BQrN77K2vJSStKPUNk36CPoazCtNpUP7BrHyMa5TMi5hbm33MncL2UxZEjIwYv0QEoO0i2VVx9ldfEfKNrwEqsOb2BjyiGOpzcCwR4B00/3YeHuHIaduoIxWTdxxfvuYu4nxzB+POojEOkCSg4SuoamBjaWrqB49dP8eUcx65p3s3tAHQBJzTC5Pomrdg1nSMVMRvT/INPmLWbOJycx50rTstEiMaLkIHHl7uw5tIXVRb9nZemfKK7cxKb+5zuMRzbDlP2DmXNoDlnJNzB++hJmfugKFlxtDBsWbuwivYmSg8RUZW0Fa94qoOjt53lz3zrWpx7geL+geahfI8yoTOPGzVMYUnM1OSPvYtpNtzDv4X5MmYJGD4mESMlBukxjcyObtq2kuPhJ/rxtJWvrd7B9YM25xeYmNyRx5e7hDD4+m+HptzBpwRKu/GQOV10FGRnhxi4iF1JykA47cHwXq9/8PStLVlBcXsI7GceoTQmGkWY2w/QjA5j2dj6ZzdcxbtJHmH7LAhYsSmb06JADF5F2KTlIVGrqT7Fu3bMUr32WlXvfYn3yPg5mBKuPpjTBjFMpvH/rOAZXzWfUsNuZ+v7byf/4AGbMgD76KRPpcfRrK5dU11DLH57+JstW/5RVA4+c24xmXKMxbW8WC4/MZFjaTUy8aglXfHQi8+bBoEHhxiwiXUPJQVp5u+Q/efx3j/KH5jVUpTWT12zcWjSTzPrryR17D9NvuY75n0shL09zCkQSlZKDAHDy1HF+/duv8uMtT1AyuJI0g+u2DmNs+cPc+NAj3Pkv/enXL+woRSRelBx6MXfn9Tf/Lz945l94NnUTp1OCmcdLXrmRObP+J/d+Zx7jx4cdpYiEQcmhFzp0dCc/+c1X+PnB59g9oI6ByXDDhrGMb/h7bv3M33PL91NJTg47ShEJk5JDL9HY3Mjzf3yMH/zp33ml/16akmB+eV+uXrWE+e//Fh/76QSGDw87ShHpLkJLDma2BPgaMA2Y7+5rWxz7CvAw0AR81t1fCiXIBLB991p++Nv/zm+qXuVoeiPDDW4rns6k9K9w1+c/zjXXJalTWURaCbPmUArcDfyoZaGZTQfuBWYAo4AVZjbZ3ZviH2LPVFt/ij88/U1+9NbPKB58lORmuGbfQD68716uXfzP/OWj2QwcGHaUItKdhZYc3H0LgLX+b+udwBPuXg/sNrMdwHygOL4R9jzr33mJx3//KE82r6UqrZlxzUnc/doCpud8nb/8yi3MmhV2hCLSU3THPofRwKoWr/dHyloxs6XAUoDc3NzYR9YNNTU18uNffoHHS39B6aAq+gLXlg1nfMXD3PjJR7jj2+la1lpELltMk4OZrQBGtHHoEXdf3tnPd/dlwDKA/Px87+zn9TRbtq7kgR/eybrBJ5hZ24cla27myiu/xce+exV5eWFHJyI9WUyTg7vf2IG3HQDGtHidEymTiMbmRr71/Qf5xrH/x4BUeODlj/KxL/+Kmx/XEFQR6RrdsVnpGeC3ZvZdgg7pScBb4YbUfZSUvc4DP1nMxgEVfGDrEK4dUMBXVryPvn3DjkxEEkmYQ1kXA48Bw4DnzWyDu9/i7pvM7PfAZqAR+IxGKgVbaX7zhw/yrSO/Y4jBfc99jKX/+iuu/0B3zO8i0tOFOVqpACi4xLFvAt+Mb0Td1/otr/Dgz+9hU0YlH9o0lBlJhTz6ynUajioiMRN1cjCzIQTNPHXAHndvjllUAkB9Yz1fX3Yf3z7yFNnN8MDTH+cjj/6CO+5JCTs0EUlw75kczGwQ8BngY0AqcAzoCww3s1XAD9z91ZhH2Qut3vSfPPTrJWztV8UdG4cyomY5//zytWRnhx2ZiPQG7dUcngR+BVzn7idbHjCzucADZjbe3X8ao/h6nbozdXz1xx/je0eWM7oBHlx+Pzd8/mc89FcpWuZCROLmPZODu9/0HsfWAeu6PKJebGXpC3zyt/eyI62ae9Zn0vfocr7+x2s0Z0FE4i4p7AAEahpq+OwPbuf6J2+jsbaah379AIuuOcgvi5QYRCQc7fU57AYcsMjXVqdEyr/n7v/R9eElvj+9s5y/+t397E49xcfeyuT0nuV8oeAarYMkIqFqr1lpXLwC6W2q6qv48k8/yo+OvcjEavjU8ocYde8yvlqQqrWQRCR07Q5lNbNk4Ffufl8c4ukVXtzwJEuf/AQHkmt4cFUmx8uW86nfXsM114QdmYhIoN0+h8js5LFmlhqHeBJaRV0Fn3z8Jj68fAkZJ2v4q599goycA/xukxKDiHQv0U6C2wW8aWbPADVnC939uzGJKgE9+/YT/M3Tn+JoUh1//WYmu95Zzh0/v4bbbgs7MhGR1qJNDjsjjyRgQOzCSTw1DTUsXXY7vy1/lVnlsLjwE5Qv+CFPbE4jKyvs6ERE2hZVcnD3rwOYWf/I61OxDCqR/J/ffp7flr/Kp1/PonRtIVc/fg33348mtIlItxZVcjCzmcCvgaGR18eBB919UwxjSwhPbHqO2dXGrqa9/Lo0nV66YZ2I9DDRToJbBnze3ce6+1jgC8CPYxdWYjhy6gjrMw4yt2wMy19WYhCRniPa5JDRcoE9d38NyIhJRAlkeemTuEHu4Q+SqrFeItKDRD1aycz+iaBpCeB+ghFM8h4K1/6G8Segf/YdYYciInJZoq05fIpgx7angaeArEiZXEJVfRWvlK/hrjJIuW5R2OGIiFyWdpNDZIb00+7+WXe/yt3nuvvn3L2iMxc2syVmtsnMms0sv0V5ppm9amanzOz7nblGmF7Y/gINNDGvbBQTFmoTBhHpWaKdId0c2finK5UCdwNvXFR+Gvgn4ItdfL24KiwrZFhNEvX73s/MmWFHIyJyeaLtczgFlJjZy1w4Q/qzHb2wu28BsIsG/Lt7DbDSzCZ29LPDVt9YzwvbnuejZc1sSL2WB8eGHZGIyOWJNjk8HXl0K2a2FFgKkNuNxom+svsVqs+cYnEZ/G7KQk14E5EeJ9pVWT/h7jdc7oeb2QpgRBuHHnH35Zf7eRdz92UEczDIz89va7+JUBSWFdK/OYX5u1J59pNqUxKRnqfd5ODuTZFO40HuXnk5H+7uN3Y8tJ6pqbmJ5VuX86EDGWxsuorps6OtnImIdB+h9TkkquL9xRytOcrit4xiFnKtdnQTkR4otD4HM1sMPEYwf+J5M9vg7rdEju0BBgKpZnYXcLO7b+7K68dKYVkhKdaHv9jWyL0s4r8oOYhID9TeHtID3b3K3X/ZxrFO9QC7ewFQcIljeZ357LC4OwVlBXyQ8Qys38ae4VeTmRl2VCIil6+9eQ6vnX1iZq9cdKywq4Pp6UqOlrCrYheLd6exO20qObOHhh2SiEiHtJccWg7CvPgvnQZoXqSwrBDDuP2V/bxxZiGz1KQkIj1Ue8nBL/G8rde9XkFZAQuz5jByXwUrm5UcRKTnaq9DOtvMPk9QSzj7nMjrYTGNrIfZc3IPGw5v4DtDPgq8TRGL+FslBxHpodpLDj/m/J7RLZ8D/CQmEfVQhWWFANy11ahLG0RZ/TSmTQs3JhGRjnrP5HB272hpX0FZATOzZzLxqU28PehqJgxKIj097KhERDom2v0c5D0cqznGyndXsnjcrVBayp8b1d8gIj2bkkMXeHbbszR7M4vr8sCd504sUnIQkR5NyaELFJQVMHbQWOa8cxQ3YxULtIeDiPRoUSUHMxtuZj81sz9GXk83s4djG1rPcKrhFC/vfJm7pt6FrVrFyVEzqGagag4i0qNFW3P4BfASMCryehvwuRjE0+O8uONF6pvqWTzlTigupmzoIvr2hYk9dqsiEZHok0OWu/8eaAZw90agKWZR9SAFZQVk9svkmppMqKxkZfNCpk+H5OSwIxMR6bhok0ONmWUSmRVtZlcDl7W3QyJqaGrg+W3Pc8eUO+iz6i0Alh9dpP4GEenxol2y+wvAM8AEM3uTYHb0R2IWVQ/x2p7XqKyvZPHUxfBMIc1DM3nz2CTuUn+DiPRwUSUHd19nZu8DphAsnbHV3c/ENLIeoGBLARkpGdw4/kYo+jInpiyEYlNntIj0eNGOVnoH+DJw2t1LlRig2ZuD7UAnfoh+1XVQVsb2zIUASg4i0uNF2+dwO9AI/N7M1pjZFzu72U9P99aBtzh06lDQpLR6NQBFLGTIEBg5MuTgREQ6Kark4O573f3b7j4X+DgwG9gd08i6uYItBfRJ6sOtk26FoiJITuaFo/OYNQtMO12ISA8X9QxpMxtrZl8GngCmEjQzdZiZLTGzTWbWbGb5LcpvMrN1ZlYS+fqBzlwnFs5uB3pD3g0M6TcEiovx2bNZs6W/mpREJCFE2+ewmmC/52RgibvPd/f/3clrlwJ3A29cVH4cuN3dZwEPAb/u5HW63JbjW9h+YnvQpNTUBKtXUz1zEdXV6m8QkcQQ7VDWB919a1de2N23ANhFbTDu/naLl5uAfmaW5u71XXn9zijYUgDAHVPugNJSOHWKndlBZ7TmOIhIInjP5GBm97v7b4DbzOy2i4+7+3djFlngHmD9pRKDmS0FlgLk5savf7xwayELRi9g9MDRUPwsAKtMyUFEEkd7zUoZka8D2nj0b+/DzWyFmZW28bgzivfOAP4V+JtLnePuy9w9393zhw2Lz66l+yr3sfbg2qBJCYLO6OHD+fP+ceTmwqBBcQlDRCSm2tsJ7keRpyvc/c2Wx8zsmvY+3N1v7EhQZpZD0MfxoLvv7MhnxMq57UCn3hUUFBfDwoWUlGrym4gkjmhHKz0WZVmnmdlg4HngHy9OSN1B4dZCpmVNY0rWFDh6FHbsoHH+IsrK1KQkIomjvT6HhcAiYJiZfb7FoYEEI5c6zMwWEySYYcDzZrbB3W8B/g6YCDxqZo9GTr/Z3Y925npdoby2nNf3vM4/XPMPQcGqVQC8O3ohjY0aqSQiiaO90UqpBH0LfQj6Gc6qopML77l7AUHT0cXl3wC+0ZnPjpXntj1Hkzedb1IqKoKUFNb6XEDJQUQSR3t9Dq8Dr5vZL9x9b5xi6rYKtxaSMzCH/FGROXvFxXDllWzY2o8+fWDq1HDjExHpKtHOc6g1s+8AM4C+ZwvdvdvNXo6V2jO1vLTjJR6+8uFgbsaZM7BmDSxdSmkpTJ4MqalhRyki0jWi7ZD+v0AZMA74OrAHWBOjmLqll3a8RF1j3fkmpY0boa4OFi2ipERNSiKSWKJNDpnu/lPgjLu/7u6fAnpNrQGCJqUhfYdw/djrg4LiYgBOzVrInj1KDiKSWKJtVjq7f8OhyEzpg8DQ2ITU/ZxpOsOzW5/l9im3k5KcEhQWFUFODqWVYwAlBxFJLNEmh2+Y2SCC7UIfIxjK+t9iFlU388beN6g4XcFdU+46X3h28ltJ8FJzHEQkkUS7TehzkaeVwA2xC6d7KiwrpF+fftwy8Zag4OBB2LsXPvc5SkogIwPy8kINUUSkS7U3Ce4xwC913N0/2+URdTPuTuHWQm6ecDPpKelBYaS/gYULKVke1BqSot4ZQ0Sk+2uv5rA2LlF0Y2sPrmV/1X6+cUOLeXnFxZCWhs+5kpISWLw4vPhERGKhvUlwv2z52szS3b02tiF1L4VlhSRbMn8x+S/OFxYVQX4+RypSKS9Xf4OIJJ5od4JbaGabCeY6YGZXmNkPYhpZN1FQVsD1Y68nMz0zKKivh3XrLuiM1kglEUk00baUfw+4BSgHcPeNwPUxiqnb2Hp8K1uObzm/dwPA+vXQ0HBu8hsoOYhI4om6G9Xd911U1NTFsXQ7rfZugAs6o0tLYfhwiNM+QyIicRPtPId9ZrYIcDNLAf4rsCV2YXUPBWUFzB05lzGDxpwvLC6GceNgxAhKStTfICKJKdqaw98CnwFGAweAOcCnYxRTt3Cw+iCrD6y+sEnJPeiMXriQpibYtElNSiKSmKKdBHccuO/sazMbQpAcvhmjuEK3vGw5AIuntUgO+/YFE+AWLmTXrmDdPSUHEUlE71lzMLMxZrbMzJ4zs4fNLMPM/hewFciOT4jhKCgrYNLQSUzLmna+sKgo+LpoEaWlwVM1K4lIImqvWelXBIvsPQbMJJgUNxqY7e7/tTMXNrMlZrbJzJrNLL9F+Xwz2xB5bIxsJxpXJ0+f5NU9r7J46uJg74azioshPR1mz6akBMxgxox4RyciEnvtNSsNdfevRZ6/ZGZLgPvcvbkLrl0K3A38qI3yfHdvNLORwEYze9bdG7vgmlF5ftvzNDY3XtikBEFymD8f+vShpATGjw/WVRIRSTTtdkib2RAzG2pmQwnmOQxq8brD3H2Lu29to7y2RSLoy3us7RQrBWUFjOw/kvmj558vrKuDt9+GhQsBtMGPiCS09moOg4B1QIu2FdZHvjowPhZBmdkC4GfAWOCBS9UazGwpsBQgNze3S65dd6aOF3e8yAOzHyDJWuTOtWuhsREWLuT0adi+HZYs6ZJLioh0O+2trZTXmQ83sxXAiDYOPeLuy9/juquBGWY2Dfilmf3R3U+3cd4yYBlAfn5+l9QwVuxaQc2ZmtZNSmc7oxcuZMsWaG5WzUFEEle0k+A6xN1v7OT7t5jZKc53hsdcQVkBg9IG8f689194oLgYJk2CrCxKXgiKlBxEJFG1t5/DboLmI2vj8NlyB77n7v/RFQGZ2ThgX6RDeiwwFdjTFZ/dnsbmRp7Z+gy3Tb6N1OTU8wfcg+Tw4Q8DQX9DamqQK0REElF7zUrjYnXhyBDVx4BhwPNmtsHdbwGuBf7RzM4AzcCnI5PwYu7Nd9+kvK78wlnRALt2wdGj5zqjS0th2jToE9N6l4hIeEL78+buBUBBG+W/Bn4d/4iCJqW05DQ+NPFDFx44u9jeokVAUHO4oddtlioivUl7M6SfNbPbI4vtXXxsvJn9DzP7VOzCix93p7CskJsm3ET/1P4XHiwqggEDYPp0KirgwAH1N4hIYmtvnsNfA9cBZWa2xsxeMLM/mdkugslr69z9ZzGPMg42HN7A3sq9rZuUIKg5LFgAycnnls1QchCRRNZen8Nh4MvAl80sDxgJ1AHbEm270IKyApIsidsn337hgepqeOcd+OpXAc5t8KM1lUQkkUW7TehwYChQDxxKtMQAwcY+1+Zey7CMi3buWbMmmNTQYmb0oEGQkxNCkCIicdLeUNY5wA8JZkofiBTnmNlJglFE6y/x1h5l54mdlBwt4d9u+bfWB892Rl99NXB+2Qxra3CviEiCaG+00i+Av4nMWD7HzK4Gfg5cEaO44qqgLBg0dcF2oGcVFcH06TB4MO7BMNaPfzy+8YmIxFt7zUoZFycGAHdfBSTMeqSFZYXMGTGHvMF5Fx5oboZVq84NYd2/Hyor1d8gIomvvZrDH83seYJ9HfZFysYADwIvxjKweDly6ghF+4r42vu/1vrgtm1w4sQF/Q2gkUoikvjaG630WTP7MHAnwSY/EPQ9PO7uL8Q6uHhYvnU5jrfdpHS2v+Gi5KCag4gkunZnSLv7H4E/xiGWUBSWFTJ+yHhmZbdRHSguhiFDYMoUIOhvyMkJikREElmHl88ws2XuvrQrgwnDo+97lOO1xy/cDvSsoqJglFJS0DVTUqJag4j0Du0NZb3Ubm8G3Nr14cTf1TlXt33g5EnYvBnuvReAM2dgyxa4+eb4xSYiEpb2ag7HgL1cuGT32aW6s2MVVLewenWwVHekv2H7dmhoUGe0iPQO7SWHXcAH3f3diw+Y2b42zk8cxcVBc9L8YB9prakkIr1Je/Mcvgdcqvv1210bSjdTVBRkggEDgKC/ITkZpk4NOS4RkTh4z+Tg7o+7+8ZLHHssNiF1A83NQbNSpEkJguQwaRL07RtiXCIicRLVaCUzu7uN4kqgxN2Pdm1I3cDmzVBVdW5mNATJYe7cEGMSEYmjqFZlBR4GfgLcF3n8GPgH4E0ze6AjFzazJWa2ycyazSy/jeO5ZnbKzL7Ykc/vlKKi4Guk5lBTE+wUqv4GEektok0OfYBp7n6Pu98DTCcYtbSAIEl0RClwN/DGJY5/l7Am3xUXw7BhMGECAJs2BcWa4yAivUW0k+DGuPuRFq+PRspOmNmZjlzY3bcAbU4+M7O7gN1ATUc+u9OKioJaQyQ2rakkIr1NtDWH18zsOTN7yMweAp6JlGUAJ7syIDPrT1Ab+XoU5y41s7VmtvbYsWNdE0B5ebDg3kWd0enpMH5811xCRKS7i7bm8BmCJqBrI69/CTzl7g7ccKk3mdkKYEQbhx5x9+WXeNvXgH9z91NtLmnRgrsvA5YB5Ofn+3ueHK1Vq4KvLTqjS0thxoxzq2iIiCS8qJKDu7uZrQQaCPoa3ookhvbed2MHYloAfMTMvg0MBprN7LS7f78Dn3X5ioqgTx/IP99HXlICt90Wl6uLiHQL0Q5l/UvgO8BrBEtnPGZmX3L3J7s6IHe/rsV1vwaciltigKAzes6coB0JOHo0eKi/QUR6k2gbSh4B5rn7Q+7+IDAf+KfOXNjMFpvZfmAh8LyZvdSZz+sSjY2tJr9p2QwR6Y2i7XNIumiyWznRJ5Y2uXsBUNDOOV/rzDUuW0kJ1Na2mvwGSg4i0rtEmxxejPzP/v9FXn8USIid4C5w0eQ3CJJDVhZkJ/YatCIiF4i2Q/pLZnYPcE2kaFnkf/6JpbgYRo6E3NxzRSUlQa2hnYFTIiIJJeqd4Nz9KeCpGMYSvuLioEkpkgmam4PZ0Q8/HHJcIiJx9p79BmZWbWZVbTyqzawqXkHGxZEjwQJKLZqU9uwJ1lVSf4OI9DbvWXNw9wHxCiR0xcXB1zY6o7Wmkoj0Nprze1ZREaSmwlVXnSs6mxxmzAgpJhGRkCg5nFVcHGzYkJZ2rqi0FMaNO7cZnIhIr6HkANDQAGvXXtDfAOdHKomI9DZKDgAbNsDp0xckh/p62LpV/Q0i0jspOUCbndFlZdDUpJqDiPROSg4QdEbn5sKoUeeKtKaSiPRmSg5wfvJbCyUlkJICkyeHFJOISIiUHPbvh3372uyMnjo1SBAiIr2NksPZ/gaNVBIROUfJobgY+vULNviJqKwMKhNKDiLSWyk5bNkSbAnaov1IndEi0ttFvSprwnrhhaCq0ILWVBKR3i60moOZLTGzTWbWbGb5LcrzzKzOzDZEHj+McSAwePAFRaWlMHDgBds6iIj0KmHWHEqBu4EftXFsp7vPiW8455WUBLUGbfAjIr1VaDUHd9/i7lvDuv6luJ9PDiIivVV37ZAeZ2Zvm9nrZnZdPC988CBUVKgzWkR6t5g2K5nZCmBEG4cecffll3jbISDX3cvNbC5QaGYz3L3VznNmthRYCpDbRR0EGqkkIhLj5ODuN3bgPfVAfeT5OjPbCUwG1rZx7jJgGUB+fr53LtqARiqJiHTDZiUzG2ZmyZHn44FJwK54Xb+kBEaOhMzMeF1RRKT7CXMo62Iz2w8sBJ43s5cih64H3jGzDcCTwN+6+4l4xaVlM0REQhzK6u4FQEEb5U8BT8U/omD/hs2b4e/+Loyri4h0H92uWSlMO3YEO8Cp5iAivZ2SQwvqjBYRCSg5tFBSAklJMH162JGIiIRLyaGF0lKYODFYwVtEpDdTcmhBI5VERAJKDhG1tUGHtPobRESUHM7ZvDlYdE81BxERJYdztKaSiMh5Sg4RJSXQty9MmBB2JCIi4VNyiCgpCYawJieHHYmISPiUHCI0UklE5DwlB+D4cTh8WMlBROQsJQfUGS0icjElB7SmkojIxZQcCGoOQ4cGm/yIiIiSA3C+M9os7EhERLqHXp8c3IOag/obRETO6/XJYe9eqK5Wf4OISEth7iG9xMw2mVmzmeVfdGy2mRVHjpeYWd9YxaGRSiIirYW2hzRQCtwN/KhloZn1AX4DPODuG80sEzgTqyA0UklEpLXQkoO7bwGw1r3ANwPvuPvGyHnlsYyjpATGjoWBA2N5FRGRniXMmsOlTAbczF4ChgFPuPu32zrRzJYCSwFyc3M7dLFZs6CDbxURSVgxTQ5mtgIY0cahR9x9+SXe1ge4FpgH1AKvmNk6d3/l4hPdfRmwDCA/P987EuNXvtKRd4mIJLaYJgd3v7EDb9sPvOHuxwHM7AXgKqBVchARkdjojkNZXwJmmVl6pHP6fcDmkGMSEelVwhzKutjM9gMLgecjfQy4ewXwXWANsAFY7+7PhxWniEhvFOZopQKg4BLHfkMwnFVERELQHZuVREQkZEoOIiLSipKDiIi0ouQgIiKtmHuH5o51O2Z2DNjbwbdnAce7MJyeQPfcO+ieE19n73esuw+7uDBhkkNnmNlad89v/8zEoXvuHXTPiS9W96tmJRERaUXJQUREWlFyCCwLO4AQ6J57B91z4ovJ/arPQUREWlHNQUREWlFyEBGRVnpVcjCzD5nZVjPbYWb/2MbxNDP7XeT4ajPLCyHMLhXFPV9vZuvNrNHMPhJGjF0pivv9vJltNrN3zOwVMxsbRpxdKYp7/lszKzGzDWa20symhxFnV2rvnlucd4+ZuZn1+KGtUXyfP2FmxyLf5w1m9leduqC794oHkAzsBMYDqcBGYPpF53wa+GHk+b3A78KOOw73nAfMBn4FfCTsmONwvzcA6ZHn/6WXfI8Htnh+B/Bi2HHH+p4j5w0A3gBWAflhxx2H7/MngO931TV7U81hPrDD3Xe5ewPwBHDnRefcCfwy8vxJ4INmZnGMsau1e8/uvsfd3wGawwiwi0Vzv6+6e23k5SogJ84xdrVo7rmqxcsMoKePQonmdxngn4F/BU7HM7gYifaeu0xvSg6jgX0tXu+PlLV5jrs3ApVAZlyii41o7jmRXO79Pgz8MaYRxV5U92xmnzGzncC3gc/GKbZYafeezewqYIwnzkZh0f5s3xNpMn3SzMZ05oK9KTmInGNm9wP5wHfCjiUe3P1xd58A/APw1bDjiSUzSyLYTfILYccSZ88Cee4+G3iZ860gHdKbksMBoGUmzYmUtXlOZP/qQUB5XKKLjWjuOZFEdb9mdiPwCHCHu9fHKbZYudzv8RPAXbEMKA7au+cBwEzgNTPbA1wNPNPDO6Xb/T67e3mLn+efAHM7c8HelBzWAJPMbJyZpRJ0OD9z0TnPAA9Fnn8E+JNHenp6qGjuOZG0e79mdiXwI4LEcDSEGLtaNPc8qcXL24DtcYwvFt7znt290t2z3D3P3fMI+pbucPe14YTbJaL5Po9s8fIOYEunrhh2L3yce/xvBbYR9Po/Ein7HwQ/OAB9gT8AO4C3gPFhxxyHe55H0H5ZQ1BL2hR2zDG+3xXAEWBD5PFM2DHH4Z7/HdgUud9XgRlhxxzre77o3Nfo4aOVovw+fyvyfd4Y+T5P7cz1tHyGiIi00pualUREJEpKDiIi0oqSg4iItKLkICIirSg5iIhIK0oOIh1gZoPN7NOR56PM7MmwYxLpShrKKtIBkeXcn3P3mWHHIhILfcIOQKSH+hdggpltIJhxPM3dZ5rZJwiWp8gAJgH/i2CJ5QeAeuBWdz9hZhOAx4FhQC3w1+5eFu+bELkUNSuJdMw/AjvdfQ7wpYuOzQTuJph9/k2g1t2vBIqBByPnLAP+3t3nAl8EfhCPoEWipZqDSNd71d2rgWozqyRYLROgBJhtZv2BRcAfWmwXkhb/MEUuTclBpOu1XOm1ucXrZoLfuSTgZKTWIdItqVlJpGOqCZaGvmwe7My228yWAFjgiq4MTqSzlBxEOsDdy4E3zayUjm0YdB/wsJltJFhJM6ZbPopcLg1lFRGRVlRzEBGRVpQcRESkFSUHERFpRclBRERaUXIQEZFWlBxERKQVJQcREWnl/wN89lsBIWTq8wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# from https://stackoverflow.com/questions/52386747/matplotlib-check-for-empty-plot\n",
    "import numpy\n",
    "time__arr = []\n",
    "lgerr_arr = []\n",
    "for i in [24, 48, 96]:\n",
    "    t, log10error, num, exact = numpy.loadtxt(fname='output_resolution_'+str(i)+'cubed.txt', delimiter=' ', unpack=True)\n",
    "    time__arr.append(t)\n",
    "    lgerr_single_dataset = []\n",
    "    if i != 24:\n",
    "        print(\"Moving from 24^3 grid to \"\n",
    "              +str(i)+\"^3 grid, logarithmic drop in numerical error should be \"+'%.2f' % (4*np.log10(i/(24.0))))\n",
    "    for log10err_onept in log10error:\n",
    "        lgerr_single_dataset.append(log10err_onept + 4*np.log10(i/(24.0)))\n",
    "    lgerr_arr.append(lgerr_single_dataset)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.set_xlabel('time')\n",
    "ax.set_ylabel('log10(|Relative Error|)')\n",
    "ax.plot(time__arr[0], lgerr_arr[0], color='b')\n",
    "ax.plot(time__arr[1], lgerr_arr[1], color='r')\n",
    "ax.plot(time__arr[2], lgerr_arr[2], color='g')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='student_exercises'></a>\n",
    "\n",
    "# Step 6: Exercises for students \\[Back to [top](#toc)\\]\n",
    "$$\\label{student_exercises}$$\n",
    "\n",
    "1. Adjust the above code to make it run twice as fast on the same numerical grid, while generating exactly the same results (stored to files above). *Bonus*: Can you make it run any faster than twice as fast (while still being written in \"pure\" Python, using `NumPy`)?\n",
    "1. How much should the (absolute value of the) relative error `|rel_error|` drop if we were to double the resolution, assuming instead 6th-order convergence? How will this adjust the `log10(|rel_error|)`?\n",
    "1. Why did we add a nonzero constant offset to the exact solution? (*Hint: Start from the definition of relative error.*)\n",
    "1. What will happen to the convergence order if we continue the simulation for a much longer time, say to $t=2$? Why? (*Hint: Convergence order will plummet!*)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 7: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-Solving_the_Scalar_Wave_Equation_with_NumPy.pdf](Tutorial-Solving_the_Scalar_Wave_Equation_with_NumPy.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": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:00.994633Z",
     "iopub.status.busy": "2021-03-07T17:25:00.993993Z",
     "iopub.status.idle": "2021-03-07T17:25:05.699479Z",
     "shell.execute_reply": "2021-03-07T17:25:05.700255Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Solving_the_Scalar_Wave_Equation_with_NumPy.tex, and\n",
      "    compiled LaTeX file to PDF file Tutorial-\n",
      "    Solving_the_Scalar_Wave_Equation_with_NumPy.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Solving_the_Scalar_Wave_Equation_with_NumPy\")"
   ]
  }
 ],
 "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": 2
}