{
 "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",
    "# 1D Three Wave `GiRaFFEfood` Initial Data for `GiRaFFE`\n",
    "\n",
    "## This module provides another initial data option for `GiRaFFE`, drawn from [this paper](https://arxiv.org/abs/1310.3274) .\n",
    "\n",
    "**Notebook Status:** <font color='green'><b> Validated </b></font>\n",
    "\n",
    "**Validation Notes:** This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented [below](#code_validation). The initial data has validated against the original `GiRaFFE`, as documented [here](Tutorial-Start_to_Finish_UnitTest-GiRaFFEfood_NRPy.ipynb).\n",
    "\n",
    "### NRPy+ Source Code for this module: [GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Three_Waves.py](../../edit/in_progress/GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Three_Waves.py)\n",
    "\n",
    "## Introduction:\n",
    "\n",
    "### Three waves:\n",
    "\n",
    "This is a flat-spacetime test representing three Alfv&eacute;n waves (one stationary, one left-going, and one right-going) with initial data \n",
    "\\begin{align}\n",
    "A_x &= 0 \\\\ \n",
    "A_y &= 3.5x H(-x) + 3.0x H(x) \\\\\n",
    "A_z &= y - 1.5x H(-x) - 3.0x H(x),\n",
    "\\end{align}\n",
    "where $H(x)$ is the Heaviside function, which generates the magnetic field\n",
    "$$\\mathbf{B}(0,x) = \\mathbf{B_a}(0,x) + \\mathbf{B_+}(0,x) + \\mathbf{B_-}(0,x)$$\n",
    "and uses the electric field\n",
    "$$\\mathbf{E}(0,x) = \\mathbf{E_a}(0,x) + \\mathbf{E_+}(0,x) + \\mathbf{E_-}(0,x),$$\n",
    "where subscripted $\\mathbf{a}$ corresponds to the stationary wave, subscripted $\\mathbf{+}$ corresponds to the right-going wave, and subscripted $\\mathbf{-}$ corresponds to the left-going wave, and where \n",
    "\\begin{align}\n",
    "\\mathbf{B_a}(0,x) &= \\left \\{ \\begin{array}{lll} (1.0,1.0,2.0) & \\mbox{if} & x<0 \\\\\n",
    "\t\t\t\t\t(1.0,1.5,2.0) & \\mbox{if} & x>0 \\end{array} \n",
    "\\right. , \\\\ \n",
    "\\mathbf{E_a}(0,x) &= \\left \\{ \\begin{array}{lll} (-1.0,1.0,0.0) & \\mbox{if} & x<0 \\\\\n",
    "\t\t\t\t\t(-1.5,1.0,0.0) & \\mbox{if} & x>0 \\end{array} \n",
    "\\right.  , \\\\\n",
    "\\mathbf{B_+}(0,x) &= \\left \\{ \\begin{array}{lll} (0.0,0.0,0.0) & \\mbox{if} & x<0 \\\\\n",
    "                                        (0.0,1.5,1.0) & \\mbox{if} & x>0 \\end{array} \n",
    "\\right. , \\\\\n",
    "\\mathbf{E_+}(0,x) &= \\left \\{ \\begin{array}{lll} (0.0,0.0,0.0) & \\mbox{if} & x<0 \\\\\n",
    "                                        (0.0,1.0,-1.5) & \\mbox{if} & x>0 \\end{array} \n",
    "\\right. , \\\\\n",
    "\\mathbf{B_-}(0,x) &= \\left \\{ \\begin{array}{lll} (0.0,0.5,1.5) & \\mbox{if} & x<0 \\\\\n",
    "                                        (0.0,0.0,0.0) & \\mbox{if} & x>0 \\end{array}\n",
    "\\right.  , \\\\\n",
    "\\mathbf{E_-}(0,x) &= \\left \\{ \\begin{array}{lll} (0.0,-1.5,0.5) & \\mbox{if} & x<0 \\\\\n",
    "                                        (0.0,0.0,0.0) & \\mbox{if} & x>0 \\end{array}\n",
    "\\right.  . \\\\\n",
    "\\end{align}\n",
    "\n",
    "For the eventual purpose of testing convergence, any quantity $Q$ evolves as $Q(t,x) = Q_a(0,x) + Q_+(0,x-t) + Q_-(0,x+t)$.\n",
    "\n",
    "See the [Tutorial-GiRaFFEfood_NRPy](Tutorial-GiRaFFEfood_NRPy.ipynb) tutorial notebook for more general detail on how this is used.\n"
   ]
  },
  {
   "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](#initializenrpy): Import core NRPy+ modules and set NRPy+ parameters\n",
    "1. [Step 2](#set_a_i): Set the vector $A_i$\n",
    "1. [Step 3](#set_vi): Calculate $v^i$ from $B^i$ and $E_i$\n",
    "1. [Step 4](#code_validation): Code Validation against `GiRaFFEfood_NRPy.GiRaFFEfood_NRPy` NRPy+ Module\n",
    "1. [Step 5](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='initializenrpy'></a>\n",
    "\n",
    "# Step 1: Import core NRPy+ modules and set NRPy+ parameters \\[Back to [top](#toc)\\]\n",
    "$$\\label{initializenrpy}$$\n",
    "\n",
    "Here, we will import the NRPy+ core modules and set the reference metric to Cartesian, set commonly used NRPy+ parameters, and set C parameters that will be set from outside the code eventually generated from these expressions. We will also set up a parameter to determine what initial data is set up, although it won't do much yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Step 0: Add NRPy's directory to the path\n",
    "# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n",
    "import os,sys\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",
    "# Step 0.a: Import the NRPy+ core modules and set the reference metric to Cartesian\n",
    "import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends\n",
    "import NRPy_param_funcs as par   # NRPy+: Parameter interface\n",
    "import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
    "import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy_Common_Functions as gfcf # Some useful functions for GiRaFFE initial data.\n",
    "\n",
    "import reference_metric as rfm   # NRPy+: Reference metric support\n",
    "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n",
    "rfm.reference_metric()\n",
    "\n",
    "# Step 1a: Set commonly used parameters.\n",
    "thismodule = \"GiRaFFEfood_NRPy_1D\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### <a id='set_a_i'></a>\n",
    "\n",
    "# Step 2: Set the vector $A_i$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{set_a_i}$$\n",
    "\n",
    "The vector potential is given as\n",
    "\\begin{align}\n",
    "A_x &= 0 \\\\ \n",
    "A_y &= 3.5x H(-x) + 3.0x H(x) \\\\\n",
    "A_z &= y - 1.5x H(-x) - 3.0x H(x),\n",
    "\\end{align}\n",
    "\n",
    "However, to take full advantage of NRPy+'s automated function generation capabilities, we want to write this without the `if` statements, replacing them with calls to `fabs()`. To do so, we will use the NRPy+ module `Min_Max_and_Piecewise_Expressions`. We thus get\n",
    "$$H(x) = \\frac{\\max(0,x)}{x}.$$\n",
    "This implementation is, of course, undefined for $x=0$; this problem is easily solved by adding a very small number (called `TINYDOUBLE` in our implementation) to the denominator (see [Tutorial-Min_Max_and_Piecewise_Expressions](Tutorial-Min_Max_and_Piecewise_Expressions.ipynb) for details on how this works). This is, conveniently, the exact implementation of the `coord_greater_bound()` function!\n",
    "\n",
    "\\begin{align}\n",
    "A_x &= 0 \\\\ \n",
    "A_y &= 3.5x H(-x) + 3.0x H(x) \\\\\n",
    "A_z &= y - 1.5x H(-x) - 3.0x H(x),\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Min_Max_and_Piecewise_Expressions as noif\n",
    "\n",
    "def Ax_TW(x,y,z, **params):\n",
    "    return sp.sympify(0)\n",
    "\n",
    "def Ay_TW(x,y,z, **params):\n",
    "    return sp.Rational(7,2)*x*noif.coord_greater_bound(-x,0) + sp.sympify(3)*x*noif.coord_greater_bound(x,0)\n",
    "\n",
    "def Az_TW(x,y,z, **params):\n",
    "    return y-sp.Rational(3,2)*x*noif.coord_greater_bound(-x,0) - sp.sympify(3)*x*noif.coord_greater_bound(x,0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='set_vi'></a>\n",
    "\n",
    "# Step 3: Calculate $v^i$ from $B^i$ and $E_i$ \\[Back to [top](#toc)\\]\n",
    "$$\\label{set_vi}$$\n",
    "\n",
    "First, we will set the three individual waves; we change all $<$ to $\\leq$ to avoid unintented behavior at $x=0$:\n",
    "\\begin{align}\n",
    "\\mathbf{B_a}(0,x) &= \\left \\{ \\begin{array}{lll} (1.0,1.0,2.0) & \\mbox{if} & x \\leq 0 \\\\\n",
    "\t\t\t\t\t(1.0,1.5,2.0) & \\mbox{if} & x>0 \\end{array} \n",
    "\\right. , \\\\ \n",
    "\\mathbf{E_a}(0,x) &= \\left \\{ \\begin{array}{lll} (-1.0,1.0,0.0) & \\mbox{if} & x \\leq 0 \\\\\n",
    "\t\t\t\t\t(-1.5,1.0,0.0) & \\mbox{if} & x>0 \\end{array} \n",
    "\\right.  , \\\\\n",
    "\\mathbf{B_+}(0,x) &= \\left \\{ \\begin{array}{lll} (0.0,0.0,0.0) & \\mbox{if} & x \\leq 0 \\\\\n",
    "                                        (0.0,1.5,1.0) & \\mbox{if} & x>0 \\end{array} \n",
    "\\right. , \\\\\n",
    "\\mathbf{E_+}(0,x) &= \\left \\{ \\begin{array}{lll} (0.0,0.0,0.0) & \\mbox{if} & x \\leq 0 \\\\\n",
    "                                        (0.0,1.0,-1.5) & \\mbox{if} & x>0 \\end{array} \n",
    "\\right. , \\\\\n",
    "\\mathbf{B_-}(0,x) &= \\left \\{ \\begin{array}{lll} (0.0,0.5,1.5) & \\mbox{if} & x \\leq 0 \\\\\n",
    "                                        (0.0,0.0,0.0) & \\mbox{if} & x>0 \\end{array}\n",
    "\\right.  , \\\\\n",
    "\\mathbf{E_-}(0,x) &= \\left \\{ \\begin{array}{lll} (0.0,-1.5,0.5) & \\mbox{if} & x \\leq 0 \\\\\n",
    "                                        (0.0,0.0,0.0) & \\mbox{if} & x>0 \\end{array}\n",
    "\\right.  . \\\\\n",
    "\\end{align}\n",
    "\n",
    "Then, we can obtain the total expressions for the magnetic and electric fields by simply adding the three waves together:\n",
    "\\begin{align}\n",
    "\\mathbf{B}(0,x) &= \\mathbf{B_a}(0,x) + \\mathbf{B_+}(0,x) + \\mathbf{B_-}(0,x) \\\\\n",
    "\\mathbf{E}(0,x) &= \\mathbf{E_a}(0,x) + \\mathbf{E_+}(0,x) + \\mathbf{E_-}(0,x)\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Step 3: Compute v^i from B^i and E_i\n",
    "def ValenciavU_func_TW(**params):\n",
    "    x = rfm.xx_to_Cart[0]\n",
    "\n",
    "    B_aU = ixp.zerorank1(DIM=3)\n",
    "    E_aU = ixp.zerorank1(DIM=3)\n",
    "    B_pU = ixp.zerorank1(DIM=3)\n",
    "    E_pU = ixp.zerorank1(DIM=3)\n",
    "    B_mU = ixp.zerorank1(DIM=3)\n",
    "    E_mU = ixp.zerorank1(DIM=3)\n",
    "\n",
    "    B_aU[0] = sp.sympify(1)\n",
    "    B_aU[1] = noif.coord_leq_bound(x,0) * sp.sympify(1) + noif.coord_greater_bound(x,0) * sp.Rational(3,2)\n",
    "    B_aU[2] = sp.sympify(2)\n",
    "\n",
    "    E_aU[0] = noif.coord_leq_bound(x,0) * sp.sympify(-1) + noif.coord_greater_bound(x,0) * sp.Rational(-3,2)\n",
    "    E_aU[1] = sp.sympify(1)\n",
    "    E_aU[2] = sp.sympify(0)\n",
    "\n",
    "    B_pU[0] = sp.sympify(0)\n",
    "    B_pU[1] = noif.coord_leq_bound(x,0) * sp.sympify(0) + noif.coord_greater_bound(x,0) * sp.Rational(3,2)\n",
    "    B_pU[2] = noif.coord_leq_bound(x,0) * sp.sympify(0) + noif.coord_greater_bound(x,0) * sp.sympify(1)\n",
    "\n",
    "    E_pU[0] = sp.sympify(0)\n",
    "    E_pU[1] = noif.coord_leq_bound(x,0) * sp.sympify(0) + noif.coord_greater_bound(x,0) * sp.sympify(1)\n",
    "    E_pU[2] = noif.coord_leq_bound(x,0) * sp.sympify(0) + noif.coord_greater_bound(x,0) * sp.Rational(-3,2)\n",
    "\n",
    "    B_mU[0] = sp.sympify(0)\n",
    "    B_mU[1] = noif.coord_leq_bound(x,0) * sp.Rational(1,2)  + noif.coord_greater_bound(x,0) * sp.sympify(0)\n",
    "    B_mU[2] = noif.coord_leq_bound(x,0) * sp.Rational(3,2) + noif.coord_greater_bound(x,0) * sp.sympify(0)\n",
    "\n",
    "    E_mU[0] = sp.sympify(0)\n",
    "    E_mU[1] = noif.coord_leq_bound(x,0) * sp.Rational(-3,2) + noif.coord_greater_bound(x,0) * sp.sympify(0)\n",
    "    E_mU[2] = noif.coord_leq_bound(x,0) * sp.Rational(1,2)  + noif.coord_greater_bound(x,0) * sp.sympify(0)\n",
    "\n",
    "    BU = ixp.zerorank1(DIM=3)\n",
    "    EU = ixp.zerorank1(DIM=3)\n",
    "    for i in range(3):\n",
    "        BU[i] = B_aU[i] + B_pU[i] + B_mU[i]\n",
    "        EU[i] = E_aU[i] + E_pU[i] + E_mU[i]\n",
    "\n",
    "    # In flat space, ED and EU are identical, so we can still use this function.\n",
    "    return gfcf.compute_ValenciavU_from_ED_and_BU(EU, BU)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='code_validation'></a>\n",
    "\n",
    "# Step 4: Code Validation against `GiRaFFEfood_NRPy/GiRaFFEfood_NRPy` NRPy+ module \\[Back to [top](#toc)\\]\n",
    "$$\\label{code_validation}$$\n",
    "\n",
    "Here, as a code validation check, we verify agreement in the SymPy expressions for the `GiRaFFE` Aligned Rotator initial data equations  we intend to use between\n",
    "1. this tutorial and \n",
    "2. the NRPy+ [`GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_1D_tests.py`](../edit/GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_1D_tests.py) module.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consistency check between GiRaFFEfood_NRPy tutorial and NRPy+ module:\n",
      "ValenciavU0 is in agreement!\n",
      "AD0 is in agreement!\n",
      "ValenciavU1 is in agreement!\n",
      "AD1 is in agreement!\n",
      "ValenciavU2 is in agreement!\n",
      "AD2 is in agreement!\n"
     ]
    }
   ],
   "source": [
    "import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy as gf\n",
    "\n",
    "A_twD = gfcf.Axyz_func_Cartesian(Ax_TW,Ay_TW,Az_TW,stagger_enable = True,)\n",
    "Valenciav_twD = ValenciavU_func_TW()\n",
    "gf.GiRaFFEfood_NRPy_generate_initial_data(ID_type = \"ThreeWaves\", stagger_enable = True)\n",
    "\n",
    "def consistency_check(quantity1,quantity2,string):\n",
    "    if quantity1-quantity2==0:\n",
    "        print(string+\" is in agreement!\")\n",
    "    else:\n",
    "        print(string+\" does not agree!\")\n",
    "        sys.exit(1)\n",
    "\n",
    "print(\"Consistency check between GiRaFFEfood_NRPy tutorial and NRPy+ module:\")\n",
    "\n",
    "for i in range(3):\n",
    "    consistency_check(Valenciav_twD[i],gf.ValenciavU[i],\"ValenciavU\"+str(i))\n",
    "    consistency_check(A_twD[i],gf.AD[i],\"AD\"+str(i))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 5: 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-GiRaFFEfood_NRPy_1D_tests.pdf](Tutorial-GiRaFFEfood_NRPy_1D_tests.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-GiRaFFEfood_NRPy-Three_Waves.tex, and compiled LaTeX file\n",
      "    to PDF file Tutorial-GiRaFFEfood_NRPy-Three_Waves.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-GiRaFFEfood_NRPy-Three_Waves\",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
}