Notebook

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

This module provides another initial data option for `GiRaFFE`, drawn from Paschalidis and Shapiro .¶

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_Three_Waves.py ¶

Introduction:¶

Three waves:¶

This is a flat-spacetime test representing three Alfvén waves (one stationary, one left-going, and one right-going) with initial data \begin{align} A_x &= 0 \\ A_y &= 3.5x H(-x) + 3.0x H(x) \\ A_z &= y - 1.5x H(-x) - 3.0x H(x), \end{align} where $H(x)$ is the Heaviside function, which generates the magnetic field $$\mathbf{B}(0,x) = \mathbf{B_a}(0,x) + \mathbf{B_+}(0,x) + \mathbf{B_-}(0,x)$$ and uses the electric field $$\mathbf{E}(0,x) = \mathbf{E_a}(0,x) + \mathbf{E_+}(0,x) + \mathbf{E_-}(0,x),$$ 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 \begin{align} \mathbf{B_a}(0,x) &= \left \{ \begin{array}{lll} (1.0,1.0,2.0) & \mbox{if} & x<0 \\ (1.0,1.5,2.0) & \mbox{if} & x>0 \end{array} \right. , \\ \mathbf{E_a}(0,x) &= \left \{ \begin{array}{lll} (-1.0,1.0,0.0) & \mbox{if} & x<0 \\ (-1.5,1.0,0.0) & \mbox{if} & x>0 \end{array} \right. , \\ \mathbf{B_+}(0,x) &= \left \{ \begin{array}{lll} (0.0,0.0,0.0) & \mbox{if} & x<0 \\ (0.0,1.5,1.0) & \mbox{if} & x>0 \end{array} \right. , \\ \mathbf{E_+}(0,x) &= \left \{ \begin{array}{lll} (0.0,0.0,0.0) & \mbox{if} & x<0 \\ (0.0,1.0,-1.5) & \mbox{if} & x>0 \end{array} \right. , \\ \mathbf{B_-}(0,x) &= \left \{ \begin{array}{lll} (0.0,0.5,1.5) & \mbox{if} & x<0 \\ (0.0,0.0,0.0) & \mbox{if} & x>0 \end{array} \right. , \\ \mathbf{E_-}(0,x) &= \left \{ \begin{array}{lll} (0.0,-1.5,0.5) & \mbox{if} & x<0 \\ (0.0,0.0,0.0) & \mbox{if} & x>0 \end{array} \right. . \\ \end{align}

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

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 &= 3.5x H(-x) + 3.0x H(x) \\ A_z &= y - 1.5x H(-x) - 3.0x H(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. We thus get $$H(x) = \frac{\max(0,x)}{x}.$$ 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 for details on how this works). This is, conveniently, the exact implementation of the coord_greater_bound() function!

\begin{align} A_x &= 0 \\ A_y &= 3.5x H(-x) + 3.0x H(x) \\ A_z &= y - 1.5x H(-x) - 3.0x H(x), \end{align}

In [2]:

import Min_Max_and_Piecewise_Expressions as noif

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

def Ay_TW(x,y,z, **params):
    return sp.Rational(7,2)*x*noif.coord_greater_bound(-x,0) + sp.sympify(3)*x*noif.coord_greater_bound(x,0)

def Az_TW(x,y,z, **params):
    return y-sp.Rational(3,2)*x*noif.coord_greater_bound(-x,0) - sp.sympify(3)*x*noif.coord_greater_bound(x,0)

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

$$\label{set_vi}$$

First, we will set the three individual waves; we change all $<$ to $\leq$ to avoid unintented behavior at $x=0$: \begin{align} \mathbf{B_a}(0,x) &= \left \{ \begin{array}{lll} (1.0,1.0,2.0) & \mbox{if} & x \leq 0 \\ (1.0,1.5,2.0) & \mbox{if} & x>0 \end{array} \right. , \\ \mathbf{E_a}(0,x) &= \left \{ \begin{array}{lll} (-1.0,1.0,0.0) & \mbox{if} & x \leq 0 \\ (-1.5,1.0,0.0) & \mbox{if} & x>0 \end{array} \right. , \\ \mathbf{B_+}(0,x) &= \left \{ \begin{array}{lll} (0.0,0.0,0.0) & \mbox{if} & x \leq 0 \\ (0.0,1.5,1.0) & \mbox{if} & x>0 \end{array} \right. , \\ \mathbf{E_+}(0,x) &= \left \{ \begin{array}{lll} (0.0,0.0,0.0) & \mbox{if} & x \leq 0 \\ (0.0,1.0,-1.5) & \mbox{if} & x>0 \end{array} \right. , \\ \mathbf{B_-}(0,x) &= \left \{ \begin{array}{lll} (0.0,0.5,1.5) & \mbox{if} & x \leq 0 \\ (0.0,0.0,0.0) & \mbox{if} & x>0 \end{array} \right. , \\ \mathbf{E_-}(0,x) &= \left \{ \begin{array}{lll} (0.0,-1.5,0.5) & \mbox{if} & x \leq 0 \\ (0.0,0.0,0.0) & \mbox{if} & x>0 \end{array} \right. . \\ \end{align}

Then, we can obtain the total expressions for the magnetic and electric fields by simply adding the three waves together: \begin{align} \mathbf{B}(0,x) &= \mathbf{B_a}(0,x) + \mathbf{B_+}(0,x) + \mathbf{B_-}(0,x) \\ \mathbf{E}(0,x) &= \mathbf{E_a}(0,x) + \mathbf{E_+}(0,x) + \mathbf{E_-}(0,x) \end{align}

In [3]:

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

    B_aU = ixp.zerorank1(DIM=3)
    E_aU = ixp.zerorank1(DIM=3)
    B_pU = ixp.zerorank1(DIM=3)
    E_pU = ixp.zerorank1(DIM=3)
    B_mU = ixp.zerorank1(DIM=3)
    E_mU = ixp.zerorank1(DIM=3)

    B_aU[0] = sp.sympify(1)
    B_aU[1] = noif.coord_leq_bound(x,0) * sp.sympify(1) + noif.coord_greater_bound(x,0) * sp.Rational(3,2)
    B_aU[2] = sp.sympify(2)

    E_aU[0] = noif.coord_leq_bound(x,0) * sp.sympify(-1) + noif.coord_greater_bound(x,0) * sp.Rational(-3,2)
    E_aU[1] = sp.sympify(1)
    E_aU[2] = sp.sympify(0)

    B_pU[0] = sp.sympify(0)
    B_pU[1] = noif.coord_leq_bound(x,0) * sp.sympify(0) + noif.coord_greater_bound(x,0) * sp.Rational(3,2)
    B_pU[2] = noif.coord_leq_bound(x,0) * sp.sympify(0) + noif.coord_greater_bound(x,0) * sp.sympify(1)

    E_pU[0] = sp.sympify(0)
    E_pU[1] = noif.coord_leq_bound(x,0) * sp.sympify(0) + noif.coord_greater_bound(x,0) * sp.sympify(1)
    E_pU[2] = noif.coord_leq_bound(x,0) * sp.sympify(0) + noif.coord_greater_bound(x,0) * sp.Rational(-3,2)

    B_mU[0] = sp.sympify(0)
    B_mU[1] = noif.coord_leq_bound(x,0) * sp.Rational(1,2)  + noif.coord_greater_bound(x,0) * sp.sympify(0)
    B_mU[2] = noif.coord_leq_bound(x,0) * sp.Rational(3,2) + noif.coord_greater_bound(x,0) * sp.sympify(0)

    E_mU[0] = sp.sympify(0)
    E_mU[1] = noif.coord_leq_bound(x,0) * sp.Rational(-3,2) + noif.coord_greater_bound(x,0) * sp.sympify(0)
    E_mU[2] = noif.coord_leq_bound(x,0) * sp.Rational(1,2)  + noif.coord_greater_bound(x,0) * sp.sympify(0)

    BU = ixp.zerorank1(DIM=3)
    EU = ixp.zerorank1(DIM=3)
    for i in range(3):
        BU[i] = B_aU[i] + B_pU[i] + B_mU[i]
        EU[i] = E_aU[i] + E_pU[i] + E_mU[i]

    # 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 [4]:

import GiRaFFEfood_NRPy.GiRaFFEfood_NRPy as gf

A_twD = gfcf.Axyz_func_Cartesian(Ax_TW,Ay_TW,Az_TW,stagger_enable = True,)
Valenciav_twD = ValenciavU_func_TW()
gf.GiRaFFEfood_NRPy_generate_initial_data(ID_type = "ThreeWaves", 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_twD[i],gf.ValenciavU[i],"ValenciavU"+str(i))
    consistency_check(A_twD[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 [5]:

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-GiRaFFEfood_NRPy-Three_Waves")

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

1D Three Wave GiRaFFEfood Initial Data for GiRaFFE¶

This module provides another initial data option for GiRaFFE, drawn from Paschalidis and Shapiro .¶

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