{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# 1D Alfven Wave `GiRaFFEfood` Initial Data for `GiRaFFE`\n", "\n", "## This module provides another initial data option for `GiRaFFE`, drawn from [Paschalidis and Shapiro](https://arxiv.org/abs/1310.3274) .\n", "\n", "**Notebook Status:** Validated \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": [ "\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": [ "\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": [ "##### \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": [ "\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": [ "\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": [ "\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\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.10.11" } }, "nbformat": 4, "nbformat_minor": 4 }