Tutorial-ScalarWaveCurvilinear.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 right-hand-side of the scalar wave equation, in ***curvilinear*** coordinates, using a reference metric formalism\n",
    "\n",
    "## Author: Zach Etienne\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). 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: [ScalarWaveCurvilinear/ScalarWaveCurvilinear_RHSs.py](../edit/ScalarWaveCurvilinear/ScalarWaveCurvilinear_RHSs.py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "0. [Preliminaries](#prelim): Reference Metrics and Picking Best Coordinate System to Solve the PDE\n",
    "1. [Example](#example): The scalar wave equation in spherical coordinates\n",
    "1. [Step 1](#contracted_christoffel): Contracted Christoffel symbols $\\hat{\\Gamma}^i = \\hat{g}^{ij}\\hat{\\Gamma}^k_{ij}$ in spherical coordinates, using NRPy+\n",
    "1. [Step 2](#rhs_scalarwave_spherical): The right-hand side of the scalar wave equation in spherical coordinates, using NRPy+\n",
    "1. [Step 3](#code_validation): Code Validation against `ScalarWave.ScalarWaveCurvilinear_RHSs`  NRPy+ Module\n",
    "1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='prelim'></a>\n",
    "\n",
    "# Preliminaries: Reference Metrics and Picking Best Coordinate System to Solve the PDE \\[Back to [top](#toc)\\]\n",
    "$$\\label{prelim}$$\n",
    "\n",
    "Recall from [NRPy+ tutorial notebook on the Cartesian scalar wave equation](Tutorial-ScalarWave.ipynb), the scalar wave equation in 3D Cartesian coordinates is given by\n",
    "\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 Cartesian coordinates in space: $u = u(t,x,y,z)$ (spatial dimension as-yet unspecified), and subject to some initial condition\n",
    "$$u(0,x,y,z) = f(x,y,z),$$\n",
    "\n",
    "with suitable (sometimes approximate) spatial boundary conditions.\n",
    "\n",
    "To simplify this equation, let's first choose units such that $c=1$. Alternative wave speeds can be constructed\n",
    "by simply rescaling the time coordinate, with the net effect being that the time $t$ is replaced with time in dimensions of space; i.e., $t\\to c t$:\n",
    "\n",
    "$$\\partial_t^2 u = \\nabla^2 u.$$\n",
    "\n",
    "As we learned in the [NRPy+ tutorial notebook on reference metrics](Tutorial-Reference_Metric.ipynb), reference metrics are a means to pick the best coordinate system for the PDE we wish to solve. However, to take advantage of reference metrics requires first that we generalize the PDE. In the case of the scalar wave equation, this involves first rewriting in [Einstein notation](https://en.wikipedia.org/wiki/Einstein_notation) (with implied summation over repeated indices) via\n",
    "\n",
    "$$(-\\partial_t^2 + \\nabla^2) u = \\eta^{\\mu\\nu} u_{,\\ \\mu\\nu} = 0,$$\n",
    "\n",
    "where $u_{,\\mu\\nu} = \\partial_\\mu \\partial_\\nu u$, and $\\eta^{\\mu\\nu}$ is the contravariant flat-space metric tensor with components $\\text{diag}(-1,1,1,1)$.\n",
    "\n",
    "Next we apply the \"comma-goes-to-semicolon rule\" and replace $\\eta^{\\mu\\nu}$ with $\\hat{g}^{\\mu\\nu}$ to generalize the scalar wave equation to an arbitrary reference metric $\\hat{g}^{\\mu\\nu}$:\n",
    "\n",
    "$$\\hat{g}^{\\mu\\nu} u_{;\\ \\mu\\nu} = \\hat{g}^{\\mu\\nu} \\hat{\\nabla}_{\\mu} \\hat{\\nabla}_{\\nu} u = 0,$$\n",
    "\n",
    "where $\\hat{\\nabla}_{\\mu}$ denotes the [covariant derivative](https://en.wikipedia.org/wiki/Covariant_derivative) with respect to the reference metric basis vectors $\\hat{x}^{\\mu}$, and $\\hat{g}^{\\mu \\nu} \\hat{\\nabla}_{\\mu} \\hat{\\nabla}_{\\nu} u$ is the covariant\n",
    "[D'Alembertian](https://en.wikipedia.org/wiki/D%27Alembert_operator) of $u$.\n",
    "\n",
    "For example, suppose we wish to model a short-wavelength wave that is nearly spherical. In this case, if we were to solve the wave equation PDE in Cartesian coordinates, we would in principle need high resolution in all three cardinal directions. If instead we chose spherical coordinates centered at the center of the wave, we might need high resolution only in the radial direction, with only a few points required in the angular directions. Thus choosing spherical coordinates would be far more computationally efficient than modeling the wave in Cartesian coordinates.\n",
    "\n",
    "Let's now expand the covariant scalar wave equation in arbitrary coordinates. Since the covariant derivative of a scalar is equivalent to its partial derivative, we have\n",
    "\\begin{align}\n",
    "0 &= \\hat{g}^{\\mu \\nu} \\hat{\\nabla}_{\\mu} \\hat{\\nabla}_{\\nu} u \\\\\n",
    "&= \\hat{g}^{\\mu \\nu} \\hat{\\nabla}_{\\mu} \\partial_{\\nu} u.\n",
    "\\end{align}\n",
    "\n",
    "$\\partial_{\\nu} u$ transforms as a one-form under covariant differentiation, so we have\n",
    "$$\\hat{\\nabla}_{\\mu} \\partial_{\\nu} u = \\partial_{\\mu} \\partial_{\\nu} u - \\hat{\\Gamma}^\\tau_{\\mu\\nu} \\partial_\\tau u,$$\n",
    "where \n",
    "\n",
    "$$\\hat{\\Gamma}^\\tau_{\\mu\\nu} = \\frac{1}{2} \\hat{g}^{\\tau\\alpha} \\left(\\partial_\\nu \\hat{g}_{\\alpha\\mu} + \\partial_\\mu \\hat{g}_{\\alpha\\nu} - \\partial_\\alpha \\hat{g}_{\\mu\\nu} \\right)$$\n",
    "are the [Christoffel symbols](https://en.wikipedia.org/wiki/Christoffel_symbols) associated with the reference metric $\\hat{g}_{\\mu\\nu}$.\n",
    "\n",
    "Then the scalar wave equation is written:\n",
    "$$0 = \\hat{g}^{\\mu \\nu} \\left( \\partial_{\\mu} \\partial_{\\nu} u - \\hat{\\Gamma}^\\tau_{\\mu\\nu} \\partial_\\tau u\\right).$$\n",
    "\n",
    "Define the contracted Christoffel symbols:\n",
    "$$\\hat{\\Gamma}^\\tau = \\hat{g}^{\\mu\\nu} \\hat{\\Gamma}^\\tau_{\\mu\\nu}.$$\n",
    "\n",
    "Then the scalar wave equation is given by\n",
    "$$0 = \\hat{g}^{\\mu \\nu} \\partial_{\\mu} \\partial_{\\nu} u - \\hat{\\Gamma}^\\tau \\partial_\\tau u.$$\n",
    "\n",
    "The reference metrics we adopt satisfy\n",
    "$$\\hat{g}^{t \\nu} = -\\delta^{t \\nu},$$\n",
    "where $\\delta^{t \\nu}$ is the [Kronecker delta](https://en.wikipedia.org/wiki/Kronecker_delta). Therefore the scalar wave equation in curvilinear coordinates can be written\n",
    "\\begin{align}\n",
    "0 &= \\hat{g}^{\\mu \\nu} \\partial_{\\mu} \\partial_{\\nu} u - \\hat{\\Gamma}^\\tau \\partial_\\tau u \\\\\n",
    "&= -\\partial_t^2 u + \\hat{g}^{i j} \\partial_{i} \\partial_{j} u - \\hat{\\Gamma}^i \\partial_i u \\\\\n",
    "\\implies \\partial_t^2 u &= \\hat{g}^{i j} \\partial_{i} \\partial_{j} u - \\hat{\\Gamma}^i \\partial_i u,\n",
    "\\end{align}\n",
    "where repeated Latin indices denote implied summation over *spatial* components only. This module implements the bottom equation for arbitrary reference metrics satisfying $\\hat{g}^{t \\nu} = -\\delta^{t \\nu}$. To gain an appreciation for what NRPy+ accomplishes automatically, let's first work out the scalar wave equation in spherical coordinates by hand:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='example'></a>\n",
    "\n",
    "# Example: The scalar wave equation in spherical coordinates \\[Back to [top](#toc)\\]\n",
    "$$\\label{example}$$\n",
    "\n",
    "For example, the spherical reference metric is written\n",
    "\n",
    "$$\\hat{g}_{\\mu\\nu} = \\begin{pmatrix}\n",
    "-1 & 0 & 0 & 0 \\\\\n",
    " 0 & 1 & 0 & 0 \\\\\n",
    " 0 & 0 & r^2 & 0 \\\\\n",
    " 0 & 0 & 0 & r^2 \\sin^2 \\theta \\\\\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "Since the inverse of a diagonal matrix is simply the inverse of the diagonal elements, we can write \n",
    "$$\\hat{g}^{\\mu\\nu} = \\begin{pmatrix}\n",
    "-1 & 0 & 0 & 0 \\\\\n",
    " 0 & 1 & 0 & 0 \\\\\n",
    " 0 & 0 & \\frac{1}{r^2} & 0 \\\\\n",
    " 0 & 0 & 0 & \\frac{1}{r^2 \\sin^2 \\theta} \\\\\n",
    "\\end{pmatrix}.$$\n",
    "\n",
    "The scalar wave equation in these coordinates can thus be written\n",
    "\\begin{align}\n",
    "0 &= \\hat{g}^{\\mu \\nu} \\partial_{\\mu} \\partial_{\\nu} u - \\hat{\\Gamma}^\\tau \\partial_\\tau u \\\\\n",
    "&= \\hat{g}^{tt} \\partial_t^2 u + \\hat{g}^{rr} \\partial_r^2 u + \\hat{g}^{\\theta\\theta} \\partial_\\theta^2 u  + \\hat{g}^{\\phi\\phi} \\partial_\\phi^2 u - \\hat{\\Gamma}^\\tau \\partial_\\tau u \\\\\n",
    "&= -\\partial_t^2 u + \\partial_r^2 u + \\frac{1}{r^2} \\partial_\\theta^2\n",
    "u + \\frac{1}{r^2 \\sin^2 \\theta} \\partial_\\phi^2 u - \\hat{\\Gamma}^\\tau \\partial_\\tau u\\\\\n",
    "\\implies \\partial_t^2 u &= \\partial_r^2 u + \\frac{1}{r^2} \\partial_\\theta^2\n",
    "u + \\frac{1}{r^2 \\sin^2 \\theta} \\partial_\\phi^2 u - \\hat{\\Gamma}^\\tau \\partial_\\tau u\n",
    "\\end{align}\n",
    "\n",
    "The contracted Christoffel symbols \n",
    "$\\hat{\\Gamma}^\\tau$ can then be computed directly from the metric $\\hat{g}_{\\mu\\nu}$.\n",
    "\n",
    "It can be shown (exercise to the reader) that the only nonzero\n",
    "components of $\\hat{\\Gamma}^\\tau$ in static spherical polar coordinates are\n",
    "given by\n",
    "\\begin{align}\n",
    "\\hat{\\Gamma}^r &= -\\frac{2}{r} \\\\\n",
    "\\hat{\\Gamma}^\\theta &= -\\frac{\\cos\\theta}{r^2 \\sin\\theta}.\n",
    "\\end{align}\n",
    "\n",
    "Thus we have found the Laplacian in spherical coordinates is simply:\n",
    "\n",
    "\\begin{align}\n",
    "\\nabla^2 u &= \n",
    "\\partial_r^2 u + \\frac{1}{r^2} \\partial_\\theta^2 u + \\frac{1}{r^2 \\sin^2 \\theta} \\partial_\\phi^2 u - \\hat{\\Gamma}^\\tau \\partial_\\tau u\\\\\n",
    "&= \\partial_r^2 u + \\frac{1}{r^2} \\partial_\\theta^2 u + \\frac{1}{r^2 \\sin^2 \\theta} \\partial_\\phi^2 u +  \\frac{2}{r} \\partial_r u + \\frac{\\cos\\theta}{r^2 \\sin\\theta} \\partial_\\theta u\n",
    "\\end{align}\n",
    "(cf. http://mathworld.wolfram.com/SphericalCoordinates.html; though note that they defined the angle $\\phi$ as $\\theta$ and $\\theta$ as $\\phi$.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='contracted_christoffel'></a>\n",
    "\n",
    "# Step 1: Contracted Christoffel symbols $\\hat{\\Gamma}^i = \\hat{g}^{ij}\\hat{\\Gamma}^k_{ij}$ in spherical coordinates, using NRPy+ \\[Back to [top](#toc)\\]\n",
    "$$\\label{contracted_christoffel}$$\n",
    "\n",
    "Let's next use NRPy+ to derive the contracted Christoffel symbols\n",
    "$$\\hat{g}^{ij} \\hat{\\Gamma}^k_{ij}$$\n",
    "in spherical coordinates, where $i\\in\\{1,2,3\\}$ and $j\\in\\{1,2,3\\}$ are spatial indices.\n",
    "\n",
    "As discussed in the [NRPy+ tutorial notebook on reference metrics](Tutorial-Reference_Metric.ipynb), several reference-metric-related quantities in spherical coordinates are computed in NRPy+ (provided the parameter **`reference_metric::CoordSystem`** is set to **`\"Spherical\"`**), including the inverse spatial spherical reference metric $\\hat{g}^{ij}$ and the Christoffel symbols from this reference metric $\\hat{\\Gamma}^{i}_{jk}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:22.880102Z",
     "iopub.status.busy": "2021-03-07T17:22:22.879343Z",
     "iopub.status.idle": "2021-03-07T17:22:23.694766Z",
     "shell.execute_reply": "2021-03-07T17:22:23.695286Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "contracted GammahatU[0]:\n",
      "-2/xx0\n",
      "\n",
      "\n",
      "\n",
      "contracted GammahatU[1]:\n",
      "-1/(xx0**2*tan(xx1))\n",
      "\n",
      "\n",
      "\n",
      "contracted GammahatU[2]:\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "import sympy as sp              # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import NRPy_param_funcs as par  # NRPy+: Parameter interface\n",
    "import grid as gri              # NRPy+: Functionality for handling numerical grids\n",
    "import indexedexp as ixp        # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import reference_metric as rfm  # NRPy+: Reference metric support\n",
    "\n",
    "# reference_metric::CoordSystem can be set to Spherical, SinhSpherical, SinhSphericalv2,\n",
    "#                           Cylindrical, SinhCylindrical, SinhCylindricalv2, etc.\n",
    "#                           See reference_metric.py and NRPy+ tutorial notebook on\n",
    "#                           reference metrics for full list and description of how\n",
    "#                           to extend.\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Spherical\")\n",
    "par.set_parval_from_str(\"grid::DIM\",3)\n",
    "\n",
    "rfm.reference_metric()\n",
    "\n",
    "contractedGammahatU = ixp.zerorank1()\n",
    "for k in range(3):\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            contractedGammahatU[k] += rfm.ghatUU[i][j] * rfm.GammahatUDD[k][i][j]\n",
    "\n",
    "for k in range(3):\n",
    "    print(\"contracted GammahatU[\"+str(k)+\"]:\")\n",
    "    print(sp.simplify(contractedGammahatU[k]))\n",
    "    # Sadly pretty_print results in garbage output in the generated PDF at the bottom of this notebook.\n",
    "#     sp.pretty_print(sp.simplify(contractedGammahatU[k]))\n",
    "    if k<2:\n",
    "        print(\"\\n\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='rhs_scalarwave_spherical'></a>\n",
    "\n",
    "# Step 2: The right-hand side of the scalar wave equation in spherical coordinates, using NRPy+ \\[Back to [top](#toc)\\]\n",
    "$$\\label{rhs_scalarwave_spherical}$$\n",
    "\n",
    "Following our [implementation of the scalar wave equation in Cartesian coordinates](Tutorial-ScalarWave.ipynb), we will introduce a new variable $v=\\partial_t u$ that will enable us to split the second time derivative into two first-order time derivatives:\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_t u &= v \\\\\n",
    "\\partial_t v &= \\hat{g}^{ij} \\partial_{i} \\partial_{j} u - \\hat{\\Gamma}^i \\partial_i u.\n",
    "\\end{align}\n",
    "\n",
    "Adding back the sound speed $c$, we have a choice of a single factor of $c$ multiplying both right-hand sides, or a factor of $c^2$ multiplying the second equation only. We'll choose the latter:\n",
    "\n",
    "\\begin{align}\n",
    "\\partial_t u &= v \\\\\n",
    "\\partial_t v &= c^2 \\left(\\hat{g}^{ij} \\partial_{i} \\partial_{j} u - \\hat{\\Gamma}^i \\partial_i u\\right).\n",
    "\\end{align}\n",
    "\n",
    "Now let's generate the C code for the finite-difference representations of the right-hand sides of the above \"time evolution\" equations for $u$ and $v$. Since the right-hand side of $\\partial_t v$ contains implied sums over $i$ and $j$ in the first term, and an implied sum over $k$ in the second term, we'll find it useful to split the right-hand side into two parts\n",
    "\n",
    "\\begin{equation}\n",
    "\\partial_t v = c^2 \\left(\n",
    "{\\underbrace {\\textstyle \\hat{g}^{ij} \\partial_{i} \\partial_{j} u}_{\\text{Part 1}}} \n",
    "{\\underbrace {\\textstyle -\\hat{\\Gamma}^i \\partial_i u}_{\\text{Part 2}}}\\right),\n",
    "\\end{equation}\n",
    "\n",
    "and perform the implied sums in two pieces:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:23.700515Z",
     "iopub.status.busy": "2021-03-07T17:22:23.699804Z",
     "iopub.status.idle": "2021-03-07T17:22:23.703494Z",
     "shell.execute_reply": "2021-03-07T17:22:23.704053Z"
    }
   },
   "outputs": [],
   "source": [
    "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",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "from outputC import lhrh         # NRPy+: Core C code output module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:23.773311Z",
     "iopub.status.busy": "2021-03-07T17:22:23.737378Z",
     "iopub.status.idle": "2021-03-07T17:22:23.879189Z",
     "shell.execute_reply": "2021-03-07T17:22:23.879729Z"
    },
    "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 double uu_dD0 = invdx0*(-2*uu_i0m1_i1_i2/3 + uu_i0m2_i1_i2/12 + 2*uu_i0p1_i1_i2/3 - uu_i0p2_i1_i2/12),\n",
      "    *    const double uu_dD1 = invdx1*(-2*uu_i0_i1m1_i2/3 + uu_i0_i1m2_i2/12 + 2*uu_i0_i1p1_i2/3 - uu_i0_i1p2_i2/12),\n",
      "    *    const double uu_dDD00 = invdx0**2*(-5*uu/2 + 4*uu_i0m1_i1_i2/3 - uu_i0m2_i1_i2/12 + 4*uu_i0p1_i1_i2/3 - uu_i0p2_i1_i2/12),\n",
      "    *    const double uu_dDD11 = invdx1**2*(-5*uu/2 + 4*uu_i0_i1m1_i2/3 - uu_i0_i1m2_i2/12 + 4*uu_i0_i1p1_i2/3 - uu_i0_i1p2_i2/12),\n",
      "    *    const double uu_dDD22 = invdx2**2*(-5*uu/2 + 4*uu_i0_i1_i2m1/3 - uu_i0_i1_i2m2/12 + 4*uu_i0_i1_i2p1/3 - uu_i0_i1_i2p2/12)]\"\n",
      "    */\n",
      "   const double uu_i0_i1_i2m2 = in_gfs[IDX4(UUGF, i0,i1,i2-2)];\n",
      "   const double uu_i0_i1_i2m1 = in_gfs[IDX4(UUGF, i0,i1,i2-1)];\n",
      "   const double uu_i0_i1m2_i2 = in_gfs[IDX4(UUGF, i0,i1-2,i2)];\n",
      "   const double uu_i0_i1m1_i2 = in_gfs[IDX4(UUGF, i0,i1-1,i2)];\n",
      "   const double uu_i0m2_i1_i2 = in_gfs[IDX4(UUGF, i0-2,i1,i2)];\n",
      "   const double uu_i0m1_i1_i2 = in_gfs[IDX4(UUGF, i0-1,i1,i2)];\n",
      "   const double uu = in_gfs[IDX4(UUGF, i0,i1,i2)];\n",
      "   const double uu_i0p1_i1_i2 = in_gfs[IDX4(UUGF, i0+1,i1,i2)];\n",
      "   const double uu_i0p2_i1_i2 = in_gfs[IDX4(UUGF, i0+2,i1,i2)];\n",
      "   const double uu_i0_i1p1_i2 = in_gfs[IDX4(UUGF, i0,i1+1,i2)];\n",
      "   const double uu_i0_i1p2_i2 = in_gfs[IDX4(UUGF, i0,i1+2,i2)];\n",
      "   const double uu_i0_i1_i2p1 = in_gfs[IDX4(UUGF, i0,i1,i2+1)];\n",
      "   const double uu_i0_i1_i2p2 = in_gfs[IDX4(UUGF, i0,i1,i2+2)];\n",
      "   const double vv = in_gfs[IDX4(VVGF, i0,i1,i2)];\n",
      "   const double FDPart1_Rational_2_3 = 2.0/3.0;\n",
      "   const double FDPart1_Rational_1_12 = 1.0/12.0;\n",
      "   const double FDPart1_Rational_5_2 = 5.0/2.0;\n",
      "   const double FDPart1_Rational_4_3 = 4.0/3.0;\n",
      "   const double FDPart1_1 = -uu_i0_i1p2_i2;\n",
      "   const double FDPart1_2 = -FDPart1_Rational_5_2*uu;\n",
      "   const double uu_dD0 = invdx0*(FDPart1_Rational_1_12*(uu_i0m2_i1_i2 - uu_i0p2_i1_i2) + FDPart1_Rational_2_3*(-uu_i0m1_i1_i2 + uu_i0p1_i1_i2));\n",
      "   const double uu_dD1 = invdx1*(FDPart1_Rational_1_12*(FDPart1_1 + uu_i0_i1m2_i2) + FDPart1_Rational_2_3*(-uu_i0_i1m1_i2 + uu_i0_i1p1_i2));\n",
      "   const double uu_dDD00 = ((invdx0)*(invdx0))*(FDPart1_2 + FDPart1_Rational_1_12*(-uu_i0m2_i1_i2 - uu_i0p2_i1_i2) + FDPart1_Rational_4_3*(uu_i0m1_i1_i2 + uu_i0p1_i1_i2));\n",
      "   const double uu_dDD11 = ((invdx1)*(invdx1))*(FDPart1_2 + FDPart1_Rational_1_12*(FDPart1_1 - uu_i0_i1m2_i2) + FDPart1_Rational_4_3*(uu_i0_i1m1_i2 + uu_i0_i1p1_i2));\n",
      "   const double uu_dDD22 = ((invdx2)*(invdx2))*(FDPart1_2 + FDPart1_Rational_1_12*(-uu_i0_i1_i2m2 - uu_i0_i1_i2p2) + FDPart1_Rational_4_3*(uu_i0_i1_i2m1 + uu_i0_i1_i2p1));\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[IDX4(UUGF, i0, i1, i2)] = vv,\n",
      "    *    rhs_gfs[IDX4(VVGF, i0, i1, i2)] = wavespeed**2*(2*uu_dD0/xx0 + uu_dD1*sin(2*xx1)/(2*xx0**2*sin(xx1)**2) + uu_dDD00 + uu_dDD11/xx0**2 + uu_dDD22/(xx0**2*sin(xx1)**2))]\"\n",
      "    */\n",
      "   const double FDPart3_0 = (1.0/((xx0)*(xx0)));\n",
      "   const double FDPart3_1 = FDPart3_0/((sin(xx1))*(sin(xx1)));\n",
      "   rhs_gfs[IDX4(UUGF, i0, i1, i2)] = vv;\n",
      "   rhs_gfs[IDX4(VVGF, i0, i1, i2)] = ((wavespeed)*(wavespeed))*(FDPart3_0*uu_dDD11 + (1.0/2.0)*FDPart3_1*uu_dD1*sin(2*xx1) + FDPart3_1*uu_dDD22 + 2*uu_dD0/xx0 + uu_dDD00);\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "# The name of this module (\"scalarwave\") is given by __name__:\n",
    "thismodule = __name__\n",
    "\n",
    "# Step 0: Read the spatial dimension parameter as DIM.\n",
    "DIM = par.parval_from_str(\"grid::DIM\")\n",
    "\n",
    "# Step 1: Set the finite differencing order to 4.\n",
    "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\",4)\n",
    "\n",
    "# Step 2a: Reset the gridfunctions list; below we define the\n",
    "#          full complement of gridfunctions needed by this\n",
    "#          tutorial. This line of code enables us to re-run this\n",
    "#          tutorial without resetting the running Python kernel.\n",
    "gri.glb_gridfcs_list = []\n",
    "# Step 2b: 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 3a: Declare the rank-1 indexed expression \\partial_{i} u,\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_dD = ixp.declarerank1(\"uu_dD\")\n",
    "\n",
    "# Step 3b: 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 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 5: Define right-hand sides for the evolution.\n",
    "uu_rhs = vv\n",
    "# Step 5b: The right-hand side of the \\partial_t v equation\n",
    "#          is given by:\n",
    "#          \\hat{g}^{ij} \\partial_i \\partial_j u - \\hat{\\Gamma}^i \\partial_i u.\n",
    "#          ^^^^^^^^^^^^ PART 1 ^^^^^^^^^^^^^^^^ ^^^^^^^^^^ PART 2 ^^^^^^^^^^^\n",
    "vv_rhs = 0\n",
    "for i in range(DIM):\n",
    "    # PART 2:\n",
    "    vv_rhs -= contractedGammahatU[i]*uu_dD[i]\n",
    "    for j in range(DIM):\n",
    "        # PART 1:\n",
    "        vv_rhs += rfm.ghatUU[i][j]*uu_dDD[i][j]\n",
    "\n",
    "vv_rhs *= wavespeed*wavespeed\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": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 3: Code Validation against `ScalarWave.ScalarWaveCurvilinear_RHSs`  NRPy+ Module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for the RHSs of the Curvilinear Scalar Wave equation (i.e., uu_rhs and vv_rhs) between\n",
    "\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [ScalarWave.ScalarWaveCurvilinear_RHSs](../edit/ScalarWaveCurvilinear/ScalarWaveCurvilinear_RHSs.py) module.\n",
    "\n",
    "By default, we analyze the RHSs in Spherical coordinates, though other coordinate systems may be chosen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:23.886982Z",
     "iopub.status.busy": "2021-03-07T17:22:23.886279Z",
     "iopub.status.idle": "2021-03-07T17:22:23.965627Z",
     "shell.execute_reply": "2021-03-07T17:22:23.966148Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between ScalarWaveCurvilinear tutorial and NRPy+ module:\n",
      "TESTS PASSED!\n"
     ]
    }
   ],
   "source": [
    "# Step 7: 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",
    "#         ScalarWaveCurvilinear/ScalarWaveCurvilinear module,\n",
    "#         to ensure consistency between the tutorial and the\n",
    "#         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 8: Call the ScalarWaveCurvilinear_RHSs() function from within the\n",
    "#         ScalarWaveCurvilinear/ScalarWaveCurvilinear_RHSs.py module,\n",
    "#         which should do exactly the same as in Steps 1-6 above.\n",
    "import ScalarWave.ScalarWaveCurvilinear_RHSs as swcrhs\n",
    "swcrhs.ScalarWaveCurvilinear_RHSs()\n",
    "\n",
    "# Step 9: Consistency check between the tutorial notebook above\n",
    "#         and the ScalarWaveCurvilinear_RHSs() function from within the\n",
    "#         ScalarWaveCurvilinear/ScalarWaveCurvilinear_RHSs.py module.\n",
    "print(\"Consistency check between ScalarWaveCurvilinear tutorial and NRPy+ module:\")\n",
    "if sp.simplify(uu_rhs - swcrhs.uu_rhs) != 0:\n",
    "    print(\"TEST FAILED: uu_ID_SphericalGaussian - swid.uu_ID = \"+str(sp.simplify(uu_rhs - swcrhs.uu_rhs))+\"\\t\\t (should be zero)\")\n",
    "    sys.exit(1)\n",
    "if sp.simplify(vv_rhs - swcrhs.vv_rhs) != 0:\n",
    "    print(\"TEST FAILED: vv_ID_SphericalGaussian - swid.vv_ID = \"+str(sp.simplify(vv_rhs - swcrhs.vv_rhs))+\"\\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 4: 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-ScalarWaveCurvilinear.pdf](Tutorial-ScalarWaveCurvilinear.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": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:22:23.971108Z",
     "iopub.status.busy": "2021-03-07T17:22:23.970267Z",
     "iopub.status.idle": "2021-03-07T17:22:27.454842Z",
     "shell.execute_reply": "2021-03-07T17:22:27.455420Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-ScalarWaveCurvilinear.tex, and compiled LaTeX file to PDF\n",
      "    file Tutorial-ScalarWaveCurvilinear.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ScalarWaveCurvilinear\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}