Tutorial-ScalarWave.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 for the Scalar Wave Equation in Cartesian Coordinates\n",
    "\n",
    "## Authors: Zach Etienne & Thiago Assumpção\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "## This module generates the C Code for the Scalarwave in Cartesian coordinates and sets up either monochromatic plane wave or spherical Gaussian [Initial Data](https://en.wikipedia.org/wiki/Initial_value_problem). \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 ([right-hand-side expressions](#code_validation1); [initial data expressions](#code_validation2)). In addition, all expressions have been validated against a trusted code (the [original SENR/NRPy+ code](https://bitbucket.org/zach_etienne/nrpy)).\n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [ScalarWave/ScalarWave_RHSs.py](../edit/ScalarWave/ScalarWave_RHSs.py)\n",
    "* [ScalarWave/InitialData.py](../edit/ScalarWave/InitialData.py)\n",
    "\n",
    "## Introduction: \n",
    "### Problem Statement\n",
    "\n",
    "We wish to 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 \\text{,}$$\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.\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\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.\n",
    "\n",
    "### The Method of Lines\n",
    "\n",
    "Once we have 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**. In [the next module](Tutorial-Start_to_Finish-ScalarWave.ipynb), 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\").\n",
    "\n",
    "### Basic Algorithm\n",
    "\n",
    "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, with NRPy+-based components highlighted in <font color='green'>green</font>. We will review how NRPy+ generates these core components in this module.\n",
    "\n",
    "1. Allocate memory for gridfunctions, including temporary storage for the RK4 time integration.\n",
    "1. <font color='green'>Set gridfunction values to initial data.</font>\n",
    "1. Evolve the system forward in time using RK4 time integration. At each RK4 substep, do the following:\n",
    "    1. <font color='green'>Evaluate scalar wave RHS expressions.</font>\n",
    "    1. Apply boundary conditions.\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\".** In the following sections we will \n",
    "\n",
    "1. Use NRPy+ to cast the scalar wave RHS expressions -- in finite difference form -- into highly efficient C code, \n",
    "    1. first in one spatial dimension with fourth-order finite differences,\n",
    "    1. and then in three spatial dimensions with tenth-order finite differences.\n",
    "1. Use NRPy+ to generate monochromatic plane-wave initial data for the scalar wave equation, where the wave propagates in an arbitrary direction.\n",
    "\n",
    "As for the $\\nabla^2 u$ term, spatial derivatives are handled in NRPy+ via [finite differencing](https://en.wikipedia.org/wiki/Finite_difference).\n",
    "\n",
    "We will sample the solution $\\{u,v\\}$ at discrete, uniformly-sampled points in space and time. For simplicity, let's assume that we consider the wave equation in one spatial dimension. Then the solution at any sampled point in space and time is given by\n",
    "$$u^n_i = u(t_n,x_i) = u(t_0 + n \\Delta t, x_0 + i \\Delta x),$$\n",
    "where $\\Delta t$ and $\\Delta x$ represent the temporal and spatial resolution, respectively. $v^n_i$ is sampled at the same points in space and time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "1. [Step 1](#initializenrpy): Initialize core NRPy+ modules\n",
    "1. [Step 2](#rhss1d): Scalar Wave RHSs in One Spatial Dimension, Fourth-Order Finite Differencing\n",
    "    1. [Step 2.a](#ccode1d): C-code output example: Scalar wave RHSs with 4th order finite difference stencils\n",
    "1. [Step 3](#rhss3d): Scalar Wave RHSs in Three Spatial Dimensions\n",
    "    1. [Step 3.a](#code_validation1): Code Validation against `ScalarWave.ScalarWave_RHSs` NRPy+ module\n",
    "    1. [Step 3.b](#ccode3d): C-code output example: Scalar wave RHSs with 10th order finite difference stencils and SIMD enabled\n",
    "1. [Step 4](#id): Setting up Initial Data for the Scalar Wave Equation\n",
    "    1. [Step 4.a](#planewave): The Monochromatic Plane-Wave Solution\n",
    "    1. [Step 4.b](#sphericalgaussian): The Spherical Gaussian Solution (*Courtesy Thiago Assumpção*)\n",
    "1. [Step 5](#code_validation2): Code Validation against `ScalarWave.InitialData` NRPy+ module\n",
    "1. [Step 6](#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 core NRPy+ modules \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "Let's start by importing all the needed modules from NRPy+:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:29.806437Z",
     "iopub.status.busy": "2021-03-07T17:22:29.805247Z",
     "iopub.status.idle": "2021-03-07T17:22:30.134159Z",
     "shell.execute_reply": "2021-03-07T17:22:30.133480Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step P1: Import needed NRPy+ core modules:\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import grid as gri               # NRPy+: Functions having to do with numerical grids\n",
    "import finite_difference as fin  # NRPy+: Finite difference C code generation module\n",
    "from outputC import lhrh         # NRPy+: Core C code output module\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rhss1d'></a>\n",
    "\n",
    "# Step 2: Scalar Wave RHSs in One Spatial Dimension \\[Back to [top](#toc)\\]\n",
    "$$\\label{rhss1d}$$\n",
    "\n",
    "To minimize complication, we will first restrict ourselves to solving the wave equation in one spatial dimension, so\n",
    "\\begin{align}\n",
    "\\partial_t u &= v \\\\\n",
    "\\partial_t v &= c^2 \\nabla^2 u \\\\\n",
    " &= c^2 \\partial_x^2 u.\n",
    "\\end{align}\n",
    "\n",
    "We will construct SymPy expressions of the right-hand sides of $u$ and $v$ using [NRPy+ finite-difference notation](Tutorial-Finite_Difference_Derivatives.ipynb) to represent the derivative, so that finite-difference C-code kernels can be easily constructed.\n",
    "\n",
    "Extension of this operator to higher spatial dimensions when using NRPy+ is straightforward, as we will see below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:30.144250Z",
     "iopub.status.busy": "2021-03-07T17:22:30.143242Z",
     "iopub.status.idle": "2021-03-07T17:22:30.455575Z",
     "shell.execute_reply": "2021-03-07T17:22:30.455020Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step P2: Define the C parameter wavespeed. The `wavespeed`\n",
    "#          variable is a proper SymPy variable, so it can be\n",
    "#          used in below expressions. In the C code, it acts\n",
    "#          just like a usual parameter, whose value is\n",
    "#          specified in the parameter file.\n",
    "thismodule = \"ScalarWave\"\n",
    "wavespeed = par.Cparameters(\"REAL\",thismodule,\"wavespeed\", 1.0)\n",
    "\n",
    "# Step 1: Set the spatial dimension parameter, and then read\n",
    "#         the parameter as DIM.\n",
    "par.set_parval_from_str(\"grid::DIM\", 1)\n",
    "DIM = par.parval_from_str(\"grid::DIM\")\n",
    "\n",
    "# Step 2: Register gridfunctions that are needed as input\n",
    "#         to the scalar wave RHS expressions.\n",
    "uu, vv = gri.register_gridfunctions(\"EVOL\", [\"uu\", \"vv\"])\n",
    "\n",
    "# Step 3: Declare the rank-2 indexed expression \\partial_{ij} u,\n",
    "#         which is symmetric about interchange of indices i and j\n",
    "#         Derivative variables like these must have an underscore\n",
    "#         in them, so the finite difference module can parse the\n",
    "#         variable name properly.\n",
    "uu_dDD = ixp.declarerank2(\"uu_dDD\", \"sym01\")\n",
    "\n",
    "# Step 4: Define right-hand sides for the evolution.\n",
    "uu_rhs = vv\n",
    "vv_rhs = 0\n",
    "for i in range(DIM):\n",
    "    vv_rhs += wavespeed*wavespeed*uu_dDD[i][i]\n",
    "\n",
    "vv_rhs = sp.simplify(vv_rhs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ccode1d'></a>\n",
    "\n",
    "## Step 2.a: C-code output example: Scalar wave RHSs with 4th order finite difference stencils \\[Back to [top](#toc)\\]\n",
    "$$\\label{ccode1d}$$\n",
    "\n",
    "As was discussed in [the finite difference section of the tutorial](Tutorial-Finite_Difference_Derivatives.ipynb), NRPy+ approximates derivatives using [finite-difference methods](https://en.wikipedia.org/wiki/Finite_difference),  the second-order derivative $\\partial_x^2$ accurate to fourth-order in uniform grid spacing $\\Delta x$ (from fitting the unique 4th-degree polynomial to 5 sample points of $u$) is given by\n",
    "\\begin{equation}\n",
    "\\left[\\partial_x^2 u(t,x)\\right]_j = \\frac{1}{(\\Delta x)^2}\n",
    "\\left(\n",
    "-\\frac{1}{12} \\left(u_{j+2} + u_{j-2}\\right) \n",
    "+ \\frac{4}{3}  \\left(u_{j+1} + u_{j-1}\\right)\n",
    "- \\frac{5}{2} u_j \\right)\n",
    "+ \\mathcal{O}\\left((\\Delta x)^4\\right).\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:30.144250Z",
     "iopub.status.busy": "2021-03-07T17:22:30.143242Z",
     "iopub.status.idle": "2021-03-07T17:22:30.455575Z",
     "shell.execute_reply": "2021-03-07T17:22:30.455020Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{\n",
      "   /*\n",
      "    * NRPy+ Finite Difference Code Generation, Step 1 of 2: Read from main memory and compute finite difference stencils:\n",
      "    */\n",
      "   /*\n",
      "    *  Original SymPy expression:\n",
      "    *  \"const double uu_dDD00 = invdx0**2*(-5*uu/2 + 4*uu_i0m1/3 - uu_i0m2/12 + 4*uu_i0p1/3 - uu_i0p2/12)\"\n",
      "    */\n",
      "   const double uu_i0m2 = in_gfs[IDX2(UUGF, i0-2)];\n",
      "   const double uu_i0m1 = in_gfs[IDX2(UUGF, i0-1)];\n",
      "   const double uu = in_gfs[IDX2(UUGF, i0)];\n",
      "   const double uu_i0p1 = in_gfs[IDX2(UUGF, i0+1)];\n",
      "   const double uu_i0p2 = in_gfs[IDX2(UUGF, i0+2)];\n",
      "   const double vv = in_gfs[IDX2(VVGF, i0)];\n",
      "   const double FDPart1_Rational_5_2 = 5.0/2.0;\n",
      "   const double FDPart1_Rational_1_12 = 1.0/12.0;\n",
      "   const double FDPart1_Rational_4_3 = 4.0/3.0;\n",
      "   const double uu_dDD00 = ((invdx0)*(invdx0))*(FDPart1_Rational_1_12*(-uu_i0m2 - uu_i0p2) + FDPart1_Rational_4_3*(uu_i0m1 + uu_i0p1) - FDPart1_Rational_5_2*uu);\n",
      "   /*\n",
      "    * NRPy+ Finite Difference Code Generation, Step 2 of 2: Evaluate SymPy expressions and write to main memory:\n",
      "    */\n",
      "   /*\n",
      "    *  Original SymPy expressions:\n",
      "    *  \"[rhs_gfs[IDX2(UUGF, i0)] = vv,\n",
      "    *    rhs_gfs[IDX2(VVGF, i0)] = uu_dDD00*wavespeed**2]\"\n",
      "    */\n",
      "   rhs_gfs[IDX2(UUGF, i0)] = vv;\n",
      "   rhs_gfs[IDX2(VVGF, i0)] = uu_dDD00*((wavespeed)*(wavespeed));\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "# Step 5: Set the finite differencing order to 4.\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", 4)\n",
    "\n",
    "# Step 6: Generate C code for scalarwave evolution equations,\n",
    "#         print output to the screen (standard out, or stdout).\n",
    "fin.FD_outputC(\"stdout\",\n",
    "               [lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"uu\"), rhs=uu_rhs),\n",
    "                lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"vv\"), rhs=vv_rhs)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Success!** Notice that indeed NRPy+ was able to compute the spatial derivative operator,\n",
    "\\begin{equation}\n",
    "\\left[\\partial_x^2 u(t,x)\\right]_j \\approx \\frac{1}{(\\Delta x)^2}\n",
    "\\left(\n",
    "-\\frac{1}{12} \\left(u_{j+2} + u_{j-2}\\right) \n",
    "+ \\frac{4}{3}  \\left(u_{j+1} + u_{j-1}\\right)\n",
    "- \\frac{5}{2} u_j \\right),\n",
    "\\end{equation}\n",
    "correctly (easier to read in the \"Original SymPy expressions\" comment block at the top of the C output.\n",
    "\n",
    "As NRPy+ is designed to generate codes in arbitrary coordinate systems, instead of sticking with Cartesian notation for 3D coordinates, $x,y,z$, we instead adopt $x_0,x_1,x_2$ for our coordinate labels. Thus you will notice the appearance of `invdx0`$=1/\\Delta x_0$, where $\\Delta x_0$ as the (uniform) grid spacing in the zeroth, or $x_0$ direction. In this case $x_0$ represents the $x$ direction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rhss3d'></a>\n",
    "\n",
    "# Step 3: Scalar Wave RHSs in Three Spatial Dimensions \\[Back to [top](#toc)\\]\n",
    "$$\\label{rhss3d}$$\n",
    "\n",
    "Let's next repeat the same process, only this time for the scalar wave equation in **3 spatial dimensions** (3D)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:30.530384Z",
     "iopub.status.busy": "2021-03-07T17:22:30.494206Z",
     "iopub.status.idle": "2021-03-07T17:22:30.807785Z",
     "shell.execute_reply": "2021-03-07T17:22:30.807120Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Define the C parameter wavespeed. The `wavespeed`\n",
    "#         variable is a proper SymPy variable, so it can be\n",
    "#         used in below expressions. In the C code, it acts\n",
    "#         just like a usual parameter, whose value is\n",
    "#         specified in the parameter file.\n",
    "wavespeed = par.Cparameters(\"REAL\", thismodule, \"wavespeed\", 1.0)\n",
    "\n",
    "# Step 2: Set the spatial dimension parameter\n",
    "#         to *FOUR* this time, and then read\n",
    "#         the parameter as DIM.\n",
    "par.set_parval_from_str(\"grid::DIM\", 3)\n",
    "DIM = par.parval_from_str(\"grid::DIM\")\n",
    "\n",
    "# Step 3a: Reset gridfunctions registered in 1D case above,\n",
    "#          to avoid NRPy+ throwing an error about double-\n",
    "#          registering gridfunctions, which is not allowed.\n",
    "gri.glb_gridfcs_list = []\n",
    "\n",
    "# Step 3b: Register gridfunctions that are needed as input\n",
    "#          to the scalar wave RHS expressions.\n",
    "uu, vv = gri.register_gridfunctions(\"EVOL\", [\"uu\", \"vv\"])\n",
    "\n",
    "# Step 4: Declare the rank-2 indexed expression \\partial_{ij} u,\n",
    "#         which is symmetric about interchange of indices i and j\n",
    "#         Derivative variables like these must have an underscore\n",
    "#         in them, so the finite difference module can parse the\n",
    "#         variable name properly.\n",
    "uu_dDD = ixp.declarerank2(\"uu_dDD\", \"sym01\")\n",
    "\n",
    "# Step 5: Define right-hand sides for the evolution.\n",
    "uu_rhs = vv\n",
    "vv_rhs = 0\n",
    "for i in range(DIM):\n",
    "    vv_rhs += wavespeed*wavespeed*uu_dDD[i][i]\n",
    "\n",
    "# Step 6: Simplify the expression for c^2 \\nabla^2 u (a.k.a., vv_rhs):\n",
    "vv_rhs = sp.simplify(vv_rhs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation1'></a>\n",
    "\n",
    "## Step 3.a:  Validate SymPy expressions against `ScalarWave.ScalarWave_RHSs` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation1}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for the RHSs of the three-spatial-dimension Scalar Wave equation (i.e., `uu_rhs` and `vv_rhs`) between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the [NRPy+ ScalarWave.ScalarWave_RHSs](../edit/ScalarWave/ScalarWave_RHSs.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:30.814167Z",
     "iopub.status.busy": "2021-03-07T17:22:30.813221Z",
     "iopub.status.idle": "2021-03-07T17:22:30.836178Z",
     "shell.execute_reply": "2021-03-07T17:22:30.835528Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between ScalarWave tutorial and NRPy+ module:\n",
      "uu_rhs - swrhs.uu_rhs = 0\t\t (should be zero)\n",
      "vv_rhs - swrhs.vv_rhs = 0\t\t (should be zero)\n"
     ]
    }
   ],
   "source": [
    "# Step 10: We already have SymPy expressions for uu_rhs and vv_rhs in\n",
    "#          terms of other SymPy variables. Even if we reset the list\n",
    "#          of NRPy+ gridfunctions, these *SymPy* expressions for\n",
    "#          uu_rhs and vv_rhs *will remain unaffected*.\n",
    "#\n",
    "#          Here, we will use the above-defined uu_rhs and vv_rhs to\n",
    "#          validate against the same expressions in the\n",
    "#          ScalarWave/ScalarWave_RHSs.py module,\n",
    "#          to ensure consistency between this tutorial\n",
    "#          (historically speaking, the tutorial was written first)\n",
    "#          and the ScalarWave_RHSs.py module itself.\n",
    "#\n",
    "# Reset the list of gridfunctions, as registering a gridfunction\n",
    "#   twice will spawn an error.\n",
    "gri.glb_gridfcs_list = []\n",
    "\n",
    "# Step 11: Call the ScalarWave_RHSs() function from within the\n",
    "#         ScalarWave/ScalarWave_RHSs.py module,\n",
    "#         which should do exactly the same as in Steps 1-10 above.\n",
    "import ScalarWave.ScalarWave_RHSs as swrhs\n",
    "swrhs.ScalarWave_RHSs()\n",
    "\n",
    "# Step 12: Consistency check between the tutorial notebook above\n",
    "#         and the ScalarWave_RHSs() function from within the\n",
    "#         ScalarWave/ScalarWave_RHSs.py module.\n",
    "print(\"Consistency check between ScalarWave tutorial and NRPy+ module:\")\n",
    "print(\"uu_rhs - swrhs.uu_rhs = \"+str(sp.simplify(uu_rhs - swrhs.uu_rhs))+\"\\t\\t (should be zero)\")\n",
    "print(\"vv_rhs - swrhs.vv_rhs = \"+str(sp.simplify(vv_rhs - swrhs.vv_rhs))+\"\\t\\t (should be zero)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='ccode3d'></a>\n",
    "\n",
    "## Step 3.b: C-code output example: Scalar wave RHSs with 10th order finite difference stencils and SIMD enabled \\[Back to [top](#toc)\\]\n",
    "$$\\label{ccode3d}$$\n",
    "\n",
    "Next we'll output the above expressions as Ccode, using the [NRPy+ finite-differencing C code kernel generation infrastructure](Tutorial-Finite_Difference_Derivatives.ipynb). This code will represent spatial derivatives as 10th-order finite differences and output the C code with [SIMD](https://en.wikipedia.org/wiki/SIMD) enabled. ([Common-subexpression elimination](https://en.wikipedia.org/wiki/Common_subexpression_elimination) is enabled by default.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:30.530384Z",
     "iopub.status.busy": "2021-03-07T17:22:30.494206Z",
     "iopub.status.idle": "2021-03-07T17:22:30.807785Z",
     "shell.execute_reply": "2021-03-07T17:22:30.807120Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{\n",
      "   /*\n",
      "    * NRPy+ Finite Difference Code Generation, Step 1 of 2: Read from main memory and compute finite difference stencils:\n",
      "    */\n",
      "   /*\n",
      "    *  Original SymPy expressions:\n",
      "    *  \"[const REAL_SIMD_ARRAY uu_dDD00 = invdx0**2*(-5269*uu/1800 + 5*uu_i0m1_i1_i2/3 - 5*uu_i0m2_i1_i2/21 + 5*uu_i0m3_i1_i2/126 - 5*uu_i0m4_i1_i2/1008 + uu_i0m5_i1_i2/3150 + 5*uu_i0p1_i1_i2/3 - 5*uu_i0p2_i1_i2/21 + 5*uu_i0p3_i1_i2/126 - 5*uu_i0p4_i1_i2/1008 + uu_i0p5_i1_i2/3150),\n",
      "    *    const REAL_SIMD_ARRAY uu_dDD11 = invdx1**2*(-5269*uu/1800 + 5*uu_i0_i1m1_i2/3 - 5*uu_i0_i1m2_i2/21 + 5*uu_i0_i1m3_i2/126 - 5*uu_i0_i1m4_i2/1008 + uu_i0_i1m5_i2/3150 + 5*uu_i0_i1p1_i2/3 - 5*uu_i0_i1p2_i2/21 + 5*uu_i0_i1p3_i2/126 - 5*uu_i0_i1p4_i2/1008 + uu_i0_i1p5_i2/3150),\n",
      "    *    const REAL_SIMD_ARRAY uu_dDD22 = invdx2**2*(-5269*uu/1800 + 5*uu_i0_i1_i2m1/3 - 5*uu_i0_i1_i2m2/21 + 5*uu_i0_i1_i2m3/126 - 5*uu_i0_i1_i2m4/1008 + uu_i0_i1_i2m5/3150 + 5*uu_i0_i1_i2p1/3 - 5*uu_i0_i1_i2p2/21 + 5*uu_i0_i1_i2p3/126 - 5*uu_i0_i1_i2p4/1008 + uu_i0_i1_i2p5/3150)]\"\n",
      "    */\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1_i2m5 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2-5)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1_i2m4 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2-4)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1_i2m3 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2-3)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1_i2m2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2-2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1_i2m1 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2-1)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1m5_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1-5,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1m4_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1-4,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1m3_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1-3,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1m2_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1-2,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1m1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1-1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0m5_i1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0-5,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0m4_i1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0-4,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0m3_i1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0-3,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0m2_i1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0-2,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0m1_i1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0-1,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0p1_i1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0+1,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0p2_i1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0+2,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0p3_i1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0+3,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0p4_i1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0+4,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0p5_i1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0+5,i1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1p1_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1+1,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1p2_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1+2,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1p3_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1+3,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1p4_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1+4,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1p5_i2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1+5,i2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1_i2p1 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2+1)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1_i2p2 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2+2)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1_i2p3 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2+3)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1_i2p4 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2+4)]);\n",
      "   const REAL_SIMD_ARRAY uu_i0_i1_i2p5 = ReadSIMD(&in_gfs[IDX4(UUGF, i0,i1,i2+5)]);\n",
      "   const REAL_SIMD_ARRAY vv = ReadSIMD(&in_gfs[IDX4(VVGF, i0,i1,i2)]);\n",
      "   const double tmpFDPart1_NegativeOne_ = -1.0;\n",
      "   const REAL_SIMD_ARRAY FDPart1_NegativeOne_ = ConstSIMD(tmpFDPart1_NegativeOne_);\n",
      "   \n",
      "   const double tmpFDPart1_Rational_1_3150 = 1.0/3150.0;\n",
      "   const REAL_SIMD_ARRAY FDPart1_Rational_1_3150 = ConstSIMD(tmpFDPart1_Rational_1_3150);\n",
      "   \n",
      "   const double tmpFDPart1_Rational_5269_1800 = 5269.0/1800.0;\n",
      "   const REAL_SIMD_ARRAY FDPart1_Rational_5269_1800 = ConstSIMD(tmpFDPart1_Rational_5269_1800);\n",
      "   \n",
      "   const double tmpFDPart1_Rational_5_1008 = 5.0/1008.0;\n",
      "   const REAL_SIMD_ARRAY FDPart1_Rational_5_1008 = ConstSIMD(tmpFDPart1_Rational_5_1008);\n",
      "   \n",
      "   const double tmpFDPart1_Rational_5_126 = 5.0/126.0;\n",
      "   const REAL_SIMD_ARRAY FDPart1_Rational_5_126 = ConstSIMD(tmpFDPart1_Rational_5_126);\n",
      "   \n",
      "   const double tmpFDPart1_Rational_5_21 = 5.0/21.0;\n",
      "   const REAL_SIMD_ARRAY FDPart1_Rational_5_21 = ConstSIMD(tmpFDPart1_Rational_5_21);\n",
      "   \n",
      "   const double tmpFDPart1_Rational_5_3 = 5.0/3.0;\n",
      "   const REAL_SIMD_ARRAY FDPart1_Rational_5_3 = ConstSIMD(tmpFDPart1_Rational_5_3);\n",
      "   \n",
      "   const REAL_SIMD_ARRAY FDPart1_0 = MulSIMD(FDPart1_Rational_5269_1800, uu);\n",
      "   const REAL_SIMD_ARRAY uu_dDD00 = MulSIMD(MulSIMD(invdx0, invdx0), FusedMulAddSIMD(FDPart1_Rational_5_126, AddSIMD(uu_i0m3_i1_i2, uu_i0p3_i1_i2), FusedMulAddSIMD(FDPart1_Rational_5_3, AddSIMD(uu_i0m1_i1_i2, uu_i0p1_i1_i2), FusedMulSubSIMD(FDPart1_Rational_1_3150, AddSIMD(uu_i0m5_i1_i2, uu_i0p5_i1_i2), FusedMulAddSIMD(FDPart1_Rational_5_1008, AddSIMD(uu_i0m4_i1_i2, uu_i0p4_i1_i2), FusedMulAddSIMD(FDPart1_Rational_5_21, AddSIMD(uu_i0m2_i1_i2, uu_i0p2_i1_i2), FDPart1_0))))));\n",
      "   const REAL_SIMD_ARRAY uu_dDD11 = MulSIMD(MulSIMD(invdx1, invdx1), FusedMulAddSIMD(FDPart1_Rational_5_126, AddSIMD(uu_i0_i1m3_i2, uu_i0_i1p3_i2), FusedMulAddSIMD(FDPart1_Rational_5_3, AddSIMD(uu_i0_i1m1_i2, uu_i0_i1p1_i2), FusedMulSubSIMD(FDPart1_Rational_1_3150, AddSIMD(uu_i0_i1m5_i2, uu_i0_i1p5_i2), FusedMulAddSIMD(FDPart1_Rational_5_1008, AddSIMD(uu_i0_i1m4_i2, uu_i0_i1p4_i2), FusedMulAddSIMD(FDPart1_Rational_5_21, AddSIMD(uu_i0_i1m2_i2, uu_i0_i1p2_i2), FDPart1_0))))));\n",
      "   const REAL_SIMD_ARRAY uu_dDD22 = MulSIMD(MulSIMD(invdx2, invdx2), FusedMulAddSIMD(FDPart1_Rational_5_126, AddSIMD(uu_i0_i1_i2m3, uu_i0_i1_i2p3), FusedMulAddSIMD(FDPart1_Rational_5_3, AddSIMD(uu_i0_i1_i2m1, uu_i0_i1_i2p1), FusedMulSubSIMD(FDPart1_Rational_1_3150, AddSIMD(uu_i0_i1_i2m5, uu_i0_i1_i2p5), FusedMulAddSIMD(FDPart1_Rational_5_1008, AddSIMD(uu_i0_i1_i2m4, uu_i0_i1_i2p4), FusedMulAddSIMD(FDPart1_Rational_5_21, AddSIMD(uu_i0_i1_i2m2, uu_i0_i1_i2p2), FDPart1_0))))));\n",
      "   /*\n",
      "    * NRPy+ Finite Difference Code Generation, Step 2 of 2: Evaluate SymPy expressions and write to main memory:\n",
      "    */\n",
      "   /*\n",
      "    *  Original SymPy expressions:\n",
      "    *  \"[const REAL_SIMD_ARRAY __RHS_exp_0 = vv,\n",
      "    *    const REAL_SIMD_ARRAY __RHS_exp_1 = wavespeed**2*(uu_dDD00 + uu_dDD11 + uu_dDD22)]\"\n",
      "    */\n",
      "   const REAL_SIMD_ARRAY __RHS_exp_0 = vv;\n",
      "   const REAL_SIMD_ARRAY __RHS_exp_1 = MulSIMD(MulSIMD(wavespeed, wavespeed), AddSIMD(uu_dDD00, AddSIMD(uu_dDD11, uu_dDD22)));\n",
      "   WriteSIMD(&rhs_gfs[IDX4(UUGF, i0, i1, i2)], __RHS_exp_0);\n",
      "   WriteSIMD(&rhs_gfs[IDX4(VVGF, i0, i1, i2)], __RHS_exp_1);\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "# Step 7: Set the finite differencing order to 10.\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", 10)\n",
    "\n",
    "# Step 8: Generate C code for scalarwave evolution equations,\n",
    "#         print output to the screen (standard out, or stdout).\n",
    "fin.FD_outputC(\"stdout\",\n",
    "               [lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"uu\"),rhs=uu_rhs),\n",
    "                lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"vv\"),rhs=vv_rhs)], params=\"enable_SIMD=True\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='id'></a>\n",
    "\n",
    "# Step 4: Setting up Initial Data for the Scalar Wave Equation \\[Back to [top](#toc)\\]\n",
    "$$\\label{id}$$\n",
    "\n",
    "<a id='planewave'></a>\n",
    "\n",
    "## Step 4.a: The Monochromatic Plane-Wave Solution \\[Back to [top](#toc)\\]\n",
    "$$\\label{planewave}$$\n",
    "\n",
    "The solution to the scalar wave equation for a monochromatic (single-wavelength) wave traveling in the $\\hat{k}$ direction is\n",
    "$$u(\\vec{x},t) = f(\\hat{k}\\cdot\\vec{x} - c t),$$\n",
    "where $\\hat{k}$ is a unit vector. We choose $f(\\hat{k}\\cdot\\vec{x} - c t)$ to take the form\n",
    "$$\n",
    "f(\\hat{k}\\cdot\\vec{x} - c t) = \\sin\\left(\\hat{k}\\cdot\\vec{x} - c t\\right) + 2,\n",
    "$$\n",
    "where we add the $+2$ to ensure that the exact solution never crosses through zero. In places where the exact solution passes through zero, the relative error (i.e., the most common error measure used to check that the numerical solution converges to the exact solution) is undefined. Also, $f(\\hat{k}\\cdot\\vec{x} - c t)$ plus a constant is still a solution to the wave equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:30.861959Z",
     "iopub.status.busy": "2021-03-07T17:22:30.861025Z",
     "iopub.status.idle": "2021-03-07T17:22:30.864015Z",
     "shell.execute_reply": "2021-03-07T17:22:30.863446Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Set parameters defined in other modules\n",
    "xx = gri.xx # Sets the Cartesian coordinates xx[0]=x; xx[1]=y; xx[2]=z\n",
    "\n",
    "# Step 2: Declare free parameters intrinsic to these initial data\n",
    "time = par.Cparameters(\"REAL\", thismodule, \"time\",0.0)\n",
    "kk = par.Cparameters(\"REAL\", thismodule, [\"kk0\", \"kk1\", \"kk2\"],[1.0,1.0,1.0])\n",
    "\n",
    "# Step 3: Normalize the k vector\n",
    "kk_norm = sp.sqrt(kk[0]**2 + kk[1]**2 + kk[2]**2)\n",
    "\n",
    "# Step 4: Compute k.x\n",
    "dot_product = sp.sympify(0)\n",
    "for i in range(DIM):\n",
    "    dot_product += xx[i]*kk[i]\n",
    "dot_product /= kk_norm\n",
    "\n",
    "# Step 5: Set initial data for uu and vv, where vv_ID = \\partial_t uu_ID.\n",
    "uu_ID_PlaneWave = sp.sin(dot_product - wavespeed*time)+2\n",
    "vv_ID_PlaneWave = sp.diff(uu_ID_PlaneWave, time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we verify that $f(\\hat{k}\\cdot\\vec{x} - c t)$ satisfies the wave equation, by computing\n",
    "$$\\left(c^2 \\nabla^2 - \\partial_t^2 \\right)\\ f\\left(\\hat{k}\\cdot\\vec{x} - c t\\right),$$\n",
    "and confirming the result is exactly zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:30.938389Z",
     "iopub.status.busy": "2021-03-07T17:22:30.902532Z",
     "iopub.status.idle": "2021-03-07T17:22:31.136590Z",
     "shell.execute_reply": "2021-03-07T17:22:31.137206Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(wavespeed**2*(sp.diff(uu_ID_PlaneWave,xx[0],2) +\n",
    "                          sp.diff(uu_ID_PlaneWave,xx[1],2) +\n",
    "                          sp.diff(uu_ID_PlaneWave,xx[2],2))\n",
    "            - sp.diff(uu_ID_PlaneWave,time,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sphericalgaussian'></a>\n",
    "\n",
    "## Step 4.b: The Spherical Gaussian Solution \\[Back to [top](#toc)\\]\n",
    "$$\\label{sphericalgaussian}$$\n",
    "\n",
    "Here we will implement the spherical Gaussian solution, consists of ingoing and outgoing wave fronts:\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",
    "where $c$ is the wavespeed, and $\\sigma$ is the width of the Gaussian (i.e., the \"standard deviation\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:31.181373Z",
     "iopub.status.busy": "2021-03-07T17:22:31.150833Z",
     "iopub.status.idle": "2021-03-07T17:22:31.184137Z",
     "shell.execute_reply": "2021-03-07T17:22:31.183590Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 1: Set parameters defined in other modules\n",
    "xx = gri.xx # Sets the Cartesian coordinates xx[0]=x; xx[1]=y; xx[2]=z\n",
    "\n",
    "# Step 2: Declare free parameters intrinsic to these initial data\n",
    "time  = par.Cparameters(\"REAL\", thismodule, \"time\",0.0)\n",
    "sigma = par.Cparameters(\"REAL\", thismodule, \"sigma\",3.0)\n",
    "\n",
    "# Step 4: Compute r\n",
    "r = sp.sympify(0)\n",
    "for i in range(DIM):\n",
    "    r += xx[i]**2\n",
    "r = sp.sqrt(r)\n",
    "\n",
    "# Step 5: Set initial data for uu and vv, where vv_ID = \\partial_t uu_ID.\n",
    "uu_ID_SphericalGaussianOUT = +(r - wavespeed*time)/r * sp.exp( -(r - wavespeed*time)**2 / (2*sigma**2) )\n",
    "uu_ID_SphericalGaussianIN  = +(r + wavespeed*time)/r * sp.exp( -(r + wavespeed*time)**2 / (2*sigma**2) )\n",
    "uu_ID_SphericalGaussian = uu_ID_SphericalGaussianOUT + uu_ID_SphericalGaussianIN + sp.sympify(1)\n",
    "vv_ID_SphericalGaussian = sp.diff(uu_ID_SphericalGaussian, time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the wave equation is linear, both the leftgoing and rightgoing waves must satisfy the wave equation, which implies that their sum also satisfies the wave equation. \n",
    "\n",
    "Next we verify that $u(r,t)$ satisfies the wave equation, by computing\n",
    "$$\\left(c^2 \\nabla^2 - \\partial_t^2 \\right)\\left\\{u_{\\rm R}(r,t)\\right\\},$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\\left(c^2 \\nabla^2 - \\partial_t^2 \\right)\\left\\{u_{\\rm L}(r,t)\\right\\},$$\n",
    "\n",
    "are separately zero. We do this because SymPy has difficulty simplifying the combined expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:31.259089Z",
     "iopub.status.busy": "2021-03-07T17:22:31.223170Z",
     "iopub.status.idle": "2021-03-07T17:22:42.951534Z",
     "shell.execute_reply": "2021-03-07T17:22:42.950760Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "print(sp.simplify(wavespeed**2*(sp.diff(uu_ID_SphericalGaussianOUT,xx[0],2) +\n",
    "                                sp.diff(uu_ID_SphericalGaussianOUT,xx[1],2) +\n",
    "                                sp.diff(uu_ID_SphericalGaussianOUT,xx[2],2))\n",
    "                - sp.diff(uu_ID_SphericalGaussianOUT,time,2)) )\n",
    "\n",
    "print(sp.simplify(wavespeed**2*(sp.diff(uu_ID_SphericalGaussianIN,xx[0],2) +\n",
    "                                sp.diff(uu_ID_SphericalGaussianIN,xx[1],2) +\n",
    "                                sp.diff(uu_ID_SphericalGaussianIN,xx[2],2))\n",
    "                    - sp.diff(uu_ID_SphericalGaussianIN,time,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation2'></a>\n",
    "\n",
    "# Step 5: Code Validation against `ScalarWave.InitialData` NRPy+ module  \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation2}$$\n",
    "\n",
    "As a code validation check, we will verify agreement in the SymPy expressions for plane-wave initial data for the Scalar Wave equation between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [ScalarWave.InitialData](../edit/ScalarWave/InitialData.py) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:42.961919Z",
     "iopub.status.busy": "2021-03-07T17:22:42.959789Z",
     "iopub.status.idle": "2021-03-07T17:22:43.007222Z",
     "shell.execute_reply": "2021-03-07T17:22:43.006498Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between ScalarWave tutorial and NRPy+ module: PlaneWave Case\n",
      "TESTS PASSED!\n",
      "Consistency check between ScalarWave tutorial and NRPy+ module: SphericalGaussian Case\n",
      "TESTS PASSED!\n"
     ]
    }
   ],
   "source": [
    "# We just defined SymPy expressions for uu_ID and vv_ID in\n",
    "# terms of other SymPy variables. Here, we will use the\n",
    "# above-defined uu_ID and vv_ID to validate against the\n",
    "# same expressions in the ScalarWave/InitialData.py\n",
    "# module, to ensure consistency between this tutorial\n",
    "# (historically speaking, the tutorial was written first)\n",
    "# and the PlaneWave ID module itself.\n",
    "#\n",
    "# Step 6: Call the InitialData(Type=\"PlaneWave\") function from within the\n",
    "#         ScalarWave/InitialData.py module,\n",
    "#         which should do exactly the same as in Steps 1-5 above.\n",
    "import ScalarWave.InitialData as swid\n",
    "swid.InitialData(Type=\"PlaneWave\")\n",
    "\n",
    "# Step 7: Consistency check between the tutorial notebook above\n",
    "#         and the PlaneWave option from within the\n",
    "#         ScalarWave/InitialData.py module.\n",
    "print(\"Consistency check between ScalarWave tutorial and NRPy+ module: PlaneWave Case\")\n",
    "if sp.simplify(uu_ID_PlaneWave - swid.uu_ID) != 0:\n",
    "    print(\"TEST FAILED: uu_ID_PlaneWave - swid.uu_ID = \"+str(sp.simplify(uu_ID_PlaneWave - swid.uu_ID))+\"\\t\\t (should be zero)\")\n",
    "    sys.exit(1)\n",
    "if sp.simplify(vv_ID_PlaneWave - swid.vv_ID) != 0:\n",
    "    print(\"TEST FAILED: vv_ID_PlaneWave - swid.vv_ID = \"+str(sp.simplify(vv_ID_PlaneWave - swid.vv_ID))+\"\\t\\t (should be zero)\")\n",
    "    sys.exit(1)\n",
    "print(\"TESTS PASSED!\")\n",
    "\n",
    "\n",
    "# Step 8: Consistency check between the tutorial notebook above\n",
    "#         and the SphericalGaussian option from within the\n",
    "#         ScalarWave/InitialData.py module.\n",
    "swid.InitialData(Type=\"SphericalGaussian\")\n",
    "print(\"Consistency check between ScalarWave tutorial and NRPy+ module: SphericalGaussian Case\")\n",
    "if sp.simplify(uu_ID_SphericalGaussian - swid.uu_ID) != 0:\n",
    "    print(\"TEST FAILED: uu_ID_SphericalGaussian - swid.uu_ID = \"+str(sp.simplify(uu_ID_SphericalGaussian - swid.uu_ID))+\"\\t\\t (should be zero)\")\n",
    "    sys.exit(1)\n",
    "if sp.simplify(vv_ID_SphericalGaussian - swid.vv_ID) != 0:\n",
    "    print(\"TEST FAILED: vv_ID_SphericalGaussian - swid.vv_ID = \"+str(sp.simplify(vv_ID_SphericalGaussian - swid.vv_ID))+\"\\t\\t (should be zero)\")\n",
    "    sys.exit(1)\n",
    "print(\"TESTS PASSED!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 6: 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-ScalarWave.pdf](Tutorial-ScalarWave.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": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:43.011416Z",
     "iopub.status.busy": "2021-03-07T17:22:43.010796Z",
     "iopub.status.idle": "2021-03-07T17:22:46.513390Z",
     "shell.execute_reply": "2021-03-07T17:22:46.514038Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ScalarWave.tex, and compiled LaTeX file to PDF file\n",
      "    Tutorial-ScalarWave.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ScalarWave\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.0rc2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}