{
"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 Alfven 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_Alfven_Wave.py](../../edit/in_progress/GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Alfven_Wave.py)\n",
"\n",
"## Introduction:\n",
"\n",
"### Alfvén Wave:\n",
"\n",
" This is a flat-spacetime test with initial data \n",
"\\begin{align}\n",
"A_x &= 0 \\\\\n",
"A_y &= \\left \\{ \\begin{array}{lll}\\gamma_\\mu x - 0.015 & \\mbox{if} & x \\leq -0.1/\\gamma_\\mu \\\\\n",
"1.15 \\gamma_\\mu x - 0.03g(x) & \\mbox{if} & -0.1/\\gamma_\\mu \\leq x \\leq 0.1/\\gamma_\\mu \\\\ \n",
"1.3 \\gamma_\\mu x - 0.015 & \\mbox{if} & x \\geq 0.1/\\gamma_\\mu \\end{array} \\right. , \\\\\n",
"A_z &= \\ y - \\gamma_\\mu (1-\\mu)x ,\n",
"\\end{align}\n",
"which generates the magnetic field in the wave frame,\n",
"\\begin{align}\n",
"B'^{x'}(x') = &\\ 1.0,\\ B'^y(x') = 1.0, \\\\\n",
"B'^z(x') = &\\ \\left \\{ \\begin{array}{lll} 1.0 & \\mbox{if} & x' \\leq -0.1/\\gamma_\\mu \\\\\n",
"\t\t\t\t1.0+0.15 f(x') & \\mbox{if} & -0.1/\\gamma_\\mu \\leq x' \\leq 0.1/\\gamma_\\mu \\\\\n",
"\t\t\t\t1.3 & \\mbox{if} & x' \\geq 0.1/\\gamma_\\mu \\end{array} \\right. .\n",
"\\end{align}\n",
"The electric field in the wave frame is then given by\n",
"$$E'^{x'}(x') = -B'^z(0,x') \\ \\ , \\ \\ E'^y(x') = 0.0 \\ \\ , \\ \\ E'^z(x') = 1.0 .$$\n",
"\n",
"These are converted to the grid frame by \n",
"\\begin{align}\n",
" B^x(0,x) = &\\ B'^{x'}(\\gamma_\\mu x) , \\\\\n",
" B^y(0,x) = &\\ \\gamma_\\mu [ B'^y(\\gamma_\\mu x) - \\mu E'^z(\\gamma_\\mu x) ] , \\\\ \n",
" B^z(0,x) = &\\ \\gamma_\\mu [ B'^z(\\gamma_\\mu x) + \\mu E'^y(\\gamma_\\mu x) ] , \n",
"\\end{align}\n",
"and\n",
"\\begin{align}\n",
" E^x(0,x) = &\\ E'^{x'}(\\gamma_\\mu x) , \\\\ \n",
" E^y(0,x) = &\\ \\gamma_\\mu [ E'^y(\\gamma_\\mu x) + \\mu B'^z(\\gamma_\\mu x) ] ,\\\\ \n",
" E^z(0,x) = &\\ \\gamma_\\mu [ E'^z(\\gamma_\\mu x) - \\mu B'^y(\\gamma_\\mu x) ],\n",
"\\end{align}\n",
"and the velocity is given by $$\\mathbf{v} = \\frac{\\mathbf{E} \\times \\mathbf{B}}{B^2}$$ in flat spacetime. Additionally, $f(x)=1+\\sin (5\\pi x)$, $-1<\\mu<1$ is the wave speed relative to the grid frame and $\\gamma_\\mu = (1-\\mu^2)^{-1/2}$, and $g(x) = \\cos (5\\pi \\gamma_\\mu x)/\\pi$.\n",
"\n",
"For the eventual purpose of testing convergence, any quantity $Q$ evolves as $Q(t,x) = Q(0,x-\\mu 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 &= \\left \\{ \\begin{array}{lll}\\gamma_\\mu x - 0.015 & \\mbox{if} & x \\leq -0.1/\\gamma_\\mu \\\\\n",
"1.15 \\gamma_\\mu x - 0.03g(x) & \\mbox{if} & -0.1/\\gamma_\\mu \\leq x \\leq 0.1/\\gamma_\\mu \\\\ \n",
"1.3 \\gamma_\\mu x - 0.015 & \\mbox{if} & x \\geq 0.1/\\gamma_\\mu \\end{array} \\right. , \\\\\n",
"A_z &= y - \\gamma_\\mu (1-\\mu)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`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"mu_AW = par.Cparameters(\"REAL\",thismodule,[\"mu_AW\"], -0.5) # The wave speed\n",
"M_PI = par.Cparameters(\"#define\",thismodule,[\"M_PI\"], \"\")\n",
"\n",
"gammamu = sp.sympify(1)/sp.sqrt(sp.sympify(1)-mu_AW**2)\n",
"bound = sp.Rational(1,10)/gammamu\n",
"def g_AW(x):\n",
" return sp.cos(sp.sympify(5)*M_PI*gammamu*x)/M_PI"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we can define the vector potential. We will rewrite $A_y$ to make use of the functions provided by `Min_Max_and_Piecewise_Expressions`. As shown below, we make sure that at each boundary, each $\\leq$ is paired with a $>$. (This choice is arbitrary, we could just as easily choose $<$ and $\\geq$.) This does not change the data since the function is continuous. However, it is necessary for the functions in `Min_Max_and_Piecewise_Expressions` to output the correct results.\n",
"\n",
"\\begin{align}\n",
"A_x &= 0 \\\\\n",
"A_y &= \\left \\{ \\begin{array}{lll}\\gamma_\\mu x - 0.015 & \\mbox{if} & x \\leq -0.1/\\gamma_\\mu \\\\\n",
"1.15 \\gamma_\\mu x - 0.03g(x) & \\mbox{if} & -0.1/\\gamma_\\mu < x \\leq 0.1/\\gamma_\\mu \\\\ \n",
"1.3 \\gamma_\\mu x - 0.015 & \\mbox{if} & x > 0.1/\\gamma_\\mu \\end{array} \\right. , \\\\\n",
"A_z &= y - \\gamma_\\mu (1-\\mu)x .\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import Min_Max_and_Piecewise_Expressions as noif\n",
"\n",
"def Ax_AW(x,y,z, **params):\n",
" return sp.sympify(0)\n",
"\n",
"def Ay_AW(x,y,z, **params):\n",
" # \\gamma_\\mu x - 0.015 if x <= -0.1/\\gamma_\\mu\n",
" # 1.15 \\gamma_\\mu x - 0.03g(x) if -0.1/\\gamma_\\mu < x <= 0.1/\\gamma_\\mu\n",
" # 1.3 \\gamma_\\mu x - 0.015 if x > 0.1/\\gamma_\\mu\n",
" Ayleft = gammamu*x - sp.Rational(15,1000)\n",
" Aycenter = sp.Rational(115,100)*gammamu*x - sp.Rational(3,100)*g_AW(x)\n",
" Ayright = sp.Rational(13,10)*gammamu*x - sp.Rational(15,1000)\n",
"\n",
" out = noif.coord_leq_bound(x,-bound)*Ayleft\\\n",
" +noif.coord_greater_bound(x,-bound)*noif.coord_leq_bound(x,bound)*Aycenter\\\n",
" +noif.coord_greater_bound(x,bound)*Ayright\n",
" return out\n",
"\n",
"def Az_AW(x,y,z, **params):\n",
" # y - \\gamma_\\mu (1-\\mu)x\n",
" return y-gammamu*(sp.sympify(1)-mu_AW)*x"
]
},
{
"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",
"Now, we will set the magnetic and electric fields that we will need to define the initial velocities. First, we need to define $$f(x)=1+\\sin (5\\pi x);$$ note that in the definition of $B^i$, we need $f(x')$ where $x'=\\gamma_\\mu x$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def f_AW(x):\n",
" xprime = gammamu*x\n",
" return 1 + sp.sin(5*M_PI*xprime)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will first set the magnetic field in the wave frame, once again rewriting $B'^z(x')$ to be compatible with `Min_Max_and_Piecewise_Expressions`:\n",
"\\begin{align}\n",
"B'^{x'}(x') = &\\ 1.0,\\ B'^y(x') = 1.0, \\\\\n",
"B'^z(x') = &\\ \\left \\{ \\begin{array}{lll} 1.0 & \\mbox{if} & x' \\leq -0.1 \\\\\n",
"\t\t\t\t1.0+0.15 f(x') & \\mbox{if} & -0.1 < x' \\leq 0.1 \\\\\n",
"\t\t\t\t1.3 & \\mbox{if} & x' > 0.1 \\end{array} \\right. .\n",
"\\end{align}\n",
"\n",
"Then, we will set the electric field in the wave frame:\n",
"\\begin{align}\n",
"E'^{x'}(x') &= -B'^z(0,x'), \\\\ \n",
"E'^y(x') &= 0.0, \\\\ \n",
"E'^z(x') &= 1.0 .\n",
"\\end{align}\n",
"\n",
"Next, we must transform the fields into the grid frame. We'll do the magnetic fields first.\n",
"\\begin{align}\n",
" B^x(0,x) = &\\ B'^{x'}(\\gamma_\\mu x) , \\\\\n",
" B^y(0,x) = &\\ \\gamma_\\mu [ B'^y(\\gamma_\\mu x) - \\mu E'^z(\\gamma_\\mu x) ] , \\\\ \n",
" B^z(0,x) = &\\ \\gamma_\\mu [ B'^z(\\gamma_\\mu x) + \\mu E'^y(\\gamma_\\mu x) ] , \n",
"\\end{align}\n",
"\n",
"And finally the electric fields:\n",
"\\begin{align}\n",
" E^x(0,x) = &\\ E'^{x'}(\\gamma_\\mu x) , \\\\ \n",
" E^y(0,x) = &\\ \\gamma_\\mu [ E'^y(\\gamma_\\mu x) + \\mu B'^z(\\gamma_\\mu x) ] ,\\\\ \n",
" E^z(0,x) = &\\ \\gamma_\\mu [ E'^z(\\gamma_\\mu x) - \\mu B'^y(\\gamma_\\mu x) ],\n",
"\\end{align}\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#Step 3: Compute v^i from B^i and E_i\n",
"def ValenciavU_func_AW(**params):\n",
" x = rfm.xx_to_Cart[0]\n",
"\n",
" Bzleft = sp.sympify(1)\n",
" Bzcenter = sp.sympify(1) + sp.Rational(15,100)*f_AW(x)\n",
" Bzright = sp.Rational(13,10)\n",
"\n",
" BpU = ixp.zerorank1()\n",
" BpU[0] = sp.sympify(1)\n",
" BpU[1] = sp.sympify(1)\n",
" BpU[2] = noif.coord_leq_bound(x,-bound)*Bzleft\\\n",
" +noif.coord_greater_bound(x,-bound)*noif.coord_leq_bound(x,bound)*Bzcenter\\\n",
" +noif.coord_greater_bound(x,bound)*Bzright\n",
"\n",
" EpU = ixp.zerorank1()\n",
" EpU[0] = -BpU[2]\n",
" EpU[1] = sp.sympify(0)\n",
" EpU[2] = sp.sympify(1)\n",
"\n",
" BU = ixp.zerorank1()\n",
" BU[0] = BpU[0]\n",
" BU[1] = gammamu*(BpU[1]-mu_AW*EpU[2])\n",
" BU[2] = gammamu*(BpU[2]+mu_AW*EpU[1])\n",
"\n",
" EU = ixp.zerorank1()\n",
" EU[0] = EpU[0]\n",
" EU[1] = gammamu*(EpU[1]+mu_AW*BpU[2])\n",
" EU[2] = gammamu*(EpU[2]-mu_AW*BpU[1])\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": 6,
"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_awD = gfcf.Axyz_func_Cartesian(Ax_AW,Ay_AW,Az_AW,stagger_enable = True,)\n",
"Valenciav_awD = ValenciavU_func_AW()\n",
"gf.GiRaFFEfood_NRPy_generate_initial_data(ID_type = \"AlfvenWave\", 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_awD[i],gf.ValenciavU[i],\"ValenciavU\"+str(i))\n",
" consistency_check(A_awD[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": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-GiRaFFEfood_NRPy-Alfven_Wave.tex, and compiled LaTeX file\n",
" to PDF file Tutorial-GiRaFFEfood_NRPy-Alfven_Wave.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-Alfven_Wave\",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
}