Tutorial-MaxwellCurvilinear.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 sides of Maxwell's equations, in ***curvilinear*** coordinates, using a reference metric formalism\n",
    "\n",
    "## Author: Ian Ruchlin\n",
    "### Formatting improvements courtesy Brandon Clark\n",
    "\n",
    "[comment]: <> (Abstract: TODO)\n",
    "\n",
    "### The following formulations of Maxwell's equations, called System I and System II, are described in [Illustrating Stability Properties of Numerical Relativity in Electrodynamics](https://arxiv.org/abs/gr-qc/0201051) by Knapp et al.\n",
    "\n",
    "**Notebook Status:** <font color='red'><b> In progress </b></font>\n",
    "\n",
    "**Validation Notes:** This module has not yet undergone validation testing. Do ***not*** use it until after appropriate validation testing has been performed.\n",
    "\n",
    "## Introduction:\n",
    "[Maxwell's equations](https://en.wikipedia.org/wiki/Maxwell%27s_equations) are subject to the Gauss' law constraint\n",
    "$$\\mathcal{C} \\equiv \\hat{D}_{i} E^{i} - 4 \\pi \\rho = 0 \\; ,$$\n",
    "where $E^{i}$ is the electric vector field, $\\hat{D}_{i}$ is the [covariant derivative](https://en.wikipedia.org/wiki/Covariant_derivative) associated with the reference metric $\\hat{\\gamma}_{i j}$ (which is taken to represent flat space), and $\\rho$ is the electric charge density. We use $\\mathcal{C}$ as a measure of numerical error. Maxwell's equations are also required to satisfy $\\hat{D}_{i} B^{i} = 0$, where $B^{i}$ is the magnetic vector field. The magnetic constraint implies that the magnetic field can be expressed as\n",
    "$$B_{i} = \\epsilon_{i j k} \\hat{D}^{j} A^{k} \\; ,$$\n",
    "where $\\epsilon_{i j k}$ is the totally antisymmetric [Levi-Civita tensor](https://en.wikipedia.org/wiki/Levi-Civita_symbol) and $A^{i}$ is the vector potential field. Together with the scalar potential $\\psi$, the electric field can be expressed in terms of the potential fields as\n",
    "$$E_{i} = -\\hat{D}_{i} \\psi - \\partial_{t} A_{i} \\; .$$\n",
    "For now, we work in vacuum, where the electric charge density and the electric current density vector both vanish ($\\rho = 0$ and $j_{i} = 0$).\n",
    "\n",
    "In addition to the Gauss constraints, the electric and magnetic fields obey two independent [electromagnetic invariants](https://en.wikipedia.org/wiki/Classification_of_electromagnetic_fields#Invariants)\n",
    "\\begin{align}\n",
    "\\mathcal{P} &\\equiv B_{i} B^{i} - E_{i} E^{i} \\; , \\\\\n",
    "\\mathcal{Q} &\\equiv E_{i} B^{i} \\; .\n",
    "\\end{align}\n",
    "In vacuum, these satisfy $\\mathcal{P} = \\mathcal{Q} = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents:\n",
    "$$\\label{toc}$$\n",
    "\n",
    "This notebook is organized as follows\n",
    "\n",
    "1. [Step 1](#sys1): System I\n",
    "1. [Step 2](#sys2): System II\n",
    "1. [Step 3](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sys1'></a>\n",
    "\n",
    "# Step 1: System I \\[Back to [top](#toc)\\]\n",
    "$$\\label{sys1}$$\n",
    "\n",
    "In terms of the above definitions, the evolution Maxwell's equations take the form\n",
    "\\begin{align}\n",
    "\\partial_{t} A_{i} &= -E_{i} - \\hat{D}_{i} \\psi \\; , \\\\\n",
    "\\partial_{t} E_{i} &= -\\hat{D}_{j} \\hat{D}^{j} A_{i} + \\hat{D}_{i} \\hat{D}_{j} A^{j}\\; , \\\\\n",
    "\\partial_{t} \\psi &= -\\hat{D}_{i} A^{i} \\; .\n",
    "\\end{align}\n",
    "Note that this coupled system contains mixed second derivatives in the second term on the right hand side of the $E^{i}$ evolution equation. We will revisit this fact when building System II.\n",
    "\n",
    "It can be shown that the Gauss constraint satisfies the evolution equation\n",
    "$$\\partial_{t} \\mathcal{C} = 0 \\; .$$\n",
    "This implies that any constraint violating numerical error remains fixed in place during the evolution. This becomes problematic when the violations grow large and spoil the physics of the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:14:52.471716Z",
     "iopub.status.busy": "2021-03-07T17:14:52.466315Z",
     "iopub.status.idle": "2021-03-07T17:14:54.475005Z",
     "shell.execute_reply": "2021-03-07T17:14:54.475489Z"
    }
   },
   "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 AD_dD00 = invdx0*(-2*AD0_i0m1_i1_i2/3 + AD0_i0m2_i1_i2/12 + 2*AD0_i0p1_i1_i2/3 - AD0_i0p2_i1_i2/12),\n",
      "    *    const double AD_dD01 = invdx1*(-2*AD0_i0_i1m1_i2/3 + AD0_i0_i1m2_i2/12 + 2*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),\n",
      "    *    const double AD_dD02 = invdx2*(-2*AD0_i0_i1_i2m1/3 + AD0_i0_i1_i2m2/12 + 2*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),\n",
      "    *    const double AD_dD10 = invdx0*(-2*AD1_i0m1_i1_i2/3 + AD1_i0m2_i1_i2/12 + 2*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),\n",
      "    *    const double AD_dD11 = invdx1*(-2*AD1_i0_i1m1_i2/3 + AD1_i0_i1m2_i2/12 + 2*AD1_i0_i1p1_i2/3 - AD1_i0_i1p2_i2/12),\n",
      "    *    const double AD_dD12 = invdx2*(-2*AD1_i0_i1_i2m1/3 + AD1_i0_i1_i2m2/12 + 2*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),\n",
      "    *    const double AD_dD20 = invdx0*(-2*AD2_i0m1_i1_i2/3 + AD2_i0m2_i1_i2/12 + 2*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),\n",
      "    *    const double AD_dD21 = invdx1*(-2*AD2_i0_i1m1_i2/3 + AD2_i0_i1m2_i2/12 + 2*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),\n",
      "    *    const double AD_dD22 = invdx2*(-2*AD2_i0_i1_i2m1/3 + AD2_i0_i1_i2m2/12 + 2*AD2_i0_i1_i2p1/3 - AD2_i0_i1_i2p2/12),\n",
      "    *    const double AD_dDD001 = invdx0*invdx1*(4*AD0_i0m1_i1m1_i2/9 - AD0_i0m1_i1m2_i2/18 - 4*AD0_i0m1_i1p1_i2/9 + AD0_i0m1_i1p2_i2/18 - AD0_i0m2_i1m1_i2/18 + AD0_i0m2_i1m2_i2/144 + AD0_i0m2_i1p1_i2/18 - AD0_i0m2_i1p2_i2/144 - 4*AD0_i0p1_i1m1_i2/9 + AD0_i0p1_i1m2_i2/18 + 4*AD0_i0p1_i1p1_i2/9 - AD0_i0p1_i1p2_i2/18 + AD0_i0p2_i1m1_i2/18 - AD0_i0p2_i1m2_i2/144 - AD0_i0p2_i1p1_i2/18 + AD0_i0p2_i1p2_i2/144),\n",
      "    *    const double AD_dDD002 = invdx0*invdx2*(4*AD0_i0m1_i1_i2m1/9 - AD0_i0m1_i1_i2m2/18 - 4*AD0_i0m1_i1_i2p1/9 + AD0_i0m1_i1_i2p2/18 - AD0_i0m2_i1_i2m1/18 + AD0_i0m2_i1_i2m2/144 + AD0_i0m2_i1_i2p1/18 - AD0_i0m2_i1_i2p2/144 - 4*AD0_i0p1_i1_i2m1/9 + AD0_i0p1_i1_i2m2/18 + 4*AD0_i0p1_i1_i2p1/9 - AD0_i0p1_i1_i2p2/18 + AD0_i0p2_i1_i2m1/18 - AD0_i0p2_i1_i2m2/144 - AD0_i0p2_i1_i2p1/18 + AD0_i0p2_i1_i2p2/144),\n",
      "    *    const double AD_dDD011 = invdx1**2*(-5*AD0/2 + 4*AD0_i0_i1m1_i2/3 - AD0_i0_i1m2_i2/12 + 4*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),\n",
      "    *    const double AD_dDD022 = invdx2**2*(-5*AD0/2 + 4*AD0_i0_i1_i2m1/3 - AD0_i0_i1_i2m2/12 + 4*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),\n",
      "    *    const double AD_dDD100 = invdx0**2*(-5*AD1/2 + 4*AD1_i0m1_i1_i2/3 - AD1_i0m2_i1_i2/12 + 4*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),\n",
      "    *    const double AD_dDD101 = invdx0*invdx1*(4*AD1_i0m1_i1m1_i2/9 - AD1_i0m1_i1m2_i2/18 - 4*AD1_i0m1_i1p1_i2/9 + AD1_i0m1_i1p2_i2/18 - AD1_i0m2_i1m1_i2/18 + AD1_i0m2_i1m2_i2/144 + AD1_i0m2_i1p1_i2/18 - AD1_i0m2_i1p2_i2/144 - 4*AD1_i0p1_i1m1_i2/9 + AD1_i0p1_i1m2_i2/18 + 4*AD1_i0p1_i1p1_i2/9 - AD1_i0p1_i1p2_i2/18 + AD1_i0p2_i1m1_i2/18 - AD1_i0p2_i1m2_i2/144 - AD1_i0p2_i1p1_i2/18 + AD1_i0p2_i1p2_i2/144),\n",
      "    *    const double AD_dDD112 = invdx1*invdx2*(4*AD1_i0_i1m1_i2m1/9 - AD1_i0_i1m1_i2m2/18 - 4*AD1_i0_i1m1_i2p1/9 + AD1_i0_i1m1_i2p2/18 - AD1_i0_i1m2_i2m1/18 + AD1_i0_i1m2_i2m2/144 + AD1_i0_i1m2_i2p1/18 - AD1_i0_i1m2_i2p2/144 - 4*AD1_i0_i1p1_i2m1/9 + AD1_i0_i1p1_i2m2/18 + 4*AD1_i0_i1p1_i2p1/9 - AD1_i0_i1p1_i2p2/18 + AD1_i0_i1p2_i2m1/18 - AD1_i0_i1p2_i2m2/144 - AD1_i0_i1p2_i2p1/18 + AD1_i0_i1p2_i2p2/144),\n",
      "    *    const double AD_dDD122 = invdx2**2*(-5*AD1/2 + 4*AD1_i0_i1_i2m1/3 - AD1_i0_i1_i2m2/12 + 4*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),\n",
      "    *    const double AD_dDD200 = invdx0**2*(-5*AD2/2 + 4*AD2_i0m1_i1_i2/3 - AD2_i0m2_i1_i2/12 + 4*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),\n",
      "    *    const double AD_dDD202 = invdx0*invdx2*(4*AD2_i0m1_i1_i2m1/9 - AD2_i0m1_i1_i2m2/18 - 4*AD2_i0m1_i1_i2p1/9 + AD2_i0m1_i1_i2p2/18 - AD2_i0m2_i1_i2m1/18 + AD2_i0m2_i1_i2m2/144 + AD2_i0m2_i1_i2p1/18 - AD2_i0m2_i1_i2p2/144 - 4*AD2_i0p1_i1_i2m1/9 + AD2_i0p1_i1_i2m2/18 + 4*AD2_i0p1_i1_i2p1/9 - AD2_i0p1_i1_i2p2/18 + AD2_i0p2_i1_i2m1/18 - AD2_i0p2_i1_i2m2/144 - AD2_i0p2_i1_i2p1/18 + AD2_i0p2_i1_i2p2/144),\n",
      "    *    const double AD_dDD211 = invdx1**2*(-5*AD2/2 + 4*AD2_i0_i1m1_i2/3 - AD2_i0_i1m2_i2/12 + 4*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),\n",
      "    *    const double AD_dDD212 = invdx1*invdx2*(4*AD2_i0_i1m1_i2m1/9 - AD2_i0_i1m1_i2m2/18 - 4*AD2_i0_i1m1_i2p1/9 + AD2_i0_i1m1_i2p2/18 - AD2_i0_i1m2_i2m1/18 + AD2_i0_i1m2_i2m2/144 + AD2_i0_i1m2_i2p1/18 - AD2_i0_i1m2_i2p2/144 - 4*AD2_i0_i1p1_i2m1/9 + AD2_i0_i1p1_i2m2/18 + 4*AD2_i0_i1p1_i2p1/9 - AD2_i0_i1p1_i2p2/18 + AD2_i0_i1p2_i2m1/18 - AD2_i0_i1p2_i2m2/144 - AD2_i0_i1p2_i2p1/18 + AD2_i0_i1p2_i2p2/144),\n",
      "    *    const double psi_dD0 = invdx0*(-2*psi_i0m1_i1_i2/3 + psi_i0m2_i1_i2/12 + 2*psi_i0p1_i1_i2/3 - psi_i0p2_i1_i2/12),\n",
      "    *    const double psi_dD1 = invdx1*(-2*psi_i0_i1m1_i2/3 + psi_i0_i1m2_i2/12 + 2*psi_i0_i1p1_i2/3 - psi_i0_i1p2_i2/12),\n",
      "    *    const double psi_dD2 = invdx2*(-2*psi_i0_i1_i2m1/3 + psi_i0_i1_i2m2/12 + 2*psi_i0_i1_i2p1/3 - psi_i0_i1_i2p2/12)]\"\n",
      "    */\n",
      "   const double psi_i0_i1_i2m2 = in_gfs[IDX4(PSIGF, i0,i1,i2-2)];\n",
      "   const double psi_i0_i1_i2m1 = in_gfs[IDX4(PSIGF, i0,i1,i2-1)];\n",
      "   const double psi_i0_i1m2_i2 = in_gfs[IDX4(PSIGF, i0,i1-2,i2)];\n",
      "   const double psi_i0_i1m1_i2 = in_gfs[IDX4(PSIGF, i0,i1-1,i2)];\n",
      "   const double psi_i0m2_i1_i2 = in_gfs[IDX4(PSIGF, i0-2,i1,i2)];\n",
      "   const double psi_i0m1_i1_i2 = in_gfs[IDX4(PSIGF, i0-1,i1,i2)];\n",
      "   const double psi_i0p1_i1_i2 = in_gfs[IDX4(PSIGF, i0+1,i1,i2)];\n",
      "   const double psi_i0p2_i1_i2 = in_gfs[IDX4(PSIGF, i0+2,i1,i2)];\n",
      "   const double psi_i0_i1p1_i2 = in_gfs[IDX4(PSIGF, i0,i1+1,i2)];\n",
      "   const double psi_i0_i1p2_i2 = in_gfs[IDX4(PSIGF, i0,i1+2,i2)];\n",
      "   const double psi_i0_i1_i2p1 = in_gfs[IDX4(PSIGF, i0,i1,i2+1)];\n",
      "   const double psi_i0_i1_i2p2 = in_gfs[IDX4(PSIGF, i0,i1,i2+2)];\n",
      "   const double ED0 = in_gfs[IDX4(ED0GF, i0,i1,i2)];\n",
      "   const double ED1 = in_gfs[IDX4(ED1GF, i0,i1,i2)];\n",
      "   const double ED2 = in_gfs[IDX4(ED2GF, i0,i1,i2)];\n",
      "   const double AD0_i0m2_i1_i2m2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2-2)];\n",
      "   const double AD0_i0m1_i1_i2m2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2-2)];\n",
      "   const double AD0_i0_i1_i2m2 = in_gfs[IDX4(AD0GF, i0,i1,i2-2)];\n",
      "   const double AD0_i0p1_i1_i2m2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2-2)];\n",
      "   const double AD0_i0p2_i1_i2m2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2-2)];\n",
      "   const double AD0_i0m2_i1_i2m1 = in_gfs[IDX4(AD0GF, i0-2,i1,i2-1)];\n",
      "   const double AD0_i0m1_i1_i2m1 = in_gfs[IDX4(AD0GF, i0-1,i1,i2-1)];\n",
      "   const double AD0_i0_i1_i2m1 = in_gfs[IDX4(AD0GF, i0,i1,i2-1)];\n",
      "   const double AD0_i0p1_i1_i2m1 = in_gfs[IDX4(AD0GF, i0+1,i1,i2-1)];\n",
      "   const double AD0_i0p2_i1_i2m1 = in_gfs[IDX4(AD0GF, i0+2,i1,i2-1)];\n",
      "   const double AD0_i0m2_i1m2_i2 = in_gfs[IDX4(AD0GF, i0-2,i1-2,i2)];\n",
      "   const double AD0_i0m1_i1m2_i2 = in_gfs[IDX4(AD0GF, i0-1,i1-2,i2)];\n",
      "   const double AD0_i0_i1m2_i2 = in_gfs[IDX4(AD0GF, i0,i1-2,i2)];\n",
      "   const double AD0_i0p1_i1m2_i2 = in_gfs[IDX4(AD0GF, i0+1,i1-2,i2)];\n",
      "   const double AD0_i0p2_i1m2_i2 = in_gfs[IDX4(AD0GF, i0+2,i1-2,i2)];\n",
      "   const double AD0_i0m2_i1m1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1-1,i2)];\n",
      "   const double AD0_i0m1_i1m1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1-1,i2)];\n",
      "   const double AD0_i0_i1m1_i2 = in_gfs[IDX4(AD0GF, i0,i1-1,i2)];\n",
      "   const double AD0_i0p1_i1m1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1-1,i2)];\n",
      "   const double AD0_i0p2_i1m1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1-1,i2)];\n",
      "   const double AD0_i0m2_i1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2)];\n",
      "   const double AD0_i0m1_i1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2)];\n",
      "   const double AD0 = in_gfs[IDX4(AD0GF, i0,i1,i2)];\n",
      "   const double AD0_i0p1_i1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2)];\n",
      "   const double AD0_i0p2_i1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2)];\n",
      "   const double AD0_i0m2_i1p1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1+1,i2)];\n",
      "   const double AD0_i0m1_i1p1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1+1,i2)];\n",
      "   const double AD0_i0_i1p1_i2 = in_gfs[IDX4(AD0GF, i0,i1+1,i2)];\n",
      "   const double AD0_i0p1_i1p1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1+1,i2)];\n",
      "   const double AD0_i0p2_i1p1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1+1,i2)];\n",
      "   const double AD0_i0m2_i1p2_i2 = in_gfs[IDX4(AD0GF, i0-2,i1+2,i2)];\n",
      "   const double AD0_i0m1_i1p2_i2 = in_gfs[IDX4(AD0GF, i0-1,i1+2,i2)];\n",
      "   const double AD0_i0_i1p2_i2 = in_gfs[IDX4(AD0GF, i0,i1+2,i2)];\n",
      "   const double AD0_i0p1_i1p2_i2 = in_gfs[IDX4(AD0GF, i0+1,i1+2,i2)];\n",
      "   const double AD0_i0p2_i1p2_i2 = in_gfs[IDX4(AD0GF, i0+2,i1+2,i2)];\n",
      "   const double AD0_i0m2_i1_i2p1 = in_gfs[IDX4(AD0GF, i0-2,i1,i2+1)];\n",
      "   const double AD0_i0m1_i1_i2p1 = in_gfs[IDX4(AD0GF, i0-1,i1,i2+1)];\n",
      "   const double AD0_i0_i1_i2p1 = in_gfs[IDX4(AD0GF, i0,i1,i2+1)];\n",
      "   const double AD0_i0p1_i1_i2p1 = in_gfs[IDX4(AD0GF, i0+1,i1,i2+1)];\n",
      "   const double AD0_i0p2_i1_i2p1 = in_gfs[IDX4(AD0GF, i0+2,i1,i2+1)];\n",
      "   const double AD0_i0m2_i1_i2p2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2+2)];\n",
      "   const double AD0_i0m1_i1_i2p2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2+2)];\n",
      "   const double AD0_i0_i1_i2p2 = in_gfs[IDX4(AD0GF, i0,i1,i2+2)];\n",
      "   const double AD0_i0p1_i1_i2p2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2+2)];\n",
      "   const double AD0_i0p2_i1_i2p2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2+2)];\n",
      "   const double AD1_i0_i1m2_i2m2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2-2)];\n",
      "   const double AD1_i0_i1m1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2-2)];\n",
      "   const double AD1_i0_i1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1,i2-2)];\n",
      "   const double AD1_i0_i1p1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2-2)];\n",
      "   const double AD1_i0_i1p2_i2m2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2-2)];\n",
      "   const double AD1_i0_i1m2_i2m1 = in_gfs[IDX4(AD1GF, i0,i1-2,i2-1)];\n",
      "   const double AD1_i0_i1m1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1-1,i2-1)];\n",
      "   const double AD1_i0_i1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1,i2-1)];\n",
      "   const double AD1_i0_i1p1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1+1,i2-1)];\n",
      "   const double AD1_i0_i1p2_i2m1 = in_gfs[IDX4(AD1GF, i0,i1+2,i2-1)];\n",
      "   const double AD1_i0m2_i1m2_i2 = in_gfs[IDX4(AD1GF, i0-2,i1-2,i2)];\n",
      "   const double AD1_i0m1_i1m2_i2 = in_gfs[IDX4(AD1GF, i0-1,i1-2,i2)];\n",
      "   const double AD1_i0_i1m2_i2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2)];\n",
      "   const double AD1_i0p1_i1m2_i2 = in_gfs[IDX4(AD1GF, i0+1,i1-2,i2)];\n",
      "   const double AD1_i0p2_i1m2_i2 = in_gfs[IDX4(AD1GF, i0+2,i1-2,i2)];\n",
      "   const double AD1_i0m2_i1m1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1-1,i2)];\n",
      "   const double AD1_i0m1_i1m1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1-1,i2)];\n",
      "   const double AD1_i0_i1m1_i2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2)];\n",
      "   const double AD1_i0p1_i1m1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1-1,i2)];\n",
      "   const double AD1_i0p2_i1m1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1-1,i2)];\n",
      "   const double AD1_i0m2_i1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1,i2)];\n",
      "   const double AD1_i0m1_i1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1,i2)];\n",
      "   const double AD1 = in_gfs[IDX4(AD1GF, i0,i1,i2)];\n",
      "   const double AD1_i0p1_i1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1,i2)];\n",
      "   const double AD1_i0p2_i1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1,i2)];\n",
      "   const double AD1_i0m2_i1p1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1+1,i2)];\n",
      "   const double AD1_i0m1_i1p1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1+1,i2)];\n",
      "   const double AD1_i0_i1p1_i2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2)];\n",
      "   const double AD1_i0p1_i1p1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1+1,i2)];\n",
      "   const double AD1_i0p2_i1p1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1+1,i2)];\n",
      "   const double AD1_i0m2_i1p2_i2 = in_gfs[IDX4(AD1GF, i0-2,i1+2,i2)];\n",
      "   const double AD1_i0m1_i1p2_i2 = in_gfs[IDX4(AD1GF, i0-1,i1+2,i2)];\n",
      "   const double AD1_i0_i1p2_i2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2)];\n",
      "   const double AD1_i0p1_i1p2_i2 = in_gfs[IDX4(AD1GF, i0+1,i1+2,i2)];\n",
      "   const double AD1_i0p2_i1p2_i2 = in_gfs[IDX4(AD1GF, i0+2,i1+2,i2)];\n",
      "   const double AD1_i0_i1m2_i2p1 = in_gfs[IDX4(AD1GF, i0,i1-2,i2+1)];\n",
      "   const double AD1_i0_i1m1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1-1,i2+1)];\n",
      "   const double AD1_i0_i1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1,i2+1)];\n",
      "   const double AD1_i0_i1p1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1+1,i2+1)];\n",
      "   const double AD1_i0_i1p2_i2p1 = in_gfs[IDX4(AD1GF, i0,i1+2,i2+1)];\n",
      "   const double AD1_i0_i1m2_i2p2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2+2)];\n",
      "   const double AD1_i0_i1m1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2+2)];\n",
      "   const double AD1_i0_i1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1,i2+2)];\n",
      "   const double AD1_i0_i1p1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2+2)];\n",
      "   const double AD1_i0_i1p2_i2p2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2+2)];\n",
      "   const double AD2_i0_i1m2_i2m2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2-2)];\n",
      "   const double AD2_i0_i1m1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2-2)];\n",
      "   const double AD2_i0m2_i1_i2m2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2-2)];\n",
      "   const double AD2_i0m1_i1_i2m2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2-2)];\n",
      "   const double AD2_i0_i1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1,i2-2)];\n",
      "   const double AD2_i0p1_i1_i2m2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2-2)];\n",
      "   const double AD2_i0p2_i1_i2m2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2-2)];\n",
      "   const double AD2_i0_i1p1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2-2)];\n",
      "   const double AD2_i0_i1p2_i2m2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2-2)];\n",
      "   const double AD2_i0_i1m2_i2m1 = in_gfs[IDX4(AD2GF, i0,i1-2,i2-1)];\n",
      "   const double AD2_i0_i1m1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1-1,i2-1)];\n",
      "   const double AD2_i0m2_i1_i2m1 = in_gfs[IDX4(AD2GF, i0-2,i1,i2-1)];\n",
      "   const double AD2_i0m1_i1_i2m1 = in_gfs[IDX4(AD2GF, i0-1,i1,i2-1)];\n",
      "   const double AD2_i0_i1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1,i2-1)];\n",
      "   const double AD2_i0p1_i1_i2m1 = in_gfs[IDX4(AD2GF, i0+1,i1,i2-1)];\n",
      "   const double AD2_i0p2_i1_i2m1 = in_gfs[IDX4(AD2GF, i0+2,i1,i2-1)];\n",
      "   const double AD2_i0_i1p1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1+1,i2-1)];\n",
      "   const double AD2_i0_i1p2_i2m1 = in_gfs[IDX4(AD2GF, i0,i1+2,i2-1)];\n",
      "   const double AD2_i0_i1m2_i2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2)];\n",
      "   const double AD2_i0_i1m1_i2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2)];\n",
      "   const double AD2_i0m2_i1_i2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2)];\n",
      "   const double AD2_i0m1_i1_i2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2)];\n",
      "   const double AD2 = in_gfs[IDX4(AD2GF, i0,i1,i2)];\n",
      "   const double AD2_i0p1_i1_i2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2)];\n",
      "   const double AD2_i0p2_i1_i2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2)];\n",
      "   const double AD2_i0_i1p1_i2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2)];\n",
      "   const double AD2_i0_i1p2_i2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2)];\n",
      "   const double AD2_i0_i1m2_i2p1 = in_gfs[IDX4(AD2GF, i0,i1-2,i2+1)];\n",
      "   const double AD2_i0_i1m1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1-1,i2+1)];\n",
      "   const double AD2_i0m2_i1_i2p1 = in_gfs[IDX4(AD2GF, i0-2,i1,i2+1)];\n",
      "   const double AD2_i0m1_i1_i2p1 = in_gfs[IDX4(AD2GF, i0-1,i1,i2+1)];\n",
      "   const double AD2_i0_i1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1,i2+1)];\n",
      "   const double AD2_i0p1_i1_i2p1 = in_gfs[IDX4(AD2GF, i0+1,i1,i2+1)];\n",
      "   const double AD2_i0p2_i1_i2p1 = in_gfs[IDX4(AD2GF, i0+2,i1,i2+1)];\n",
      "   const double AD2_i0_i1p1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1+1,i2+1)];\n",
      "   const double AD2_i0_i1p2_i2p1 = in_gfs[IDX4(AD2GF, i0,i1+2,i2+1)];\n",
      "   const double AD2_i0_i1m2_i2p2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2+2)];\n",
      "   const double AD2_i0_i1m1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2+2)];\n",
      "   const double AD2_i0m2_i1_i2p2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2+2)];\n",
      "   const double AD2_i0m1_i1_i2p2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2+2)];\n",
      "   const double AD2_i0_i1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1,i2+2)];\n",
      "   const double AD2_i0p1_i1_i2p2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2+2)];\n",
      "   const double AD2_i0p2_i1_i2p2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2+2)];\n",
      "   const double AD2_i0_i1p1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2+2)];\n",
      "   const double AD2_i0_i1p2_i2p2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2+2)];\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_4_9 = 4.0/9.0;\n",
      "   const double FDPart1_Rational_1_18 = 1.0/18.0;\n",
      "   const double FDPart1_Rational_1_144 = 1.0/144.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 = -AD0_i0_i1_i2p2;\n",
      "   const double FDPart1_9 = -AD0*FDPart1_Rational_5_2;\n",
      "   const double FDPart1_12 = -AD1*FDPart1_Rational_5_2;\n",
      "   const double FDPart1_14 = -AD2*FDPart1_Rational_5_2;\n",
      "   const double AD_dD00 = invdx0*(FDPart1_Rational_1_12*(AD0_i0m2_i1_i2 - AD0_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD0_i0m1_i1_i2 + AD0_i0p1_i1_i2));\n",
      "   const double AD_dD01 = invdx1*(FDPart1_Rational_1_12*(AD0_i0_i1m2_i2 - AD0_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));\n",
      "   const double AD_dD02 = invdx2*(FDPart1_Rational_1_12*(AD0_i0_i1_i2m2 + FDPart1_1) + FDPart1_Rational_2_3*(-AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));\n",
      "   const double AD_dD10 = invdx0*(FDPart1_Rational_1_12*(AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));\n",
      "   const double AD_dD11 = invdx1*(FDPart1_Rational_1_12*(AD1_i0_i1m2_i2 - AD1_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD1_i0_i1m1_i2 + AD1_i0_i1p1_i2));\n",
      "   const double AD_dD12 = invdx2*(FDPart1_Rational_1_12*(AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));\n",
      "   const double AD_dD20 = invdx0*(FDPart1_Rational_1_12*(AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));\n",
      "   const double AD_dD21 = invdx1*(FDPart1_Rational_1_12*(AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));\n",
      "   const double AD_dD22 = invdx2*(FDPart1_Rational_1_12*(AD2_i0_i1_i2m2 - AD2_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD2_i0_i1_i2m1 + AD2_i0_i1_i2p1));\n",
      "   const double AD_dDD001 = invdx0*invdx1*(FDPart1_Rational_1_144*(AD0_i0m2_i1m2_i2 - AD0_i0m2_i1p2_i2 - AD0_i0p2_i1m2_i2 + AD0_i0p2_i1p2_i2) + FDPart1_Rational_1_18*(-AD0_i0m1_i1m2_i2 + AD0_i0m1_i1p2_i2 - AD0_i0m2_i1m1_i2 + AD0_i0m2_i1p1_i2 + AD0_i0p1_i1m2_i2 - AD0_i0p1_i1p2_i2 + AD0_i0p2_i1m1_i2 - AD0_i0p2_i1p1_i2) + FDPart1_Rational_4_9*(AD0_i0m1_i1m1_i2 - AD0_i0m1_i1p1_i2 - AD0_i0p1_i1m1_i2 + AD0_i0p1_i1p1_i2));\n",
      "   const double AD_dDD002 = invdx0*invdx2*(FDPart1_Rational_1_144*(AD0_i0m2_i1_i2m2 - AD0_i0m2_i1_i2p2 - AD0_i0p2_i1_i2m2 + AD0_i0p2_i1_i2p2) + FDPart1_Rational_1_18*(-AD0_i0m1_i1_i2m2 + AD0_i0m1_i1_i2p2 - AD0_i0m2_i1_i2m1 + AD0_i0m2_i1_i2p1 + AD0_i0p1_i1_i2m2 - AD0_i0p1_i1_i2p2 + AD0_i0p2_i1_i2m1 - AD0_i0p2_i1_i2p1) + FDPart1_Rational_4_9*(AD0_i0m1_i1_i2m1 - AD0_i0m1_i1_i2p1 - AD0_i0p1_i1_i2m1 + AD0_i0p1_i1_i2p1));\n",
      "   const double AD_dDD011 = ((invdx1)*(invdx1))*(FDPart1_9 + FDPart1_Rational_1_12*(-AD0_i0_i1m2_i2 - AD0_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));\n",
      "   const double AD_dDD022 = ((invdx2)*(invdx2))*(FDPart1_9 + FDPart1_Rational_1_12*(-AD0_i0_i1_i2m2 + FDPart1_1) + FDPart1_Rational_4_3*(AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));\n",
      "   const double AD_dDD100 = ((invdx0)*(invdx0))*(FDPart1_12 + FDPart1_Rational_1_12*(-AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));\n",
      "   const double AD_dDD101 = invdx0*invdx1*(FDPart1_Rational_1_144*(AD1_i0m2_i1m2_i2 - AD1_i0m2_i1p2_i2 - AD1_i0p2_i1m2_i2 + AD1_i0p2_i1p2_i2) + FDPart1_Rational_1_18*(-AD1_i0m1_i1m2_i2 + AD1_i0m1_i1p2_i2 - AD1_i0m2_i1m1_i2 + AD1_i0m2_i1p1_i2 + AD1_i0p1_i1m2_i2 - AD1_i0p1_i1p2_i2 + AD1_i0p2_i1m1_i2 - AD1_i0p2_i1p1_i2) + FDPart1_Rational_4_9*(AD1_i0m1_i1m1_i2 - AD1_i0m1_i1p1_i2 - AD1_i0p1_i1m1_i2 + AD1_i0p1_i1p1_i2));\n",
      "   const double AD_dDD112 = invdx1*invdx2*(FDPart1_Rational_1_144*(AD1_i0_i1m2_i2m2 - AD1_i0_i1m2_i2p2 - AD1_i0_i1p2_i2m2 + AD1_i0_i1p2_i2p2) + FDPart1_Rational_1_18*(-AD1_i0_i1m1_i2m2 + AD1_i0_i1m1_i2p2 - AD1_i0_i1m2_i2m1 + AD1_i0_i1m2_i2p1 + AD1_i0_i1p1_i2m2 - AD1_i0_i1p1_i2p2 + AD1_i0_i1p2_i2m1 - AD1_i0_i1p2_i2p1) + FDPart1_Rational_4_9*(AD1_i0_i1m1_i2m1 - AD1_i0_i1m1_i2p1 - AD1_i0_i1p1_i2m1 + AD1_i0_i1p1_i2p1));\n",
      "   const double AD_dDD122 = ((invdx2)*(invdx2))*(FDPart1_12 + FDPart1_Rational_1_12*(-AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));\n",
      "   const double AD_dDD200 = ((invdx0)*(invdx0))*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));\n",
      "   const double AD_dDD202 = invdx0*invdx2*(FDPart1_Rational_1_144*(AD2_i0m2_i1_i2m2 - AD2_i0m2_i1_i2p2 - AD2_i0p2_i1_i2m2 + AD2_i0p2_i1_i2p2) + FDPart1_Rational_1_18*(-AD2_i0m1_i1_i2m2 + AD2_i0m1_i1_i2p2 - AD2_i0m2_i1_i2m1 + AD2_i0m2_i1_i2p1 + AD2_i0p1_i1_i2m2 - AD2_i0p1_i1_i2p2 + AD2_i0p2_i1_i2m1 - AD2_i0p2_i1_i2p1) + FDPart1_Rational_4_9*(AD2_i0m1_i1_i2m1 - AD2_i0m1_i1_i2p1 - AD2_i0p1_i1_i2m1 + AD2_i0p1_i1_i2p1));\n",
      "   const double AD_dDD211 = ((invdx1)*(invdx1))*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));\n",
      "   const double AD_dDD212 = invdx1*invdx2*(FDPart1_Rational_1_144*(AD2_i0_i1m2_i2m2 - AD2_i0_i1m2_i2p2 - AD2_i0_i1p2_i2m2 + AD2_i0_i1p2_i2p2) + FDPart1_Rational_1_18*(-AD2_i0_i1m1_i2m2 + AD2_i0_i1m1_i2p2 - AD2_i0_i1m2_i2m1 + AD2_i0_i1m2_i2p1 + AD2_i0_i1p1_i2m2 - AD2_i0_i1p1_i2p2 + AD2_i0_i1p2_i2m1 - AD2_i0_i1p2_i2p1) + FDPart1_Rational_4_9*(AD2_i0_i1m1_i2m1 - AD2_i0_i1m1_i2p1 - AD2_i0_i1p1_i2m1 + AD2_i0_i1p1_i2p1));\n",
      "   const double psi_dD0 = invdx0*(FDPart1_Rational_1_12*(psi_i0m2_i1_i2 - psi_i0p2_i1_i2) + FDPart1_Rational_2_3*(-psi_i0m1_i1_i2 + psi_i0p1_i1_i2));\n",
      "   const double psi_dD1 = invdx1*(FDPart1_Rational_1_12*(psi_i0_i1m2_i2 - psi_i0_i1p2_i2) + FDPart1_Rational_2_3*(-psi_i0_i1m1_i2 + psi_i0_i1p1_i2));\n",
      "   const double psi_dD2 = invdx2*(FDPart1_Rational_1_12*(psi_i0_i1_i2m2 - psi_i0_i1_i2p2) + FDPart1_Rational_2_3*(-psi_i0_i1_i2m1 + psi_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(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0,\n",
      "    *    rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = (AD0 + AD_dD00*xx0 + AD_dDD101 - 2*(AD0*xx0 + AD_dD11)/xx0)/xx0**2 - (AD_dD00*xx0 - AD_dD11/xx0 + AD_dDD011 - (AD0*xx0 + AD_dD11)/xx0)/xx0**2 + (AD0*sin(xx1)**2 + AD_dD00*xx0*sin(xx1)**2 + AD_dD10*sin(2*xx1)/2 + AD_dDD202 - 2*(AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/xx0)/(xx0**2*sin(xx1)**2) - (AD_dD00*xx0*sin(xx1)**2 - AD_dD22/xx0 + AD_dDD022 + (-AD1/xx0 + AD_dD01)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/xx0)/(xx0**2*sin(xx1)**2),\n",
      "    *    rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1,\n",
      "    *    rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1/xx0**2 + AD_dD10/xx0 + AD_dDD001 - AD_dDD100 - (-AD1/xx0 + AD_dD01)/xx0 - (-AD_dD22*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD122 + xx0*(-AD1/xx0 + AD_dD10)*sin(xx1)**2 + (AD0*xx0 + AD_dD11)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)*sin(2*xx1)/(2*sin(xx1)**2))/(xx0**2*sin(xx1)**2) + (2*AD0*xx0*sin(xx1)*cos(xx1) + AD1*cos(2*xx1) + AD_dD01*xx0*sin(xx1)**2 + AD_dD11*sin(2*xx1)/2 + AD_dDD212 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)*sin(2*xx1)/sin(xx1)**2)/(xx0**2*sin(xx1)**2),\n",
      "    *    rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2,\n",
      "    *    rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2/xx0**2 + AD_dD20/xx0 + AD_dDD002 - AD_dDD200 - (-AD2/xx0 + AD_dD02)/xx0 + (AD_dD02*xx0 + AD_dDD112 - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD12)*sin(2*xx1)/(2*sin(xx1)**2) - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/(2*sin(xx1)**2))/xx0**2 - (AD2*(-cos(2*xx1)/sin(xx1)**2 + sin(2*xx1)*cos(xx1)/sin(xx1)**3) - AD_dD21*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD211 + xx0*(-AD2/xx0 + AD_dD20) - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/(2*sin(xx1)**2))/xx0**2,\n",
      "    *    rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -AD_dD00 - (AD0*xx0 + AD_dD11)/xx0**2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/(xx0**2*sin(xx1)**2)]\"\n",
      "    */\n",
      "   const double FDPart3_0 = (1.0/((xx0)*(xx0)));\n",
      "   const double FDPart3_1 = AD_dD00*xx0;\n",
      "   const double FDPart3_2 = (1.0/(xx0));\n",
      "   const double FDPart3_3 = AD0*xx0;\n",
      "   const double FDPart3_4 = AD_dD11 + FDPart3_3;\n",
      "   const double FDPart3_6 = sin(xx1);\n",
      "   const double FDPart3_7 = ((FDPart3_6)*(FDPart3_6));\n",
      "   const double FDPart3_9 = sin(2*xx1);\n",
      "   const double FDPart3_10 = (1.0/2.0)*FDPart3_9;\n",
      "   const double FDPart3_12 = AD1*FDPart3_10 + AD_dD22 + FDPart3_3*FDPart3_7;\n",
      "   const double FDPart3_14 = (1.0/(FDPart3_7));\n",
      "   const double FDPart3_15 = FDPart3_0*FDPart3_14;\n",
      "   const double FDPart3_16 = -AD1*FDPart3_2;\n",
      "   const double FDPart3_18 = cos(2*xx1);\n",
      "   const double FDPart3_20 = cos(xx1);\n",
      "   const double FDPart3_22 = FDPart3_10*FDPart3_14;\n",
      "   const double FDPart3_23 = -AD2*FDPart3_2;\n",
      "   const double FDPart3_24 = -AD2*FDPart3_22;\n",
      "   const double FDPart3_25 = -FDPart3_22*(AD_dD21 + FDPart3_24);\n",
      "   rhs_gfs[IDX4(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0;\n",
      "   rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = FDPart3_0*(AD0 + AD_dDD101 + FDPart3_1 - 2*FDPart3_2*FDPart3_4) - FDPart3_0*(-AD_dD11*FDPart3_2 + AD_dDD011 + FDPart3_1 - FDPart3_2*FDPart3_4) + FDPart3_15*(AD0*FDPart3_7 + AD_dD10*FDPart3_10 + AD_dDD202 + FDPart3_1*FDPart3_7 - 2*FDPart3_12*FDPart3_2) - FDPart3_15*(-AD_dD22*FDPart3_2 + AD_dDD022 + FDPart3_1*FDPart3_7 + FDPart3_10*(AD_dD01 + FDPart3_16) - FDPart3_12*FDPart3_2);\n",
      "   rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1;\n",
      "   rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1*FDPart3_0 + AD_dD10*FDPart3_2 + AD_dDD001 - AD_dDD100 - FDPart3_15*(-AD_dD22*FDPart3_22 + AD_dDD122 - FDPart3_10*FDPart3_12*FDPart3_14 + FDPart3_10*FDPart3_4 + FDPart3_7*xx0*(AD_dD10 + FDPart3_16)) + FDPart3_15*(AD1*FDPart3_18 + AD_dD01*FDPart3_7*xx0 + AD_dD11*FDPart3_10 + AD_dDD212 - FDPart3_12*FDPart3_14*FDPart3_9 + 2*FDPart3_20*FDPart3_3*FDPart3_6) - FDPart3_2*(AD_dD01 + FDPart3_16);\n",
      "   rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2;\n",
      "   rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2*FDPart3_0 + AD_dD20*FDPart3_2 + AD_dDD002 - AD_dDD200 + FDPart3_0*(AD_dD02*xx0 + AD_dDD112 - FDPart3_22*(AD_dD12 + FDPart3_24) + FDPart3_25) - FDPart3_0*(AD2*(-FDPart3_14*FDPart3_18 + FDPart3_20*FDPart3_9/((FDPart3_6)*(FDPart3_6)*(FDPart3_6))) - AD_dD21*FDPart3_22 + AD_dDD211 + FDPart3_25 + xx0*(AD_dD20 + FDPart3_23)) - FDPart3_2*(AD_dD02 + FDPart3_23);\n",
      "   rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -AD_dD00 - FDPart3_0*FDPart3_4 - FDPart3_12*FDPart3_15;\n",
      "}\n"
     ]
    }
   ],
   "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\n",
    "\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",
    "# The name of this module (\"maxwell\") 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 2: Register gridfunctions that are needed as input.\n",
    "psi = gri.register_gridfunctions(\"EVOL\", [\"psi\"])\n",
    "\n",
    "# Step 3a: Declare the rank-1 indexed expressions E_{i}, A_{i},\n",
    "#          and \\partial_{i} \\psi. Derivative variables like these\n",
    "#          must have an underscore in them, so the finite\n",
    "#          difference module can parse the variable name properly.\n",
    "ED = ixp.register_gridfunctions_for_single_rank1(\"EVOL\", \"ED\")\n",
    "AD = ixp.register_gridfunctions_for_single_rank1(\"EVOL\", \"AD\")\n",
    "psi_dD = ixp.declarerank1(\"psi_dD\")\n",
    "\n",
    "# Step 3b: Declare the rank-2 indexed expression \\partial_{j} A_{i},\n",
    "#          which is not symmetric in its indices.\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",
    "AD_dD = ixp.declarerank2(\"AD_dD\", \"nosym\")\n",
    "\n",
    "# Step 3c: Declare the rank-3 indexed expression \\partial_{jk} A_{i},\n",
    "#          which is symmetric in the two {jk} indices.\n",
    "AD_dDD = ixp.declarerank3(\"AD_dDD\", \"sym12\")\n",
    "\n",
    "# Step 4: Calculate first and second covariant derivatives, and the\n",
    "#         necessary contractions.\n",
    "# First covariant derivative\n",
    "# D_{j} A_{i} = A_{i,j} - \\Gamma^{k}_{ij} A_{k}\n",
    "AD_dHatD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        AD_dHatD[i][j] = AD_dD[i][j]\n",
    "        for k in range(DIM):\n",
    "            AD_dHatD[i][j] -= rfm.GammahatUDD[k][i][j] * AD[k]\n",
    "\n",
    "# Second covariant derivative\n",
    "# D_{k} D_{j} A_{i} = \\partial_{k} D_{j} A_{i} - \\Gamma^{l}_{jk} D_{l} A_{i}\n",
    "#                    - \\Gamma^{l}_{ik} D_{j} A_{l}\n",
    "#                   =  A_{i,jk}\n",
    "#                    - \\Gamma^{l}_{ij,k} A_{l}\n",
    "#                    - \\Gamma^{l}_{ij}   A_{l,k}\n",
    "#                    - \\Gamma^{l}_{jk}   A_{i;\\hat{l}}\n",
    "#                    - \\Gamma^{l}_{ik}   A_{l;\\hat{j}}\n",
    "AD_dHatDD = ixp.zerorank3()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        for k in range(DIM):\n",
    "            AD_dHatDD[i][j][k] = AD_dDD[i][j][k]\n",
    "            for l in range(DIM):\n",
    "                AD_dHatDD[i][j][k] += - rfm.GammahatUDDdD[l][i][j][k]  * AD[l]  \\\n",
    "                                      - rfm.GammahatUDD[l][i][j] *    AD_dD[l][k] \\\n",
    "                                      - rfm.GammahatUDD[l][j][k] * AD_dHatD[i][l] \\\n",
    "                                      - rfm.GammahatUDD[l][i][k] * AD_dHatD[l][j]\n",
    "\n",
    "# Covariant divergence\n",
    "# D_{i} A^{i} = ghat^{ij} D_{j} A_{i}\n",
    "DivA = 0\n",
    "# Gradient of covariant divergence\n",
    "# DivA_dD_{i} = ghat^{jk} A_{k;\\hat{j}\\hat{i}}\n",
    "DivA_dD = ixp.zerorank1()\n",
    "# Covariant Laplacian\n",
    "# LapAD_{i} = ghat^{jk} A_{i;\\hat{j}\\hat{k}}\n",
    "LapAD = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        DivA += rfm.ghatUU[i][j] * AD_dHatD[i][j]\n",
    "        for k in range(DIM):\n",
    "            DivA_dD[i] += rfm.ghatUU[j][k] * AD_dHatDD[k][j][i]\n",
    "            LapAD[i]   += rfm.ghatUU[j][k] * AD_dHatDD[i][j][k]\n",
    "\n",
    "# Step 5: Define right-hand sides for the evolution.\n",
    "AD_rhs = ixp.zerorank1()\n",
    "ED_rhs = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    AD_rhs[i] = -ED[i] - psi_dD[i]\n",
    "    ED_rhs[i] = -LapAD[i] + DivA_dD[i]\n",
    "psi_rhs = -DivA\n",
    "\n",
    "# Step 6: Generate C code for System I Maxwell's evolution equations,\n",
    "#         print output to the screen (standard out, or stdout).\n",
    "lhrh_list = []\n",
    "for i in range(DIM):\n",
    "    lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"AD\" + str(i)), rhs=AD_rhs[i]))\n",
    "    lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"ED\" + str(i)), rhs=ED_rhs[i]))\n",
    "lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"psi\"), rhs=psi_rhs))\n",
    "\n",
    "fin.FD_outputC(\"stdout\", lhrh_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='sys2'></a>\n",
    "\n",
    "# Step 2: System II \\[Back to [top](#toc)\\]\n",
    "$$\\label{sys2}$$\n",
    "\n",
    "Define the auxiliary variable\n",
    "$$\\Gamma \\equiv \\hat{D}_{i} A^{i} \\; .$$\n",
    "Substituting this into Maxwell's equations yields the system\n",
    "\\begin{align}\n",
    "\\partial_{t} A_{i} &= -E_{i} - \\hat{D}_{i} \\psi \\; , \\\\\n",
    "\\partial_{t} E_{i} &= -\\hat{D}_{j} \\hat{D}^{j} A_{i} + \\hat{D}_{i} \\Gamma \\; , \\\\\n",
    "\\partial_{t} \\psi &= -\\Gamma \\; , \\\\\n",
    "\\partial_{t} \\Gamma &= -\\hat{D}_{i} \\hat{D}^{i} \\psi \\; .\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "\n",
    "It can be shown that the Gauss constraint now satisfies the wave equation\n",
    "$$\\partial_{t}^{2} \\mathcal{C} = \\hat{D}_{i} \\hat{D}^{i} \\mathcal{C} \\; .$$\n",
    "Thus, any constraint violation introduced by numerical error propagates away at the speed of light. This property increases the stability of of the simulation, compared to System I above. A similar trick is used in the [BSSN formulation](Tutorial-BSSNCurvilinear.ipynb) of Einstein's equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:14:54.549764Z",
     "iopub.status.busy": "2021-03-07T17:14:54.513988Z",
     "iopub.status.idle": "2021-03-07T17:14:55.234470Z",
     "shell.execute_reply": "2021-03-07T17:14:55.233887Z"
    }
   },
   "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 AD_dD00 = invdx0*(-2*AD0_i0m1_i1_i2/3 + AD0_i0m2_i1_i2/12 + 2*AD0_i0p1_i1_i2/3 - AD0_i0p2_i1_i2/12),\n",
      "    *    const double AD_dD01 = invdx1*(-2*AD0_i0_i1m1_i2/3 + AD0_i0_i1m2_i2/12 + 2*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),\n",
      "    *    const double AD_dD02 = invdx2*(-2*AD0_i0_i1_i2m1/3 + AD0_i0_i1_i2m2/12 + 2*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),\n",
      "    *    const double AD_dD10 = invdx0*(-2*AD1_i0m1_i1_i2/3 + AD1_i0m2_i1_i2/12 + 2*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),\n",
      "    *    const double AD_dD11 = invdx1*(-2*AD1_i0_i1m1_i2/3 + AD1_i0_i1m2_i2/12 + 2*AD1_i0_i1p1_i2/3 - AD1_i0_i1p2_i2/12),\n",
      "    *    const double AD_dD12 = invdx2*(-2*AD1_i0_i1_i2m1/3 + AD1_i0_i1_i2m2/12 + 2*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),\n",
      "    *    const double AD_dD20 = invdx0*(-2*AD2_i0m1_i1_i2/3 + AD2_i0m2_i1_i2/12 + 2*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),\n",
      "    *    const double AD_dD21 = invdx1*(-2*AD2_i0_i1m1_i2/3 + AD2_i0_i1m2_i2/12 + 2*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),\n",
      "    *    const double AD_dD22 = invdx2*(-2*AD2_i0_i1_i2m1/3 + AD2_i0_i1_i2m2/12 + 2*AD2_i0_i1_i2p1/3 - AD2_i0_i1_i2p2/12),\n",
      "    *    const double AD_dDD000 = invdx0**2*(-5*AD0/2 + 4*AD0_i0m1_i1_i2/3 - AD0_i0m2_i1_i2/12 + 4*AD0_i0p1_i1_i2/3 - AD0_i0p2_i1_i2/12),\n",
      "    *    const double AD_dDD011 = invdx1**2*(-5*AD0/2 + 4*AD0_i0_i1m1_i2/3 - AD0_i0_i1m2_i2/12 + 4*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),\n",
      "    *    const double AD_dDD022 = invdx2**2*(-5*AD0/2 + 4*AD0_i0_i1_i2m1/3 - AD0_i0_i1_i2m2/12 + 4*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),\n",
      "    *    const double AD_dDD100 = invdx0**2*(-5*AD1/2 + 4*AD1_i0m1_i1_i2/3 - AD1_i0m2_i1_i2/12 + 4*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),\n",
      "    *    const double AD_dDD111 = invdx1**2*(-5*AD1/2 + 4*AD1_i0_i1m1_i2/3 - AD1_i0_i1m2_i2/12 + 4*AD1_i0_i1p1_i2/3 - AD1_i0_i1p2_i2/12),\n",
      "    *    const double AD_dDD122 = invdx2**2*(-5*AD1/2 + 4*AD1_i0_i1_i2m1/3 - AD1_i0_i1_i2m2/12 + 4*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),\n",
      "    *    const double AD_dDD200 = invdx0**2*(-5*AD2/2 + 4*AD2_i0m1_i1_i2/3 - AD2_i0m2_i1_i2/12 + 4*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),\n",
      "    *    const double AD_dDD211 = invdx1**2*(-5*AD2/2 + 4*AD2_i0_i1m1_i2/3 - AD2_i0_i1m2_i2/12 + 4*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),\n",
      "    *    const double AD_dDD222 = invdx2**2*(-5*AD2/2 + 4*AD2_i0_i1_i2m1/3 - AD2_i0_i1_i2m2/12 + 4*AD2_i0_i1_i2p1/3 - AD2_i0_i1_i2p2/12),\n",
      "    *    const double Gamma_dD0 = invdx0*(-2*Gamma_i0m1_i1_i2/3 + Gamma_i0m2_i1_i2/12 + 2*Gamma_i0p1_i1_i2/3 - Gamma_i0p2_i1_i2/12),\n",
      "    *    const double Gamma_dD1 = invdx1*(-2*Gamma_i0_i1m1_i2/3 + Gamma_i0_i1m2_i2/12 + 2*Gamma_i0_i1p1_i2/3 - Gamma_i0_i1p2_i2/12),\n",
      "    *    const double Gamma_dD2 = invdx2*(-2*Gamma_i0_i1_i2m1/3 + Gamma_i0_i1_i2m2/12 + 2*Gamma_i0_i1_i2p1/3 - Gamma_i0_i1_i2p2/12),\n",
      "    *    const double psi_dD0 = invdx0*(-2*psi_i0m1_i1_i2/3 + psi_i0m2_i1_i2/12 + 2*psi_i0p1_i1_i2/3 - psi_i0p2_i1_i2/12),\n",
      "    *    const double psi_dD1 = invdx1*(-2*psi_i0_i1m1_i2/3 + psi_i0_i1m2_i2/12 + 2*psi_i0_i1p1_i2/3 - psi_i0_i1p2_i2/12),\n",
      "    *    const double psi_dD2 = invdx2*(-2*psi_i0_i1_i2m1/3 + psi_i0_i1_i2m2/12 + 2*psi_i0_i1_i2p1/3 - psi_i0_i1_i2p2/12),\n",
      "    *    const double psi_dDD00 = invdx0**2*(-5*psi/2 + 4*psi_i0m1_i1_i2/3 - psi_i0m2_i1_i2/12 + 4*psi_i0p1_i1_i2/3 - psi_i0p2_i1_i2/12),\n",
      "    *    const double psi_dDD11 = invdx1**2*(-5*psi/2 + 4*psi_i0_i1m1_i2/3 - psi_i0_i1m2_i2/12 + 4*psi_i0_i1p1_i2/3 - psi_i0_i1p2_i2/12),\n",
      "    *    const double psi_dDD22 = invdx2**2*(-5*psi/2 + 4*psi_i0_i1_i2m1/3 - psi_i0_i1_i2m2/12 + 4*psi_i0_i1_i2p1/3 - psi_i0_i1_i2p2/12)]\"\n",
      "    */\n",
      "   const double psi_i0_i1_i2m2 = in_gfs[IDX4(PSIGF, i0,i1,i2-2)];\n",
      "   const double psi_i0_i1_i2m1 = in_gfs[IDX4(PSIGF, i0,i1,i2-1)];\n",
      "   const double psi_i0_i1m2_i2 = in_gfs[IDX4(PSIGF, i0,i1-2,i2)];\n",
      "   const double psi_i0_i1m1_i2 = in_gfs[IDX4(PSIGF, i0,i1-1,i2)];\n",
      "   const double psi_i0m2_i1_i2 = in_gfs[IDX4(PSIGF, i0-2,i1,i2)];\n",
      "   const double psi_i0m1_i1_i2 = in_gfs[IDX4(PSIGF, i0-1,i1,i2)];\n",
      "   const double psi = in_gfs[IDX4(PSIGF, i0,i1,i2)];\n",
      "   const double psi_i0p1_i1_i2 = in_gfs[IDX4(PSIGF, i0+1,i1,i2)];\n",
      "   const double psi_i0p2_i1_i2 = in_gfs[IDX4(PSIGF, i0+2,i1,i2)];\n",
      "   const double psi_i0_i1p1_i2 = in_gfs[IDX4(PSIGF, i0,i1+1,i2)];\n",
      "   const double psi_i0_i1p2_i2 = in_gfs[IDX4(PSIGF, i0,i1+2,i2)];\n",
      "   const double psi_i0_i1_i2p1 = in_gfs[IDX4(PSIGF, i0,i1,i2+1)];\n",
      "   const double psi_i0_i1_i2p2 = in_gfs[IDX4(PSIGF, i0,i1,i2+2)];\n",
      "   const double ED0 = in_gfs[IDX4(ED0GF, i0,i1,i2)];\n",
      "   const double ED1 = in_gfs[IDX4(ED1GF, i0,i1,i2)];\n",
      "   const double ED2 = in_gfs[IDX4(ED2GF, i0,i1,i2)];\n",
      "   const double AD0_i0_i1_i2m2 = in_gfs[IDX4(AD0GF, i0,i1,i2-2)];\n",
      "   const double AD0_i0_i1_i2m1 = in_gfs[IDX4(AD0GF, i0,i1,i2-1)];\n",
      "   const double AD0_i0_i1m2_i2 = in_gfs[IDX4(AD0GF, i0,i1-2,i2)];\n",
      "   const double AD0_i0_i1m1_i2 = in_gfs[IDX4(AD0GF, i0,i1-1,i2)];\n",
      "   const double AD0_i0m2_i1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2)];\n",
      "   const double AD0_i0m1_i1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2)];\n",
      "   const double AD0 = in_gfs[IDX4(AD0GF, i0,i1,i2)];\n",
      "   const double AD0_i0p1_i1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2)];\n",
      "   const double AD0_i0p2_i1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2)];\n",
      "   const double AD0_i0_i1p1_i2 = in_gfs[IDX4(AD0GF, i0,i1+1,i2)];\n",
      "   const double AD0_i0_i1p2_i2 = in_gfs[IDX4(AD0GF, i0,i1+2,i2)];\n",
      "   const double AD0_i0_i1_i2p1 = in_gfs[IDX4(AD0GF, i0,i1,i2+1)];\n",
      "   const double AD0_i0_i1_i2p2 = in_gfs[IDX4(AD0GF, i0,i1,i2+2)];\n",
      "   const double AD1_i0_i1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1,i2-2)];\n",
      "   const double AD1_i0_i1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1,i2-1)];\n",
      "   const double AD1_i0_i1m2_i2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2)];\n",
      "   const double AD1_i0_i1m1_i2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2)];\n",
      "   const double AD1_i0m2_i1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1,i2)];\n",
      "   const double AD1_i0m1_i1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1,i2)];\n",
      "   const double AD1 = in_gfs[IDX4(AD1GF, i0,i1,i2)];\n",
      "   const double AD1_i0p1_i1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1,i2)];\n",
      "   const double AD1_i0p2_i1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1,i2)];\n",
      "   const double AD1_i0_i1p1_i2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2)];\n",
      "   const double AD1_i0_i1p2_i2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2)];\n",
      "   const double AD1_i0_i1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1,i2+1)];\n",
      "   const double AD1_i0_i1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1,i2+2)];\n",
      "   const double AD2_i0_i1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1,i2-2)];\n",
      "   const double AD2_i0_i1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1,i2-1)];\n",
      "   const double AD2_i0_i1m2_i2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2)];\n",
      "   const double AD2_i0_i1m1_i2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2)];\n",
      "   const double AD2_i0m2_i1_i2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2)];\n",
      "   const double AD2_i0m1_i1_i2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2)];\n",
      "   const double AD2 = in_gfs[IDX4(AD2GF, i0,i1,i2)];\n",
      "   const double AD2_i0p1_i1_i2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2)];\n",
      "   const double AD2_i0p2_i1_i2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2)];\n",
      "   const double AD2_i0_i1p1_i2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2)];\n",
      "   const double AD2_i0_i1p2_i2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2)];\n",
      "   const double AD2_i0_i1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1,i2+1)];\n",
      "   const double AD2_i0_i1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1,i2+2)];\n",
      "   const double Gamma_i0_i1_i2m2 = in_gfs[IDX4(GAMMAGF, i0,i1,i2-2)];\n",
      "   const double Gamma_i0_i1_i2m1 = in_gfs[IDX4(GAMMAGF, i0,i1,i2-1)];\n",
      "   const double Gamma_i0_i1m2_i2 = in_gfs[IDX4(GAMMAGF, i0,i1-2,i2)];\n",
      "   const double Gamma_i0_i1m1_i2 = in_gfs[IDX4(GAMMAGF, i0,i1-1,i2)];\n",
      "   const double Gamma_i0m2_i1_i2 = in_gfs[IDX4(GAMMAGF, i0-2,i1,i2)];\n",
      "   const double Gamma_i0m1_i1_i2 = in_gfs[IDX4(GAMMAGF, i0-1,i1,i2)];\n",
      "   const double Gamma = in_gfs[IDX4(GAMMAGF, i0,i1,i2)];\n",
      "   const double Gamma_i0p1_i1_i2 = in_gfs[IDX4(GAMMAGF, i0+1,i1,i2)];\n",
      "   const double Gamma_i0p2_i1_i2 = in_gfs[IDX4(GAMMAGF, i0+2,i1,i2)];\n",
      "   const double Gamma_i0_i1p1_i2 = in_gfs[IDX4(GAMMAGF, i0,i1+1,i2)];\n",
      "   const double Gamma_i0_i1p2_i2 = in_gfs[IDX4(GAMMAGF, i0,i1+2,i2)];\n",
      "   const double Gamma_i0_i1_i2p1 = in_gfs[IDX4(GAMMAGF, i0,i1,i2+1)];\n",
      "   const double Gamma_i0_i1_i2p2 = in_gfs[IDX4(GAMMAGF, i0,i1,i2+2)];\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 = -AD0_i0_i1p2_i2;\n",
      "   const double FDPart1_9 = ((invdx0)*(invdx0));\n",
      "   const double FDPart1_10 = -AD0*FDPart1_Rational_5_2;\n",
      "   const double FDPart1_11 = ((invdx1)*(invdx1));\n",
      "   const double FDPart1_12 = ((invdx2)*(invdx2));\n",
      "   const double FDPart1_13 = -AD1*FDPart1_Rational_5_2;\n",
      "   const double FDPart1_14 = -AD2*FDPart1_Rational_5_2;\n",
      "   const double FDPart1_18 = -FDPart1_Rational_5_2*psi;\n",
      "   const double AD_dD00 = invdx0*(FDPart1_Rational_1_12*(AD0_i0m2_i1_i2 - AD0_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD0_i0m1_i1_i2 + AD0_i0p1_i1_i2));\n",
      "   const double AD_dD01 = invdx1*(FDPart1_Rational_1_12*(AD0_i0_i1m2_i2 + FDPart1_1) + FDPart1_Rational_2_3*(-AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));\n",
      "   const double AD_dD02 = invdx2*(FDPart1_Rational_1_12*(AD0_i0_i1_i2m2 - AD0_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));\n",
      "   const double AD_dD10 = invdx0*(FDPart1_Rational_1_12*(AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));\n",
      "   const double AD_dD11 = invdx1*(FDPart1_Rational_1_12*(AD1_i0_i1m2_i2 - AD1_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD1_i0_i1m1_i2 + AD1_i0_i1p1_i2));\n",
      "   const double AD_dD12 = invdx2*(FDPart1_Rational_1_12*(AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));\n",
      "   const double AD_dD20 = invdx0*(FDPart1_Rational_1_12*(AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));\n",
      "   const double AD_dD21 = invdx1*(FDPart1_Rational_1_12*(AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));\n",
      "   const double AD_dD22 = invdx2*(FDPart1_Rational_1_12*(AD2_i0_i1_i2m2 - AD2_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD2_i0_i1_i2m1 + AD2_i0_i1_i2p1));\n",
      "   const double AD_dDD000 = FDPart1_9*(FDPart1_10 + FDPart1_Rational_1_12*(-AD0_i0m2_i1_i2 - AD0_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD0_i0m1_i1_i2 + AD0_i0p1_i1_i2));\n",
      "   const double AD_dDD011 = FDPart1_11*(FDPart1_10 + FDPart1_Rational_1_12*(-AD0_i0_i1m2_i2 + FDPart1_1) + FDPart1_Rational_4_3*(AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));\n",
      "   const double AD_dDD022 = FDPart1_12*(FDPart1_10 + FDPart1_Rational_1_12*(-AD0_i0_i1_i2m2 - AD0_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));\n",
      "   const double AD_dDD100 = FDPart1_9*(FDPart1_13 + FDPart1_Rational_1_12*(-AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));\n",
      "   const double AD_dDD111 = FDPart1_11*(FDPart1_13 + FDPart1_Rational_1_12*(-AD1_i0_i1m2_i2 - AD1_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD1_i0_i1m1_i2 + AD1_i0_i1p1_i2));\n",
      "   const double AD_dDD122 = FDPart1_12*(FDPart1_13 + FDPart1_Rational_1_12*(-AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));\n",
      "   const double AD_dDD200 = FDPart1_9*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));\n",
      "   const double AD_dDD211 = FDPart1_11*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));\n",
      "   const double AD_dDD222 = FDPart1_12*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0_i1_i2m2 - AD2_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD2_i0_i1_i2m1 + AD2_i0_i1_i2p1));\n",
      "   const double Gamma_dD0 = invdx0*(FDPart1_Rational_1_12*(Gamma_i0m2_i1_i2 - Gamma_i0p2_i1_i2) + FDPart1_Rational_2_3*(-Gamma_i0m1_i1_i2 + Gamma_i0p1_i1_i2));\n",
      "   const double Gamma_dD1 = invdx1*(FDPart1_Rational_1_12*(Gamma_i0_i1m2_i2 - Gamma_i0_i1p2_i2) + FDPart1_Rational_2_3*(-Gamma_i0_i1m1_i2 + Gamma_i0_i1p1_i2));\n",
      "   const double Gamma_dD2 = invdx2*(FDPart1_Rational_1_12*(Gamma_i0_i1_i2m2 - Gamma_i0_i1_i2p2) + FDPart1_Rational_2_3*(-Gamma_i0_i1_i2m1 + Gamma_i0_i1_i2p1));\n",
      "   const double psi_dD0 = invdx0*(FDPart1_Rational_1_12*(psi_i0m2_i1_i2 - psi_i0p2_i1_i2) + FDPart1_Rational_2_3*(-psi_i0m1_i1_i2 + psi_i0p1_i1_i2));\n",
      "   const double psi_dD1 = invdx1*(FDPart1_Rational_1_12*(psi_i0_i1m2_i2 - psi_i0_i1p2_i2) + FDPart1_Rational_2_3*(-psi_i0_i1m1_i2 + psi_i0_i1p1_i2));\n",
      "   const double psi_dD2 = invdx2*(FDPart1_Rational_1_12*(psi_i0_i1_i2m2 - psi_i0_i1_i2p2) + FDPart1_Rational_2_3*(-psi_i0_i1_i2m1 + psi_i0_i1_i2p1));\n",
      "   const double psi_dDD00 = FDPart1_9*(FDPart1_18 + FDPart1_Rational_1_12*(-psi_i0m2_i1_i2 - psi_i0p2_i1_i2) + FDPart1_Rational_4_3*(psi_i0m1_i1_i2 + psi_i0p1_i1_i2));\n",
      "   const double psi_dDD11 = FDPart1_11*(FDPart1_18 + FDPart1_Rational_1_12*(-psi_i0_i1m2_i2 - psi_i0_i1p2_i2) + FDPart1_Rational_4_3*(psi_i0_i1m1_i2 + psi_i0_i1p1_i2));\n",
      "   const double psi_dDD22 = FDPart1_12*(FDPart1_18 + FDPart1_Rational_1_12*(-psi_i0_i1_i2m2 - psi_i0_i1_i2p2) + FDPart1_Rational_4_3*(psi_i0_i1_i2m1 + psi_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(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0,\n",
      "    *    rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = -AD_dDD000 + Gamma_dD0 - (AD_dD00*xx0 - AD_dD11/xx0 + AD_dDD011 - (AD0*xx0 + AD_dD11)/xx0)/xx0**2 - (AD_dD00*xx0*sin(xx1)**2 - AD_dD22/xx0 + AD_dDD022 + (-AD1/xx0 + AD_dD01)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/xx0)/(xx0**2*sin(xx1)**2),\n",
      "    *    rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1,\n",
      "    *    rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1/xx0**2 + AD_dD10/xx0 - AD_dDD100 + Gamma_dD1 + (-AD1/xx0 + AD_dD10)/xx0 - (AD_dD01*xx0 + AD_dDD111 + xx0*(-AD1/xx0 + AD_dD01) + xx0*(-AD1/xx0 + AD_dD10))/xx0**2 - (-AD_dD22*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD122 + xx0*(-AD1/xx0 + AD_dD10)*sin(xx1)**2 + (AD0*xx0 + AD_dD11)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)*sin(2*xx1)/(2*sin(xx1)**2))/(xx0**2*sin(xx1)**2),\n",
      "    *    rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2,\n",
      "    *    rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2/xx0**2 + AD_dD20/xx0 - AD_dDD200 + Gamma_dD2 + (-AD2/xx0 + AD_dD20)/xx0 - (AD2*(-cos(2*xx1)/sin(xx1)**2 + sin(2*xx1)*cos(xx1)/sin(xx1)**3) - AD_dD21*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD211 + xx0*(-AD2/xx0 + AD_dD20) - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/(2*sin(xx1)**2))/xx0**2 - (AD_dD02*xx0*sin(xx1)**2 + AD_dD12*sin(2*xx1)/2 + AD_dDD222 + xx0*(-AD2/xx0 + AD_dD02)*sin(xx1)**2 + xx0*(-AD2/xx0 + AD_dD20)*sin(xx1)**2 + (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD12)*sin(2*xx1)/2 + (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/2)/(xx0**2*sin(xx1)**2),\n",
      "    *    rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -Gamma,\n",
      "    *    rhs_gfs[IDX4(GAMMAGF, i0, i1, i2)] = -psi_dDD00 - (psi_dD0*xx0 + psi_dDD11)/xx0**2 - (psi_dD0*xx0*sin(xx1)**2 + psi_dD1*sin(2*xx1)/2 + psi_dDD22)/(xx0**2*sin(xx1)**2)]\"\n",
      "    */\n",
      "   const double FDPart3_0 = (1.0/((xx0)*(xx0)));\n",
      "   const double FDPart3_2 = (1.0/(xx0));\n",
      "   const double FDPart3_4 = AD0*xx0 + AD_dD11;\n",
      "   const double FDPart3_5 = sin(xx1);\n",
      "   const double FDPart3_6 = ((FDPart3_5)*(FDPart3_5));\n",
      "   const double FDPart3_7 = -AD1*FDPart3_2;\n",
      "   const double FDPart3_10 = sin(2*xx1);\n",
      "   const double FDPart3_11 = (1.0/2.0)*FDPart3_10;\n",
      "   const double FDPart3_12 = AD0*FDPart3_6*xx0 + AD1*FDPart3_11 + AD_dD22;\n",
      "   const double FDPart3_13 = (1.0/(FDPart3_6));\n",
      "   const double FDPart3_14 = FDPart3_0*FDPart3_13;\n",
      "   const double FDPart3_16 = xx0*(AD_dD10 + FDPart3_7);\n",
      "   const double FDPart3_17 = FDPart3_11*FDPart3_13;\n",
      "   const double FDPart3_18 = -AD2*FDPart3_2;\n",
      "   const double FDPart3_20 = xx0*(AD_dD20 + FDPart3_18);\n",
      "   const double FDPart3_21 = -AD2*FDPart3_17;\n",
      "   const double FDPart3_22 = FDPart3_11*(AD_dD21 + FDPart3_21);\n",
      "   rhs_gfs[IDX4(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0;\n",
      "   rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = -AD_dDD000 - FDPart3_0*(AD_dD00*xx0 - AD_dD11*FDPart3_2 + AD_dDD011 - FDPart3_2*FDPart3_4) - FDPart3_14*(AD_dD00*FDPart3_6*xx0 - AD_dD22*FDPart3_2 + AD_dDD022 + FDPart3_11*(AD_dD01 + FDPart3_7) - FDPart3_12*FDPart3_2) + Gamma_dD0;\n",
      "   rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1;\n",
      "   rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1*FDPart3_0 + AD_dD10*FDPart3_2 - AD_dDD100 - FDPart3_0*(AD_dD01*xx0 + AD_dDD111 + FDPart3_16 + xx0*(AD_dD01 + FDPart3_7)) - FDPart3_14*(-AD_dD22*FDPart3_17 + AD_dDD122 + FDPart3_11*FDPart3_4 - FDPart3_12*FDPart3_17 + FDPart3_16*FDPart3_6) + FDPart3_2*(AD_dD10 + FDPart3_7) + Gamma_dD1;\n",
      "   rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2;\n",
      "   rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2*FDPart3_0 + AD_dD20*FDPart3_2 - AD_dDD200 - FDPart3_0*(AD2*(FDPart3_10*cos(xx1)/((FDPart3_5)*(FDPart3_5)*(FDPart3_5)) - FDPart3_13*cos(2*xx1)) - AD_dD21*FDPart3_17 + AD_dDD211 - FDPart3_13*FDPart3_22 + FDPart3_20) - FDPart3_14*(AD_dD02*FDPart3_6*xx0 + AD_dD12*FDPart3_11 + AD_dDD222 + FDPart3_11*(AD_dD12 + FDPart3_21) + FDPart3_20*FDPart3_6 + FDPart3_22 + FDPart3_6*xx0*(AD_dD02 + FDPart3_18)) + FDPart3_2*(AD_dD20 + FDPart3_18) + Gamma_dD2;\n",
      "   rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -Gamma;\n",
      "   rhs_gfs[IDX4(GAMMAGF, i0, i1, i2)] = -FDPart3_0*(psi_dD0*xx0 + psi_dDD11) - FDPart3_14*(FDPart3_11*psi_dD1 + FDPart3_6*psi_dD0*xx0 + psi_dDD22) - psi_dDD00;\n",
      "}\n"
     ]
    }
   ],
   "source": [
    "# We inherit here all of the definitions from System I, above\n",
    "\n",
    "# Step 7a: Register the scalar auxiliary variable \\Gamma\n",
    "Gamma = gri.register_gridfunctions(\"EVOL\", [\"Gamma\"])\n",
    "\n",
    "# Step 7b: Declare the ordinary gradient \\partial_{i} \\Gamma\n",
    "Gamma_dD = ixp.declarerank1(\"Gamma_dD\")\n",
    "\n",
    "# Step 8a: Construct the second covariant derivative of the scalar \\psi\n",
    "# \\psi_{;\\hat{i}\\hat{j}} = \\psi_{,i;\\hat{j}}\n",
    "#                        = \\psi_{,ij} - \\Gamma^{k}_{ij} \\psi_{,k}\n",
    "psi_dDD = ixp.declarerank2(\"psi_dDD\", \"sym01\")\n",
    "psi_dHatDD = ixp.zerorank2()\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        psi_dHatDD[i][j] = psi_dDD[i][j]\n",
    "        for k in range(DIM):\n",
    "            psi_dHatDD[i][j] += - rfm.GammahatUDD[k][i][j] * psi_dD[k]\n",
    "\n",
    "# Step 8b: Construct the covariant Laplacian of \\psi\n",
    "# Lappsi = ghat^{ij} D_{j} D_{i} \\psi\n",
    "Lappsi = 0\n",
    "for i in range(DIM):\n",
    "    for j in range(DIM):\n",
    "        Lappsi += rfm.ghatUU[i][j] * psi_dHatDD[i][j]\n",
    "\n",
    "# Step 9: Define right-hand sides for the evolution.\n",
    "AD_rhs = ixp.zerorank1()\n",
    "ED_rhs = ixp.zerorank1()\n",
    "for i in range(DIM):\n",
    "    AD_rhs[i] = -ED[i] - psi_dD[i]\n",
    "    ED_rhs[i] = -LapAD[i] + Gamma_dD[i]\n",
    "psi_rhs = -Gamma\n",
    "Gamma_rhs = -Lappsi\n",
    "\n",
    "# Step 10: Generate C code for System II Maxwell's evolution equations,\n",
    "#          print output to the screen (standard out, or stdout).\n",
    "lhrh_list = []\n",
    "for i in range(DIM):\n",
    "    lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"AD\" + str(i)), rhs=AD_rhs[i]))\n",
    "    lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"ED\" + str(i)), rhs=ED_rhs[i]))\n",
    "lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"psi\"), rhs=psi_rhs))\n",
    "lhrh_list.append(lhrh(lhs=gri.gfaccess(\"rhs_gfs\", \"Gamma\"), rhs=Gamma_rhs))\n",
    "\n",
    "fin.FD_outputC(\"stdout\", lhrh_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 3: 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-MaxwellCurvilinear.pdf](Tutorial-MaxwellCurvilinear.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": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:14:55.238869Z",
     "iopub.status.busy": "2021-03-07T17:14:55.238152Z",
     "iopub.status.idle": "2021-03-07T17:14:58.527854Z",
     "shell.execute_reply": "2021-03-07T17:14:58.527060Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-MaxwellCurvilinear.tex, and compiled LaTeX file to PDF\n",
      "    file Tutorial-MaxwellCurvilinear.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-MaxwellCurvilinear\")"
   ]
  }
 ],
 "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
}