Notebook

1D Alfven Wave `GiRaFFEfood` Initial Data for `GiRaFFE`¶

This module provides another initial data option for `GiRaFFE`, drawn from this paper .¶

Notebook Status: Validated

Validation Notes: This tutorial notebook has been confirmed to be self-consistent with its corresponding NRPy+ module, as documented below. The initial data has validated against the original GiRaFFE, as documented here.

NRPy+ Source Code for this module: GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Alfven_Wave.py ¶

Introduction:¶

Alfvén Wave:¶

This is a flat-spacetime test with initial data \begin{align} A_x &= 0 \\ A_y &= \left \{ \begin{array}{lll}\gamma_\mu x - 0.015 & \mbox{if} & x \leq -0.1/\gamma_\mu \\ 1.15 \gamma_\mu x - 0.03g(x) & \mbox{if} & -0.1/\gamma_\mu \leq x \leq 0.1/\gamma_\mu \\ 1.3 \gamma_\mu x - 0.015 & \mbox{if} & x \geq 0.1/\gamma_\mu \end{array} \right. , \\ A_z &= \ y - \gamma_\mu (1-\mu)x , \end{align} which generates the magnetic field in the wave frame, \begin{align} B'^{x'}(x') = &\ 1.0,\ B'^y(x') = 1.0, \\ B'^z(x') = &\ \left \{ \begin{array}{lll} 1.0 & \mbox{if} & x' \leq -0.1/\gamma_\mu \\ 1.0+0.15 f(x') & \mbox{if} & -0.1/\gamma_\mu \leq x' \leq 0.1/\gamma_\mu \\ 1.3 & \mbox{if} & x' \geq 0.1/\gamma_\mu \end{array} \right. . \end{align} The electric field in the wave frame is then given by $$E'^{x'}(x') = -B'^z(0,x') \ \ , \ \ E'^y(x') = 0.0 \ \ , \ \ E'^z(x') = 1.0 .$$

These are converted to the grid frame by \begin{align} B^x(0,x) = &\ B'^{x'}(\gamma_\mu x) , \\ B^y(0,x) = &\ \gamma_\mu [ B'^y(\gamma_\mu x) - \mu E'^z(\gamma_\mu x) ] , \\ B^z(0,x) = &\ \gamma_\mu [ B'^z(\gamma_\mu x) + \mu E'^y(\gamma_\mu x) ] , \end{align} and \begin{align} E^x(0,x) = &\ E'^{x'}(\gamma_\mu x) , \\ E^y(0,x) = &\ \gamma_\mu [ E'^y(\gamma_\mu x) + \mu B'^z(\gamma_\mu x) ] ,\\ E^z(0,x) = &\ \gamma_\mu [ E'^z(\gamma_\mu x) - \mu B'^y(\gamma_\mu x) ], \end{align} 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$.

For the eventual purpose of testing convergence, any quantity $Q$ evolves as $Q(t,x) = Q(0,x-\mu t)$

See the Tutorial-GiRaFFEfood_NRPy tutorial notebook for more general detail on how this is used.

Table of Contents:¶

$$\label{toc}$$

This notebook is organized as follows

Step 1: Import core NRPy+ modules and set NRPy+ parameters
Step 2: Set the vector $A_i$
Step 3: Calculate $v^i$ from $B^i$ and $E_i$
Step 4: Code Validation against GiRaFFEfood_NRPy.GiRaFFEfood_NRPy NRPy+ Module
Step 5: Output this notebook to $\LaTeX$-formatted PDF file

Step 1: Import core NRPy+ modules and set NRPy+ parameters [Back to top]¶

$$\label{initializenrpy}$$

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.

In [1]:

# Step 0: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os,sys
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)

# Step 0.a: Import the NRPy+ core modules and set the reference metric to Cartesian
import sympy as sp               # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par   # NRPy+: Parameter interface
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy_Common_Functions as gfcf # Some useful functions for GiRaFFE initial data.

import reference_metric as rfm   # NRPy+: Reference metric support
par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()

# Step 1a: Set commonly used parameters.
thismodule = "GiRaFFEfood_NRPy_1D"

¶

Step 2: Set the vector $A_i$ [Back to top]¶

$$\label{set_a_i}$$

