in_progress/Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-Afield_flux.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",
    "# Start-to-Finish Example: Unit Testing `GiRaFFE_NRPy`: Induction Equation\n",
    "\n",
    "## Author: Patrick Nelson\n",
    "\n",
    "## This module validates the routine to compute the flux of $\\epsilon_{ijk} v^j B^k$ for `GiRaFFE`.\n",
    "\n",
    "**Notebook Status:** <font color='green'><b>Validated</b></font>\n",
    "\n",
    "**Validation Notes:** This module validates the routines in [Tutorial-GiRaFFE_NRPy-Afield_flux_handwritten](Tutorial-GiRaFFE_NRPy-Afield_flux_handwritten.ipynb).\n",
    "\n",
    "### NRPy+ Source Code for this module: \n",
    "* [GiRaFFE_NRPy/GiRaFFE_NRPy_Afield_flux_handwritten.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_Afield_flux_handwritten.py) [\\[**tutorial**\\]](Tutorial-GiRaFFE_NRPy-Afield_flux_handwritten.ipynb) Generates the symbolic expressions needed to compute the flux term of the induction equation right-hand side.\n",
    "\n",
    "## Introduction:\n",
    "\n",
    "This notebook validates our algorithm to compute the flux of $\\epsilon_{ijk} v^j B^k$ through cell faces, which contributes to the right-hand side of the evolution equation for $A_i$, for use in `GiRaFFE_NRPy`. Because the original `GiRaFFE` used staggered grids and we do not, we can not trivially do a direct comparison to the old code. Instead, we will compare the numerical results with the expected analytic results. \n",
    "\n",
    "It is, in general, good coding practice to unit test functions individually to verify that they produce the expected and intended output. Here, we expect our functions to produce the correct cross product between the velocity and magnetic field in an arbitrary spacetime. To that end, we will choose functions with relatively simple analytic forms.\n",
    "\n",
    "<font color='yellow'><b>TERRENCE SAYS: wording is confusing here</b></font>\n",
    "\n",
    "In this test, we will generate analytic forms for the magnetic field and three-velocity. We will need to do this in a $7 \\times 7 \\times 7$ cube in order to run the PPM routine on the data to generate the inputs compute the non-gauge terms of the RHS of the induction equation, since it would be impractical to usefully spoof the left- and right-face values for the HLLE solver. We care about this here, because we are comparing against an analytic expression and not the old code.\n",
    "\n",
    "When this notebook is run, the difference between the approximate and exact right-hand sides will be output to text files that can be found in the same directory as this notebook. These will be read in in [Step 3](#convergence), and used there to confirm second order convergence of the algorithm. "
   ]
  },
  {
   "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](#setup): Set up core functions and parameters for unit testing the A2B algorithm\n",
    "    1. [Step 1.a](#metric) Set analytic 3-metric data\n",
    "    1. [Step 1.b](#magnetic) Set analytic magnetic field\n",
    "    1. [Step 1.c](#velocity) Set analytic Valencia three-velocity\n",
    "    1. [Step 1.d](#functions) Generate C functions to write the test data\n",
    "    1. [Step 1.e](#free_parameters) Set free parameters in the code\n",
    "    1. [Step 1.f](#module) Generate `GiRaFFE_NRPy` Files\n",
    "    1. [Step 1.g](#electric_field) Calculate the $E_i$ field contribution to $\\partial_t A_i$\n",
    "    1. [Step 1.h](#exact_flux) Calculate the *exact* flux of $\\epsilon_{ijk} v^j B^k$\n",
    "1. [Step 2](#mainc): `Afield_flux_unit_test.c`: The Main C Code\n",
    "    1. [Step 2.a](#compile_run): Compile and run the code\n",
    "1. [Step 3](#convergence): Code validation: Verify that relative error in numerical solution converges to zero at the expected order\n",
    "1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='setup'></a>\n",
    "\n",
    "# Step 1: Set up core functions and parameters for unit testing the A2B algorithm \\[Back to [top](#toc)\\]\n",
    "\n",
    "$$\\label{setup}$$\n",
    "\n",
    "We'll start by appending the relevant paths to `sys.path` so that we can access sympy modules in other places. Then, we'll import NRPy+ core functionality and set up a directory in which to carry out our test. We must also set the desired finite differencing order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.040703Z",
     "iopub.status.busy": "2020-12-17T02:17:01.040227Z",
     "iopub.status.idle": "2020-12-17T02:17:01.223614Z",
     "shell.execute_reply": "2020-12-17T02:17:01.223072Z"
    }
   },
   "outputs": [],
   "source": [
    "import os, sys, shutil        # Standard Python modules for multiplatform OS-level functions\n",
    "# First, we'll add the parent directory to the list of directories Python will check for modules.\n",
    "nrpy_dir_path = os.path.join(\"..\")\n",
    "if nrpy_dir_path not in sys.path:\n",
    "    sys.path.append(nrpy_dir_path)\n",
    "nrpy_dir_path = os.path.join(\"..\",\"..\")\n",
    "if nrpy_dir_path not in sys.path:\n",
    "    sys.path.append(nrpy_dir_path)\n",
    "\n",
    "from outputC import outCfunction, outputC, lhrh # NRPy+: Core C code output module\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import finite_difference as fin  # NRPy+: Finite difference C code generation module\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import grid as gri               # NRPy+: Functions having to do with 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",
    "import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface\n",
    "\n",
    "Ccodesdir = \"Start-to-Finish-UnitTests/Afield_flux_UnitTest/\"\n",
    "\n",
    "# First remove C code output directory if it exists\n",
    "# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
    "shutil.rmtree(Ccodesdir, ignore_errors=True)\n",
    "# Then create a fresh directory\n",
    "cmd.mkdir(Ccodesdir)\n",
    "\n",
    "outdir = os.path.join(Ccodesdir,\"output/\")\n",
    "cmd.mkdir(outdir)\n",
    "\n",
    "thismodule = \"Start_to_Finish_UnitTest-GiRaFFE_NRPy-Afield_flux\"\n",
    "Use_Shock_Data = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='metric'></a>\n",
    "\n",
    "## Step 1.a: Set analytic 3-metric data \\[Back to [top](#toc)\\]\n",
    "$$\\label{metric}$$\n",
    "\n",
    "First, we'll set some arbitrary functions for the three-metric $\\gamma_{ij}$, the shift vector $\\beta^i$, and the lapse function $\\alpha$. We will need these to compute the characteristic speeds in each direction. In addition, we will need the metric determinant to exactly compute the electric field at the test point.\n",
    "\n",
    "We will use this basic form for $\\gamma_{ij}$. We will thus write a function with the following arbitrary equations. \n",
    "\n",
    "\\begin{align}\n",
    "\\gamma_{xx} &= ax^3 + by^3 + cz^3 + dy^2 + ez^2 + 1 \\\\\n",
    "\\gamma_{yy} &= gx^3 + hy^3 + lz^3 + mx^2 + nz^2 + 1 \\\\\n",
    "\\gamma_{zz} &= px^3 + qy^3 + rz^3 + sx^2 + ty^2 + 1. \\\\\n",
    "\\gamma_{xy} &= \\frac{1}{10} \\exp\\left(-\\left((x-b)^2+(y-c)^2+(z-d)^2\\right)\\right) \\\\\n",
    "\\gamma_{xz} &= \\frac{1}{10} \\exp\\left(-\\left((x-g)^2+(y-h)^2+(z-l)^2\\right)\\right) \\\\\n",
    "\\gamma_{yz} &= \\frac{1}{10} \\exp\\left(-\\left((x-n)^2+(y-o)^2+(z-p)^2\\right)\\right), \\\\\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.236343Z",
     "iopub.status.busy": "2020-12-17T02:17:01.235780Z",
     "iopub.status.idle": "2020-12-17T02:17:01.349146Z",
     "shell.execute_reply": "2020-12-17T02:17:01.349583Z"
    }
   },
   "outputs": [],
   "source": [
    "a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t = par.Cparameters(\"REAL\",thismodule,[\"a\",\"b\",\"c\",\"d\",\"e\",\"f\",\"g\",\"h\",\"l\",\"m\",\"n\",\"o\",\"p\",\"q\",\"r\",\"s\",\"t\"],1e300)\n",
    "# Note that the above default value allows us to set these directly in the C code\n",
    "M_PI  = par.Cparameters(\"#define\",thismodule,[\"M_PI\"], \"\")\n",
    "\n",
    "\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n",
    "rfm.reference_metric()\n",
    "x = rfm.xx_to_Cart[0]\n",
    "y = rfm.xx_to_Cart[1]\n",
    "z = rfm.xx_to_Cart[2]\n",
    "\n",
    "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gammaDD\",\"sym01\",DIM=3)\n",
    "betaU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"betaU\",DIM=3)\n",
    "alpha = gri.register_gridfunctions(\"AUXEVOL\",\"alpha\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "\\gamma_{xx} &= ax^3 + by^3 + cz^3 + dy^2 + ez^2 + 1 \\\\\n",
    "\\gamma_{yy} &= gx^3 + hy^3 + lz^3 + mx^2 + nz^2 + 1 \\\\\n",
    "\\gamma_{zz} &= px^3 + qy^3 + rz^3 + sx^2 + ty^2 + 1. \\\\\n",
    "\\gamma_{xy} &= \\frac{1}{10} \\exp\\left(-\\left((x-b)^2+(y-c)^2+(z-d)^2\\right)\\right) \\\\\n",
    "\\gamma_{xz} &= \\frac{1}{10} \\exp\\left(-\\left((x-g)^2+(y-h)^2+(z-l)^2\\right)\\right) \\\\\n",
    "\\gamma_{yz} &= \\frac{1}{10} \\exp\\left(-\\left((x-n)^2+(y-o)^2+(z-p)^2\\right)\\right), \\\\\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.236343Z",
     "iopub.status.busy": "2020-12-17T02:17:01.235780Z",
     "iopub.status.idle": "2020-12-17T02:17:01.349146Z",
     "shell.execute_reply": "2020-12-17T02:17:01.349583Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function calculate_metric_gfs() to file Start-to-Finish-UnitTests/Afield_flux_UnitTest/calculate_metric_gfs.h\n"
     ]
    }
   ],
   "source": [
    "gammaDD[0][0] = a*x**3 + b*y**3 + c*z**3 + d*y**2 + e*z**2 + sp.sympify(1)\n",
    "gammaDD[1][1] = g*x**3 + h*y**3 + l*z**3 + m*x**2 + n*z**2 + sp.sympify(1)\n",
    "gammaDD[2][2] = p*x**3 + q*y**3 + r*z**3 + s*x**2 + t*y**2 + sp.sympify(1)\n",
    "gammaDD[0][1] = sp.Rational(1,10) * sp.exp(-((x-b)**2 + (y-c)**2 + (z-d)**2))\n",
    "gammaDD[0][2] = sp.Rational(1,10) * sp.exp(-((x-g)**2 + (y-h)**2 + (z-l)**2))\n",
    "gammaDD[1][2] = sp.Rational(1,10) * sp.exp(-((x-n)**2 + (y-o)**2 + (z-p)**2))\n",
    "\n",
    "betaU[0] = sp.sympify(0)#(sp.sympify(2)/M_PI) * sp.atan(a*x + b*y + c*z)\n",
    "betaU[1] = sp.sympify(0)#(sp.sympify(2)/M_PI) * sp.atan(b*x + c*y + a*z)\n",
    "betaU[2] = sp.sympify(0)#(sp.sympify(2)/M_PI) * sp.atan(c*x + a*y + b*z)\n",
    "\n",
    "alpha = sp.sympify(1)#sp.sympify(1) - sp.sympify(1) / (sp.sympify(2) + x**2 + y**2 + z**2)\n",
    "\n",
    "metric_gfs_to_print = [\n",
    "                       lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD00\"),rhs=gammaDD[0][0]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD01\"),rhs=gammaDD[0][1]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD02\"),rhs=gammaDD[0][2]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD11\"),rhs=gammaDD[1][1]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD12\"),rhs=gammaDD[1][2]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"gammaDD22\"),rhs=gammaDD[2][2]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"betaU0\"),rhs=betaU[0]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"betaU1\"),rhs=betaU[1]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"betaU2\"),rhs=betaU[2]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"aux_gfs\",\"alpha\"),rhs=alpha),\n",
    "                      ]\n",
    "\n",
    "desc = \"Calculate the metric gridfunctions\"\n",
    "name = \"calculate_metric_gfs\"\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"const paramstruct *restrict params,REAL *restrict xx[3],REAL *restrict auxevol_gfs\",\n",
    "    body     = fin.FD_outputC(\"returnstring\",metric_gfs_to_print,params=\"outCverbose=False\"),\n",
    "    loopopts=\"AllPoints,Read_xxs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='velocity'></a>\n",
    "\n",
    "## Step 1.b: Valencia three-velocity \\[Back to [top](#toc)\\]\n",
    "$$\\label{velocity}$$\n",
    "\n",
    "Here, we'll generate some functions for the velocity. Let's choose arctangents, since those have asymptotes that can be easily manipulated to prevent accidentally setting superluminal speeds. \n",
    "\\begin{align}\n",
    "\\bar{v}^x &= \\frac{1}{5\\pi} \\arctan(a(x-k) + by + cz) \\\\\n",
    "\\bar{v}^y &= \\frac{1}{5\\pi} \\arctan(bx + c(y-k) + az) \\\\\n",
    "\\bar{v}^z &= \\frac{1}{5\\pi} \\arctan(cx + ay + b(z-k)) \\\\\n",
    "\\end{align}\n",
    "We offset by $k$ so we have smooth data in the region of interest when `Use_Shock_Data` is false.\n",
    "\n",
    "If we want to add a jump at the origin, we can simply add $\\max(0,5x)$ to the argument of the arctangent. This will add a shock in the $x$-direction. The maximum will be described without the use of if statements using the `Min_Max_and_Piecewise_Expressions` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.363176Z",
     "iopub.status.busy": "2020-12-17T02:17:01.353086Z",
     "iopub.status.idle": "2020-12-17T02:17:01.368735Z",
     "shell.execute_reply": "2020-12-17T02:17:01.369072Z"
    }
   },
   "outputs": [],
   "source": [
    "import Min_Max_and_Piecewise_Expressions as noif\n",
    "\n",
    "offset = sp.sympify(5)\n",
    "args = ixp.zerorank1()\n",
    "args[0] = a*(x-offset) + b*y + c*z\n",
    "args[1] = b*x + c*(y-offset) + a*z\n",
    "args[2] = c*x + a*y + b*(z-offset)\n",
    "if Use_Shock_Data:\n",
    "    for i in range(3):\n",
    "        args[i] += noif.max_noif(0,5*x)\n",
    "\n",
    "ValenciavU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"ValenciavU\")\n",
    "for i in range(3):\n",
    "    ValenciavU[i] = (sp.Rational(1,5)/M_PI)*sp.atan(args[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='magnetic'></a>\n",
    "\n",
    "## Step 1.c: Magnetic field \\[Back to [top](#toc)\\]\n",
    "$$\\label{magnetic}$$\n",
    "\n",
    "We'll also need some functions for the magnetic field. Exponentials sound fun.\n",
    "\\begin{align}\n",
    "B^x &= \\exp\\left(\\frac{1}{10}(ey+fz)\\right) \\\\\n",
    "B^y &= \\exp\\left(\\frac{1}{10}(fz+dx)\\right) \\\\\n",
    "B^z &= \\exp\\left(\\frac{1}{10}(dx+ey)\\right) \\\\\n",
    "\\end{align}\n",
    "In this case, we'll add $\\max{0,5x}$ to the field to add the jump."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.376791Z",
     "iopub.status.busy": "2020-12-17T02:17:01.376260Z",
     "iopub.status.idle": "2020-12-17T02:17:01.378116Z",
     "shell.execute_reply": "2020-12-17T02:17:01.377592Z"
    }
   },
   "outputs": [],
   "source": [
    "BU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"BU\")\n",
    "BU[0] = sp.exp(sp.Rational(1,10)*(e*y+f*z))\n",
    "BU[1] = sp.exp(sp.Rational(1,10)*(f*z+d*x))\n",
    "BU[2] = sp.exp(sp.Rational(1,10)*(d*x+e*y))\n",
    "if Use_Shock_Data:\n",
    "    for i in range(3):\n",
    "        BU[i] += noif.max_noif(0,5*x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='functions'></a>\n",
    "\n",
    "## Step 1.d: Generate C functions to write the test data \\[Back to [top](#toc)\\]\n",
    "$$\\label{functions}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.398205Z",
     "iopub.status.busy": "2020-12-17T02:17:01.394874Z",
     "iopub.status.idle": "2020-12-17T02:17:01.425858Z",
     "shell.execute_reply": "2020-12-17T02:17:01.425330Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function calculate_BU() to file Start-to-Finish-UnitTests/Afield_flux_UnitTest/calculate_BU.h\n",
      "Output C function calculate_ValenciavU() to file Start-to-Finish-UnitTests/Afield_flux_UnitTest/calculate_ValenciavU.h\n"
     ]
    }
   ],
   "source": [
    "BU_to_print = [\n",
    "                lhrh(lhs=gri.gfaccess(\"out_gfs\",\"BU0\"),rhs=BU[0]),\n",
    "                lhrh(lhs=gri.gfaccess(\"out_gfs\",\"BU1\"),rhs=BU[1]),\n",
    "                lhrh(lhs=gri.gfaccess(\"out_gfs\",\"BU2\"),rhs=BU[2]),\n",
    "               ]\n",
    "\n",
    "desc = \"Calculate sample magnetic field data\"\n",
    "name = \"calculate_BU\"\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs\",\n",
    "    body     = fin.FD_outputC(\"returnstring\",BU_to_print,params=\"outCverbose=False\"),\n",
    "    loopopts=\"AllPoints,Read_xxs\")\n",
    "\n",
    "ValenciavU_to_print = [\n",
    "                       lhrh(lhs=gri.gfaccess(\"out_gfs\",\"ValenciavU0\"),rhs=ValenciavU[0]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"out_gfs\",\"ValenciavU1\"),rhs=ValenciavU[1]),\n",
    "                       lhrh(lhs=gri.gfaccess(\"out_gfs\",\"ValenciavU2\"),rhs=ValenciavU[2]),\n",
    "                      ]\n",
    "\n",
    "desc = \"Calculate sample velocity data\"\n",
    "name = \"calculate_ValenciavU\"\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs\",\n",
    "    body     = fin.FD_outputC(\"returnstring\",ValenciavU_to_print,params=\"outCverbose=False\"),\n",
    "    loopopts=\"AllPoints,Read_xxs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='free_parameters'></a>\n",
    "\n",
    "## Step 1.e: Set free parameters in the code \\[Back to [top](#toc)\\]\n",
    "$$\\label{free_parameters}$$\n",
    "\n",
    "We also need to create the files that interact with NRPy's C parameter interface. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.429081Z",
     "iopub.status.busy": "2020-12-17T02:17:01.428556Z",
     "iopub.status.idle": "2020-12-17T02:17:01.431013Z",
     "shell.execute_reply": "2020-12-17T02:17:01.430536Z"
    }
   },
   "outputs": [],
   "source": [
    "# Step 3.d.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n",
    "NGHOSTS  = par.Cparameters(\"int\",thismodule,[\"NGHOSTS\"], 3)\n",
    "# Step 3.d.ii: Set free_parameters.h\n",
    "with open(os.path.join(Ccodesdir,\"free_parameters.h\"),\"w\") as file:\n",
    "    file.write(\"\"\"\n",
    "// Override parameter defaults with values based on command line arguments and NGHOSTS.\n",
    "// We'll use this grid. It has one point and one ghost zone.\n",
    "params.NGHOSTS = 3;\n",
    "params.Nxx0 = atoi(argv[1]);\n",
    "params.Nxx1 = atoi(argv[2]);\n",
    "params.Nxx2 = atoi(argv[3]);\n",
    "params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*params.NGHOSTS;\n",
    "params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*params.NGHOSTS;\n",
    "params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*params.NGHOSTS;\n",
    "// Step 0d: Set up space and time coordinates\n",
    "// Step 0d.i: Declare \\Delta x^i=dxx{0,1,2} and invdxx{0,1,2}, as well as xxmin[3] and xxmax[3]:\n",
    "const REAL xxmin[3] = {-0.1,-0.1,-0.1};\n",
    "const REAL xxmax[3] = { 0.1, 0.1, 0.1};\n",
    "\n",
    "params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx0);\n",
    "params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx1);\n",
    "params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx2);\n",
    "//printf(\"dxx0,dxx1,dxx2 = %.5e,%.5e,%.5e\\\\n\",params.dxx0,params.dxx1,params.dxx2);\n",
    "params.invdx0 = 1.0 / params.dxx0;\n",
    "params.invdx1 = 1.0 / params.dxx1;\n",
    "params.invdx2 = 1.0 / params.dxx2;\n",
    "\\n\"\"\")\n",
    "\n",
    "# Generates declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n",
    "par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='module'></a>\n",
    "\n",
    "## Step 1.f: Generate `GiRaFFE_NRPy` Files \\[Back to [top](#toc)\\]\n",
    "$$\\label{module}$$\n",
    "\n",
    "Here, we generate the functions we want to test by calling the function found [here](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_PPM.py) and documented in [this tutorial](Tutorial-Start_to_Finish-GiRaFFE_NRPy-PPM.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.433683Z",
     "iopub.status.busy": "2020-12-17T02:17:01.433277Z",
     "iopub.status.idle": "2020-12-17T02:17:01.435087Z",
     "shell.execute_reply": "2020-12-17T02:17:01.435423Z"
    }
   },
   "outputs": [],
   "source": [
    "subdir = \"PPM\"\n",
    "cmd.mkdir(os.path.join(Ccodesdir,subdir))\n",
    "import GiRaFFE_NRPy.GiRaFFE_NRPy_PPM as PPM\n",
    "PPM.GiRaFFE_NRPy_PPM(os.path.join(Ccodesdir,subdir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='electric_field'></a>\n",
    "\n",
    "## Step 1.g: Calculate the $E_i$ field contribution to $\\partial_t A_i$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{electric_field}$$\n",
    "\n",
    "Here, we generate the functions necessary to calculate the electric flux on the cell faces, which is the algorithm we are specifically trying to test. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.440450Z",
     "iopub.status.busy": "2020-12-17T02:17:01.440016Z",
     "iopub.status.idle": "2020-12-17T02:17:01.566425Z",
     "shell.execute_reply": "2020-12-17T02:17:01.565935Z"
    }
   },
   "outputs": [],
   "source": [
    "import GiRaFFE_NRPy.GiRaFFE_NRPy_Afield_flux_handwritten as Af\n",
    "\n",
    "# We will pass values of the gridfunction on the cell faces into the function. This requires us\n",
    "# to declare them as C parameters in NRPy+. We will denote this with the _face infix/suffix.\n",
    "alpha_face = gri.register_gridfunctions(\"AUXEVOL\",\"alpha_face\")\n",
    "gamma_faceDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gamma_faceDD\",\"sym01\")\n",
    "beta_faceU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"beta_faceU\")\n",
    "\n",
    "# We'll need some more gridfunctions, now, to represent the reconstructions of BU and ValenciavU\n",
    "# on the right and left faces\n",
    "Valenciav_rU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"Valenciav_rU\",DIM=3)\n",
    "B_rU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"B_rU\",DIM=3)\n",
    "Valenciav_lU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"Valenciav_lU\",DIM=3)\n",
    "B_lU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"B_lU\",DIM=3)\n",
    "\n",
    "AD = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"AD\",DIM=3)\n",
    "\n",
    "subdir = \"RHSs\"\n",
    "cmd.mkdir(os.path.join(Ccodesdir,subdir))\n",
    "Af.GiRaFFE_NRPy_Afield_flux(os.path.join(Ccodesdir,subdir))\n",
    "\n",
    "subdir = \"FCVAL\"\n",
    "cmd.mkdir(os.path.join(Ccodesdir,subdir))\n",
    "# We also need to interpolate the metric values to cell faces, which is easily done with this module:\n",
    "import GiRaFFE_NRPy.GiRaFFE_NRPy_Metric_Face_Values as FCVAL\n",
    "FCVAL.GiRaFFE_NRPy_FCVAL(os.path.join(Ccodesdir,subdir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='exact_flux'></a>\n",
    "\n",
    "## Step 1.h: Calculate the *exact* flux of $\\epsilon_{ijk} v^j B^k$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{exact_flux}$$\n",
    "\n",
    "Here, we generate a function to analytically calculate the electric flux on the cell faces for comparison. We'll need to import the Levi-Civita tensor for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.572214Z",
     "iopub.status.busy": "2020-12-17T02:17:01.571668Z",
     "iopub.status.idle": "2020-12-17T02:17:01.728683Z",
     "shell.execute_reply": "2020-12-17T02:17:01.728147Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Output C function calculate_E_exactD() to file Start-to-Finish-UnitTests/Afield_flux_UnitTest/calculate_E_exactD.h\n"
     ]
    }
   ],
   "source": [
    "import GRHD.equations as gh\n",
    "gh.compute_sqrtgammaDET(gammaDD)\n",
    "LeviCivitaDDD = ixp.LeviCivitaTensorDDD_dim3_rank3(gh.sqrtgammaDET)\n",
    "\n",
    "driftvU = ixp.zerorank1()\n",
    "for i in range(3):\n",
    "    # TODO: don't use _face metric gridfunctions once we're in curved space!\n",
    "    driftvU[i] = alpha*ValenciavU[i]-betaU[i]\n",
    "\n",
    "A_rhsD = ixp.zerorank1()\n",
    "for i in range(3):\n",
    "    for j in range(3):\n",
    "        for k in range(3):\n",
    "            A_rhsD[i] += LeviCivitaDDD[i][j][k]*driftvU[j]*BU[k]\n",
    "\n",
    "A_rhsD_to_print = [\n",
    "                   lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"AD0\"),rhs=A_rhsD[0]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"AD1\"),rhs=A_rhsD[1]),\n",
    "                   lhrh(lhs=gri.gfaccess(\"rhs_gfs\",\"AD2\"),rhs=A_rhsD[2]),\n",
    "                  ]\n",
    "\n",
    "desc = \"Calculate analytic electric field, part of the right-hand side of AD\"\n",
    "name = \"calculate_E_exactD\"\n",
    "outCfunction(\n",
    "    outfile  = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n",
    "    params   =\"const paramstruct *params,REAL *xx[3],const REAL *auxevol_gfs,REAL *rhs_gfs\",\n",
    "    body     = fin.FD_outputC(\"returnstring\",A_rhsD_to_print,params=\"outCverbose=False\"),\n",
    "    loopopts =\"AllPoints,Read_xxs\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='mainc'></a>\n",
    "\n",
    "# Step 2: `Afield_flux_unit_test.c`: The Main C Code \\[Back to [top](#toc)\\]\n",
    "$$\\label{mainc}$$\n",
    "\n",
    "Now that we have our vector potential and analytic magnetic field to compare against, we will start writing our unit test. We'll also import common C functionality, define `REAL`, the number of ghost zones, and the faces, and set the standard macros for NRPy+ style memory access. After the preliminaries and setting data on our grid, we compare the output between `calculate_E_exactD()` and `calculate_E_field_flat_all_in_one()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.731745Z",
     "iopub.status.busy": "2020-12-17T02:17:01.731232Z",
     "iopub.status.idle": "2020-12-17T02:17:01.733685Z",
     "shell.execute_reply": "2020-12-17T02:17:01.733238Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing Start-to-Finish-UnitTests/Afield_flux_UnitTest//Afield_flux_unit_test.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile $Ccodesdir/Afield_flux_unit_test.c\n",
    "\n",
    "// These are common packages that we are likely to need.\n",
    "#include \"stdio.h\"\n",
    "#include \"stdlib.h\"\n",
    "#include \"math.h\"\n",
    "#include \"string.h\" // Needed for strncmp, etc.\n",
    "#include \"stdint.h\" // Needed for Windows GCC 6.x compatibility\n",
    "#include <time.h>   // Needed to set a random seed.\n",
    "\n",
    "#define REAL double\n",
    "#include \"declare_Cparameters_struct.h\"\n",
    "\n",
    "REAL a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t;\n",
    "\n",
    "// Standard NRPy+ memory access:\n",
    "#define IDX3S(i,j,k) ( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * (k) ) )\n",
    "#define IDX4S(g,i,j,k) \\\n",
    "( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )\n",
    "#define IDX4ptS(g,idx) ( (idx) + (Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2) * (g) )\n",
    "\n",
    "#define LOOP_REGION(i0min,i0max, i1min,i1max, i2min,i2max) \\\n",
    "  for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll now define the gridfunction names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.736891Z",
     "iopub.status.busy": "2020-12-17T02:17:01.736420Z",
     "iopub.status.idle": "2020-12-17T02:17:01.738208Z",
     "shell.execute_reply": "2020-12-17T02:17:01.738550Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to Start-to-Finish-UnitTests/Afield_flux_UnitTest//Afield_flux_unit_test.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $Ccodesdir/Afield_flux_unit_test.c\n",
    "\n",
    "// Let's also #define the NRPy+ gridfunctions\n",
    "#define VALENCIAVU0GF 0\n",
    "#define VALENCIAVU1GF 1\n",
    "#define VALENCIAVU2GF 2\n",
    "#define BU0GF 3\n",
    "#define BU1GF 4\n",
    "#define BU2GF 5\n",
    "#define VALENCIAV_RU0GF 6\n",
    "#define VALENCIAV_RU1GF 7\n",
    "#define VALENCIAV_RU2GF 8\n",
    "#define B_RU0GF 9\n",
    "#define B_RU1GF 10\n",
    "#define B_RU2GF 11\n",
    "#define VALENCIAV_LU0GF 12\n",
    "#define VALENCIAV_LU1GF 13\n",
    "#define VALENCIAV_LU2GF 14\n",
    "#define B_LU0GF 15\n",
    "#define B_LU1GF 16\n",
    "#define B_LU2GF 17\n",
    "#define GAMMA_FACEDD00GF 18\n",
    "#define GAMMA_FACEDD01GF 19\n",
    "#define GAMMA_FACEDD02GF 20\n",
    "#define GAMMA_FACEDD11GF 21\n",
    "#define GAMMA_FACEDD12GF 22\n",
    "#define GAMMA_FACEDD22GF 23\n",
    "#define BETA_FACEU0GF 24\n",
    "#define BETA_FACEU1GF 25\n",
    "#define BETA_FACEU2GF 26\n",
    "#define ALPHA_FACEGF 27\n",
    "#define GAMMADD00GF 28\n",
    "#define GAMMADD01GF 29\n",
    "#define GAMMADD02GF 30\n",
    "#define GAMMADD11GF 31\n",
    "#define GAMMADD12GF 32\n",
    "#define GAMMADD22GF 33\n",
    "#define BETAU0GF 34\n",
    "#define BETAU1GF 35\n",
    "#define BETAU2GF 36\n",
    "#define ALPHAGF 37\n",
    "#define GAMMA_FACEUU00GF 38\n",
    "#define GAMMA_FACEUU01GF 39\n",
    "#define GAMMA_FACEUU02GF 40\n",
    "#define GAMMA_FACEUU11GF 41\n",
    "#define GAMMA_FACEUU12GF 42\n",
    "#define GAMMA_FACEUU22GF 43\n",
    "#define GAMMAUU00GF 44\n",
    "#define GAMMAUU01GF 45\n",
    "#define GAMMAUU02GF 46\n",
    "#define GAMMAUU11GF 47\n",
    "#define GAMMAUU12GF 48\n",
    "#define GAMMAUU22GF 49\n",
    "#define NUM_AUXEVOL_GFS 50\n",
    "\n",
    "#define AD0GF 0\n",
    "#define AD1GF 1\n",
    "#define AD2GF 2\n",
    "#define NUM_EVOL_GFS 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we'll handle the different A2B codes. There are several things to do here. First, we'll add `#include`s to the C code so that we have access to the functions we want to test. We must also create a directory and copy the files to that directory. We will choose to do this in the subfolder `A2B` relative to this tutorial.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.741301Z",
     "iopub.status.busy": "2020-12-17T02:17:01.740895Z",
     "iopub.status.idle": "2020-12-17T02:17:01.742828Z",
     "shell.execute_reply": "2020-12-17T02:17:01.743162Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to Start-to-Finish-UnitTests/Afield_flux_UnitTest//Afield_flux_unit_test.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $Ccodesdir/Afield_flux_unit_test.c\n",
    "\n",
    "// Some specific definitions needed for this file\n",
    "typedef struct __gf_and_gz_struct__ {\n",
    "  REAL *gf;\n",
    "  int gz_lo[4],gz_hi[4];\n",
    "} gf_and_gz_struct;\n",
    "\n",
    "const int VX=0,VY=1,VZ=2,BX=3,BY=4,BZ=5;\n",
    "const int NUM_RECONSTRUCT_GFS = 6;\n",
    "\n",
    "#include \"PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c\"\n",
    "#include \"PPM/loop_defines_reconstruction_NRPy.h\"\n",
    "#include \"FCVAL/interpolate_metric_gfs_to_cell_faces.h\"\n",
    "\n",
    "#include \"calculate_metric_gfs.h\"\n",
    "#include \"calculate_BU.h\"\n",
    "#include \"calculate_ValenciavU.h\"\n",
    "#include \"calculate_E_exactD.h\"\n",
    "\n",
    "// These are the functions we want to test.\n",
    "#include \"RHSs/calculate_E_field_flat_all_in_one.h\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we'll write the main method. First, we'll set up the grid. In this test, we cannot use only one point. To properly perform reconstruction and have a point at the center of the grid that's ideal for convergence testing, we'll need at least a $9 \\times 9 \\times 9$ grid; for the higher resolution, we'll double the interior points in one direction, yielding a $11 \\times 11 \\times 11$ grid. Then, we'll write the magnetic fields and drift velocity. After that, we'll calculate the electric field contribution to the $A_i$ RHS in two ways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.746761Z",
     "iopub.status.busy": "2020-12-17T02:17:01.744855Z",
     "iopub.status.idle": "2020-12-17T02:17:01.749221Z",
     "shell.execute_reply": "2020-12-17T02:17:01.748812Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Appending to Start-to-Finish-UnitTests/Afield_flux_UnitTest//Afield_flux_unit_test.c\n"
     ]
    }
   ],
   "source": [
    "%%writefile -a $Ccodesdir/Afield_flux_unit_test.c\n",
    "\n",
    "int main(int argc, const char *argv[]) {\n",
    "    paramstruct params;\n",
    "#include \"set_Cparameters_default.h\"\n",
    "\n",
    "    // Step 0c: Set free parameters, overwriting Cparameters defaults\n",
    "    //          by hand or with command-line input, as desired.\n",
    "#include \"free_parameters.h\"\n",
    "#include \"set_Cparameters-nopointer.h\"\n",
    "\n",
    "    // We'll define our grid slightly different from how we normally would. We let our outermost\n",
    "    // ghostzones coincide with xxmin and xxmax instead of the interior of the grid. This means\n",
    "    // that the ghostzone points will have identical positions so we can do convergence tests of them.    // Step 0d.ii: Set up uniform coordinate grids\n",
    "    REAL *xx[3];\n",
    "    xx[0] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS0);\n",
    "    xx[1] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS1);\n",
    "    xx[2] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS2);\n",
    "    for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] = xxmin[0] + (j-NGHOSTS)*dxx0;\n",
    "    for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + (j-NGHOSTS)*dxx1;\n",
    "    for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + (j-NGHOSTS)*dxx2;\n",
    "\n",
    "    //for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) printf(\"x[%d] = %.5e\\n\",j,xx[0][j]);\n",
    "\n",
    "    // This is the array to which we'll write the NRPy+ variables.\n",
    "    REAL *auxevol_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);\n",
    "    REAL *rhs_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);\n",
    "    REAL *rhs_exact_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);\n",
    "\n",
    "    a = (double)(rand()%20-10);\n",
    "    b = (double)(rand()%20-10);\n",
    "    c = (double)(rand()%20-10);\n",
    "    d = (double)(rand()%20-10);\n",
    "    e = (double)(rand()%20-10);\n",
    "    f = (double)(rand()%20-10);\n",
    "    g = (double)(rand()%20-10);\n",
    "    h = (double)(rand()%20-10);\n",
    "    l = (double)(rand()%20-10);\n",
    "    m = (double)(rand()%20-10);\n",
    "    n = (double)(rand()%20-10);\n",
    "    o = (double)(rand()%20-10);\n",
    "    p = (double)(rand()%20-10);\n",
    "    q = (double)(rand()%20-10);\n",
    "    r = (double)(rand()%20-10);\n",
    "    s = (double)(rand()%20-10);\n",
    "    t = (double)(rand()%20-10);\n",
    "    //printf(\"a,b,c,d,e,f = %f,%f,%f,%f,%f,%f\\n\",a,b,c,d,e,f);\n",
    "\n",
    "    // First, calculate the test metric on our grid:\n",
    "    calculate_metric_gfs(&params,xx,auxevol_gfs);\n",
    "\n",
    "    // Calculate the initial data and the exact solution.\n",
    "    calculate_BU(&params,xx,auxevol_gfs);\n",
    "    calculate_ValenciavU(&params,xx,auxevol_gfs);\n",
    "    //for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] -= 0.5*dxx0;\n",
    "    calculate_E_exactD(&params,xx,auxevol_gfs,rhs_exact_gfs);\n",
    "\n",
    "\n",
    "    // Set up the pointers and constants for the reconstruction procedure.\n",
    "    gf_and_gz_struct in_prims[NUM_RECONSTRUCT_GFS], out_prims_r[NUM_RECONSTRUCT_GFS], out_prims_l[NUM_RECONSTRUCT_GFS];\n",
    "\n",
    "    const int Nxxp2NG012 = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;\n",
    "    REAL temporary[Nxxp2NG012];\n",
    "\n",
    "    int ww=0;\n",
    "    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*VALENCIAVU0GF;\n",
    "      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU0GF;\n",
    "      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU0GF;\n",
    "    ww++;\n",
    "    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*VALENCIAVU1GF;\n",
    "      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU1GF;\n",
    "      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU1GF;\n",
    "    ww++;\n",
    "    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*VALENCIAVU2GF;\n",
    "      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU2GF;\n",
    "      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU2GF;\n",
    "    ww++;\n",
    "    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*BU0GF;\n",
    "      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU0GF;\n",
    "      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU0GF;\n",
    "    ww++;\n",
    "    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*BU1GF;\n",
    "      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU1GF;\n",
    "      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU1GF;\n",
    "    ww++;\n",
    "    in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*BU2GF;\n",
    "      out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU2GF;\n",
    "      out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU2GF;\n",
    "    ww++;\n",
    "\n",
    "    // Prims are defined AT ALL GRIDPOINTS, so we set the # of ghostzones to zero:\n",
    "    for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { in_prims[i].gz_lo[j]=0; in_prims[i].gz_hi[j]=0; }\n",
    "    // Left/right variables are not yet defined, yet we set the # of gz's to zero by default:\n",
    "    for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { out_prims_r[i].gz_lo[j]=0; out_prims_r[i].gz_hi[j]=0; }\n",
    "    for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { out_prims_l[i].gz_lo[j]=0; out_prims_l[i].gz_hi[j]=0; }\n",
    "\n",
    "   int which_prims_to_reconstruct[6] = {VX, VY, VZ, BX, BY, BZ};\n",
    "   int num_prims_to_reconstruct = 6;\n",
    "\n",
    "    // In each direction, perform the PPM reconstruction procedure.\n",
    "    // Then, add the fluxes to the RHS as appropriate.\n",
    "    for(int flux_dirn=0;flux_dirn<3;flux_dirn++) {\n",
    "        // This function is housed in the file: \"reconstruct_set_of_prims_PPM_GRFFE_NRPy.c\"\n",
    "        interpolate_metric_gfs_to_cell_faces(&params,auxevol_gfs,flux_dirn+1);\n",
    "        reconstruct_set_of_prims_PPM_GRFFE_NRPy(&params, auxevol_gfs, flux_dirn+1, num_prims_to_reconstruct,\n",
    "                                                which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);\n",
    "        for(int count=0;count<=1;count++) {\n",
    "            // This function is written to be general, using notation that matches the forward permutation added to AD2,\n",
    "            // i.e., [F_HLL^x(B^y)]_z corresponding to flux_dirn=0, count=1.\n",
    "            // The SIGN parameter is necessary because\n",
    "            // -E_z(x_i,y_j,z_k) = 0.25 ( [F_HLL^x(B^y)]_z(i+1/2,j,k)+[F_HLL^x(B^y)]_z(i-1/2,j,k)\n",
    "            //                           -[F_HLL^y(B^x)]_z(i,j+1/2,k)-[F_HLL^y(B^x)]_z(i,j-1/2,k) )\n",
    "            // Note the negative signs on the reversed permutation terms!\n",
    "\n",
    "            // By cyclically permuting with flux_dirn, we\n",
    "            // get contributions to the other components, and by incrementing count, we get the backward permutations:\n",
    "            // Let's suppose flux_dirn = 0. Then we will need to update Ay (count=0) and Az (count=1):\n",
    "            //     flux_dirn=count=0 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (0+1+0)%3=AD1GF <- Updating Ay!\n",
    "            //        (flux_dirn)%3 = (0)%3 = 0               Vx\n",
    "            //        (flux_dirn-count+2)%3 = (0-0+2)%3 = 2   Vz .  Inputs Vx, Vz -> SIGN = -1 ; 2.0*((REAL)count)-1.0=-1 check!\n",
    "            //     flux_dirn=0,count=1 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (0+1+1)%3=AD2GF <- Updating Az!\n",
    "            //        (flux_dirn)%3 = (0)%3 = 0               Vx\n",
    "            //        (flux_dirn-count+2)%3 = (0-1+2)%3 = 1   Vy .  Inputs Vx, Vy -> SIGN = +1 ; 2.0*((REAL)count)-1.0=2-1=+1 check!\n",
    "            // Let's suppose flux_dirn = 1. Then we will need to update Az (count=0) and Ax (count=1):\n",
    "            //     flux_dirn=1,count=0 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (1+1+0)%3=AD2GF <- Updating Az!\n",
    "            //        (flux_dirn)%3 = (1)%3 = 1               Vy\n",
    "            //        (flux_dirn-count+2)%3 = (1-0+2)%3 = 0   Vx .  Inputs Vy, Vx -> SIGN = -1 ; 2.0*((REAL)count)-1.0=-1 check!\n",
    "            //     flux_dirn=count=1 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (1+1+1)%3=AD0GF <- Updating Ax!\n",
    "            //        (flux_dirn)%3 = (1)%3 = 1               Vy\n",
    "            //        (flux_dirn-count+2)%3 = (1-1+2)%3 = 2   Vz .  Inputs Vy, Vz -> SIGN = +1 ; 2.0*((REAL)count)-1.0=2-1=+1 check!\n",
    "            // Let's suppose flux_dirn = 2. Then we will need to update Ax (count=0) and Ay (count=1):\n",
    "            //     flux_dirn=2,count=0 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (2+1+0)%3=AD0GF <- Updating Ax!\n",
    "            //        (flux_dirn)%3 = (2)%3 = 2               Vz\n",
    "            //        (flux_dirn-count+2)%3 = (2-0+2)%3 = 1   Vy .  Inputs Vz, Vy -> SIGN = -1 ; 2.0*((REAL)count)-1.0=-1 check!\n",
    "            //     flux_dirn=2,count=1 -> AD0GF+(flux_dirn+1+count)%3 = AD0GF + (2+1+1)%3=AD1GF <- Updating Ay!\n",
    "            //        (flux_dirn)%3 = (2)%3 = 2               Vz\n",
    "            //        (flux_dirn-count+2)%3 = (2-1+2)%3 = 0   Vx .  Inputs Vz, Vx -> SIGN = +1 ; 2.0*((REAL)count)-1.0=2-1=+1 check!\n",
    "            calculate_E_field_flat_all_in_one(&params,\n",
    "              &auxevol_gfs[IDX4ptS(VALENCIAV_RU0GF+(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(VALENCIAV_RU0GF+(flux_dirn-count+2)%3, 0)],\n",
    "              &auxevol_gfs[IDX4ptS(VALENCIAV_LU0GF+(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(VALENCIAV_LU0GF+(flux_dirn-count+2)%3, 0)],\n",
    "              &auxevol_gfs[IDX4ptS(B_RU0GF        +(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(B_RU0GF        +(flux_dirn-count+2)%3, 0)],\n",
    "              &auxevol_gfs[IDX4ptS(B_LU0GF        +(flux_dirn)%3, 0)],&auxevol_gfs[IDX4ptS(B_LU0GF        +(flux_dirn-count+2)%3, 0)],\n",
    "              &auxevol_gfs[IDX4ptS(B_RU0GF        +(flux_dirn-count+2)%3, 0)],\n",
    "              &auxevol_gfs[IDX4ptS(B_LU0GF        +(flux_dirn-count+2)%3, 0)],\n",
    "              &auxevol_gfs[IDX4ptS(GAMMA_FACEDD00GF,0)],&auxevol_gfs[IDX4ptS(GAMMA_FACEDD01GF,0)],&auxevol_gfs[IDX4ptS(GAMMA_FACEDD02GF,0)],\n",
    "              &auxevol_gfs[IDX4ptS(GAMMA_FACEDD11GF,0)],&auxevol_gfs[IDX4ptS(GAMMA_FACEDD12GF,0)],&auxevol_gfs[IDX4ptS(GAMMA_FACEDD22GF,0)],\n",
    "              &auxevol_gfs[IDX4ptS(BETA_FACEU0GF+flux_dirn,0)],&auxevol_gfs[IDX4ptS(BETA_FACEU0GF+(flux_dirn-count+2)%3,0)],&auxevol_gfs[IDX4ptS(ALPHA_FACEGF,0)],\n",
    "              &rhs_gfs[IDX4ptS(AD0GF+(flux_dirn+1+count)%3,0)], 2.0*((REAL)count)-1.0, flux_dirn);\n",
    "        }\n",
    "    }\n",
    "\n",
    "    char filename[100];\n",
    "    sprintf(filename,\"out%d-numer.txt\",Nxx0);\n",
    "    FILE *out2D = fopen(filename, \"w\");\n",
    "    // We print the difference between approximate and exact numbers.\n",
    "    int i0 = Nxx_plus_2NGHOSTS0/2;\n",
    "    int i1 = Nxx_plus_2NGHOSTS1/2;\n",
    "    int i2 = Nxx_plus_2NGHOSTS2/2;\n",
    "\n",
    "    printf(\"Numerical: %.15e\\n\",rhs_gfs[IDX4S(AD0GF,i0,i1,i2)]);\n",
    "    printf(\"Analytic:  %.15e\\n\",rhs_exact_gfs[IDX4S(AD0GF,i0,i1,i2)]);\n",
    "    printf(\"Numerical: %.15e\\n\",rhs_gfs[IDX4S(AD1GF,i0,i1,i2)]);\n",
    "    printf(\"Analytic:  %.15e\\n\",rhs_exact_gfs[IDX4S(AD1GF,i0,i1,i2)]);\n",
    "    printf(\"Numerical: %.15e\\n\",rhs_gfs[IDX4S(AD2GF,i0,i1,i2)]);\n",
    "    printf(\"Analytic:  %.15e\\n\\n\",rhs_exact_gfs[IDX4S(AD2GF,i0,i1,i2)]);\n",
    "\n",
    "    //printf(\"i0,i1,i2 = %d,%d,%d\\n\",i0,i1,i2);\n",
    "    fprintf(out2D,\"%.16e\\t%.16e\\t%.16e\\t %e %e %e\\n\",\n",
    "            rhs_exact_gfs[IDX4S(AD0GF,i0,i1,i2)]-rhs_gfs[IDX4S(AD0GF,i0,i1,i2)],\n",
    "            rhs_exact_gfs[IDX4S(AD1GF,i0,i1,i2)]-rhs_gfs[IDX4S(AD1GF,i0,i1,i2)],\n",
    "            rhs_exact_gfs[IDX4S(AD2GF,i0,i1,i2)]-rhs_gfs[IDX4S(AD2GF,i0,i1,i2)],\n",
    "            xx[0][i0],xx[1][i1],xx[2][i2]\n",
    "            );\n",
    "    fclose(out2D);\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='compile_run'></a>\n",
    "\n",
    "## Step 2.a: Compile and run the code\n",
    "\n",
    "$$\\label{compile_run}$$\n",
    "\n",
    "Now that we have our file, we can compile it and run the executable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:01.754217Z",
     "iopub.status.busy": "2020-12-17T02:17:01.753819Z",
     "iopub.status.idle": "2020-12-17T02:17:02.611865Z",
     "shell.execute_reply": "2020-12-17T02:17:02.611413Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compiling executable...\n",
      "(EXEC): Executing `gcc -Ofast -fopenmp -march=native -funroll-loops Start-to-Finish-UnitTests/Afield_flux_UnitTest/Afield_flux_unit_test.c -o Start-to-Finish-UnitTests/Afield_flux_UnitTest/output/Afield_flux_unit_test -lm`...\n",
      "(BENCH): Finished executing in 0.6064903736114502 seconds.\n",
      "Finished compilation.\n",
      "(EXEC): Executing `taskset -c 0,1,2,3,4,5 ./Afield_flux_unit_test 2 2 2`...\n",
      "(BENCH): Finished executing in 0.20341801643371582 seconds.\n",
      "(EXEC): Executing `taskset -c 0,1,2,3,4,5 ./Afield_flux_unit_test 4 4 4`...\n",
      "(BENCH): Finished executing in 0.20840764045715332 seconds.\n"
     ]
    }
   ],
   "source": [
    "cmd.C_compile(os.path.join(Ccodesdir,\"Afield_flux_unit_test.c\"), os.path.join(outdir,\"Afield_flux_unit_test\"))\n",
    "\n",
    "# Change to output directory\n",
    "os.chdir(outdir)\n",
    "\n",
    "cmd.Execute(os.path.join(\"Afield_flux_unit_test\"),\"2 2 2\",\"out2.txt\")\n",
    "# To do a convergence test, we'll also need a second grid with twice the resolution.\n",
    "cmd.Execute(os.path.join(\"Afield_flux_unit_test\"),\"4 4 4\",\"out4.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='convergence'></a>\n",
    "\n",
    "# Step 3: Code validation: Verify that relative error in numerical solution converges to zero at the expected order \\[Back to [top](#toc)\\]\n",
    "$$\\label{convergence}$$\n",
    "\n",
    "For testing purposes, we have only checked these algorithms on a small grid. By construction, we have only guaranteed ourselves output from the functions we are testing at a point, so we will simply print the convergence order at that point after processing our outputs below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:02.617734Z",
     "iopub.status.busy": "2020-12-17T02:17:02.617125Z",
     "iopub.status.idle": "2020-12-17T02:17:02.786296Z",
     "shell.execute_reply": "2020-12-17T02:17:02.786723Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The following quantities converge at the listed order (should be ~1 for Shock data, ~2 otherwise):\n",
      "A_rhsD0: 2.0234555321265675\n",
      "A_rhsD1: 1.0081097780531285\n",
      "A_rhsD2: 1.1291244848636108\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "Data1 = np.loadtxt(\"out2-numer.txt\")\n",
    "Data2 = np.loadtxt(\"out4-numer.txt\")\n",
    "\"\"\n",
    "def conv_order(lowres,highres):\n",
    "    order = np.log2(np.abs(lowres/highres))\n",
    "    if Use_Shock_Data:\n",
    "        if order < 0.7 or order > 2.3:\n",
    "            sys.exit(1)\n",
    "    else:\n",
    "        if order < 1.7 or order > 2.3:\n",
    "            sys.exit(1)\n",
    "    return order\n",
    "\n",
    "\n",
    "print(\"The following quantities converge at the listed order (should be ~1 for Shock data, ~2 otherwise):\")\n",
    "print(\"A_rhsD0: \"+str(conv_order(Data1[0],Data2[0])))\n",
    "print(\"A_rhsD1: \"+str(conv_order(Data1[1],Data2[1])))\n",
    "print(\"A_rhsD2: \"+str(conv_order(Data1[2],Data2[2])))"
   ]
  },
  {
   "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",
    "\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-Start_to_Finish_UnitTest-GiRaFFE_NRPy-A2B.pdf](Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-A2B.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": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2020-12-17T02:17:02.789397Z",
     "iopub.status.busy": "2020-12-17T02:17:02.788891Z",
     "iopub.status.idle": "2020-12-17T02:17:05.646368Z",
     "shell.execute_reply": "2020-12-17T02:17:05.646845Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-Afield_flux.tex, and\n",
      "    compiled LaTeX file to PDF file Tutorial-Start_to_Finish_UnitTest-\n",
      "    GiRaFFE_NRPy-Afield_flux.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "\n",
    "# Change to NRPy directory\n",
    "os.chdir(\"../../../\")\n",
    "\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-Afield_flux\",location_of_template_file=os.path.join(\"..\"))"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}