{
"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é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
}