The vector potential is given as \begin{align} A_x &= 0 \\ A_y &= \left \{ \begin{array}{lll}\gamma_\mu x - 0.015 & \mbox{if} & x \leq -0.1/\gamma_\mu \\ 1.15 \gamma_\mu x - 0.03g(x) & \mbox{if} & -0.1/\gamma_\mu \leq x \leq 0.1/\gamma_\mu \\ 1.3 \gamma_\mu x - 0.015 & \mbox{if} & x \geq 0.1/\gamma_\mu \end{array} \right. , \\ A_z &= y - \gamma_\mu (1-\mu)x . \end{align}

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.

In [2]:

mu_AW = par.Cparameters("REAL",thismodule,["mu_AW"], -0.5) # The wave speed
M_PI  = par.Cparameters("#define",thismodule,["M_PI"], "")

gammamu = sp.sympify(1)/sp.sqrt(sp.sympify(1)-mu_AW**2)
bound = sp.Rational(1,10)/gammamu
def g_AW(x):
    return sp.cos(sp.sympify(5)*M_PI*gammamu*x)/M_PI

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.

\begin{align} A_x &= 0 \\ A_y &= \left \{ \begin{array}{lll}\gamma_\mu x - 0.015 & \mbox{if} & x \leq -0.1/\gamma_\mu \\ 1.15 \gamma_\mu x - 0.03g(x) & \mbox{if} & -0.1/\gamma_\mu < x \leq 0.1/\gamma_\mu \\ 1.3 \gamma_\mu x - 0.015 & \mbox{if} & x > 0.1/\gamma_\mu \end{array} \right. , \\ A_z &= y - \gamma_\mu (1-\mu)x . \end{align}

In [3]:

import Min_Max_and_Piecewise_Expressions as noif

def Ax_AW(x,y,z, **params):
    return sp.sympify(0)

def Ay_AW(x,y,z, **params):
    # \gamma_\mu x - 0.015 if x <= -0.1/\gamma_\mu
    # 1.15 \gamma_\mu x - 0.03g(x) if -0.1/\gamma_\mu < x <= 0.1/\gamma_\mu
    # 1.3 \gamma_\mu x - 0.015 if x > 0.1/\gamma_\mu
    Ayleft = gammamu*x - sp.Rational(15,1000)
    Aycenter = sp.Rational(115,100)*gammamu*x - sp.Rational(3,100)*g_AW(x)
    Ayright = sp.Rational(13,10)*gammamu*x - sp.Rational(15,1000)

    out = noif.coord_leq_bound(x,-bound)*Ayleft\
         +noif.coord_greater_bound(x,-bound)*noif.coord_leq_bound(x,bound)*Aycenter\
         +noif.coord_greater_bound(x,bound)*Ayright
    return out

def Az_AW(x,y,z, **params):
    # y - \gamma_\mu (1-\mu)x
    return y-gammamu*(sp.sympify(1)-mu_AW)*x

Step 3: Calculate $v^i$ from $B^i$ and $E_i$ [Back to top]¶

$$\label{set_vi}$$

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$.

In [4]:

def f_AW(x):
    xprime = gammamu*x
    return 1 + sp.sin(5*M_PI*xprime)

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: \begin{align} B'^{x'}(x') = &\ 1.0,\ B'^y(x') = 1.0, \\ B'^z(x') = &\ \left \{ \begin{array}{lll} 1.0 & \mbox{if} & x' \leq -0.1 \\ 1.0+0.15 f(x') & \mbox{if} & -0.1 < x' \leq 0.1 \\ 1.3 & \mbox{if} & x' > 0.1 \end{array} \right. . \end{align}

Then, we will set the electric field in the wave frame: \begin{align} E'^{x'}(x') &= -B'^z(0,x'), \\ E'^y(x') &= 0.0, \\ E'^z(x') &= 1.0 . \end{align}

Next, we must transform the fields into the grid frame. We'll do the magnetic fields first. \begin{align} B^x(0,x) = &\ B'^{x'}(\gamma_\mu x) , \\ B^y(0,x) = &\ \gamma_\mu [ B'^y(\gamma_\mu x) - \mu E'^z(\gamma_\mu x) ] , \\ B^z(0,x) = &\ \gamma_\mu [ B'^z(\gamma_\mu x) + \mu E'^y(\gamma_\mu x) ] , \end{align}

And finally the electric fields: \begin{align} E^x(0,x) = &\ E'^{x'}(\gamma_\mu x) , \\ E^y(0,x) = &\ \gamma_\mu [ E'^y(\gamma_\mu x) + \mu B'^z(\gamma_\mu x) ] ,\\ E^z(0,x) = &\ \gamma_\mu [ E'^z(\gamma_\mu x) - \mu B'^y(\gamma_\mu x) ], \end{align}

In [5]:

#Step 3: Compute v^i from B^i and E_i
def ValenciavU_func_AW(**params):
    x = rfm.xx_to_Cart[0]

    Bzleft = sp.sympify(1)
    Bzcenter = sp.sympify(1) + sp.Rational(15,100)*f_AW(x)
    Bzright = sp.Rational(13,10)

    BpU = ixp.zerorank1()
    BpU[0] = sp.sympify(1)
    BpU[1] = sp.sympify(1)
    BpU[2] = noif.coord_leq_bound(x,-bound)*Bzleft\
            +noif.coord_greater_bound(x,-bound)*noif.coord_leq_bound(x,bound)*Bzcenter\
            +noif.coord_greater_bound(x,bound)*Bzright

    EpU = ixp.zerorank1()
    EpU[0] = -BpU[2]
    EpU[1] = sp.sympify(0)
    EpU[2] = sp.sympify(1)

    BU = ixp.zerorank1()
    BU[0] = BpU[0]
    BU[1] = gammamu*(BpU[1]-mu_AW*EpU[2])
    BU[2] = gammamu*(BpU[2]+mu_AW*EpU[1])

    EU = ixp.zerorank1()
    EU[0] = EpU[0]
    EU[1] = gammamu*(EpU[1]+mu_AW*BpU[2])
    EU[2] = gammamu*(EpU[2]-mu_AW*BpU[1])

    # In flat space, ED and EU are identical, so we can still use this function.
    return gfcf.compute_ValenciavU_from_ED_and_BU(EU, BU)

Step 4: Code Validation against `GiRaFFEfood_NRPy.GiRaFFEfood_NRPy` NRPy+ Module [Back to top]¶

$$\label{code_validation}$$

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

this tutorial and
the NRPy+ GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_1D_tests.py module.

In [6]:

import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy as gf

A_awD = gfcf.Axyz_func_Cartesian(Ax_AW,Ay_AW,Az_AW,stagger_enable = True,)
Valenciav_awD = ValenciavU_func_AW()
gf.GiRaFFEfood_NRPy_generate_initial_data(ID_type = "AlfvenWave", stagger_enable = True)

def consistency_check(quantity1,quantity2,string):
    if quantity1-quantity2==0:
        print(string+" is in agreement!")
    else:
        print(string+" does not agree!")
        sys.exit(1)

print("Consistency check between GiRaFFEfood_NRPy tutorial and NRPy+ module:")

for i in range(3):
    consistency_check(Valenciav_awD[i],gf.ValenciavU[i],"ValenciavU"+str(i))
    consistency_check(A_awD[i],gf.AD[i],"AD"+str(i))

Consistency check between GiRaFFEfood_NRPy tutorial and NRPy+ module:
ValenciavU0 is in agreement!
AD0 is in agreement!
ValenciavU1 is in agreement!
AD1 is in agreement!
ValenciavU2 is in agreement!
AD2 is in agreement!

Step 5: Output this notebook to $\LaTeX$-formatted PDF file [Back to top]¶

$$\label{latex_pdf_output}$$

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 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.)

In [7]:

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-GiRaFFEfood_NRPy-Alfven_Wave",location_of_template_file=os.path.join(".."))

Created Tutorial-GiRaFFEfood_NRPy-Alfven_Wave.tex, and compiled LaTeX file
    to PDF file Tutorial-GiRaFFEfood_NRPy-Alfven_Wave.pdf

1D Alfven Wave GiRaFFEfood Initial Data for GiRaFFE¶

This module provides another initial data option for GiRaFFE, drawn from this paper .¶

NRPy+ Source Code for this module: GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Alfven_Wave.py¶