#!/usr/bin/env python
# coding: utf-8

# <script async src="https://www.googletagmanager.com/gtag/js?id=UA-59152712-8"></script>
# <script>
#   window.dataLayer = window.dataLayer || [];
#   function gtag(){dataLayer.push(arguments);}
#   gtag('js', new Date());
# 
#   gtag('config', 'UA-59152712-8');
# </script>
# 
# # 1D Alfven Wave `GiRaFFEfood` Initial Data for `GiRaFFE`
# 
# ## This module provides another initial data option for `GiRaFFE`, drawn from [Paschalidis and Shapiro](https://arxiv.org/abs/1310.3274) .
# 
# **Notebook Status:** <font color='green'><b> Validated </b></font>
# 
# **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).
# 
# ### NRPy+ Source Code for this module: [GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Alfven_Wave.py](../../edit/in_progress/GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_Alfven_Wave.py)
# 
# ## Introduction:
# 
# ### Alfv&eacute;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-GiRaFFEfood_NRPy.ipynb) tutorial notebook for more general detail on how this is used.
# 

# <a id='toc'></a>
# 
# # Table of Contents:
# $$\label{toc}$$
# 
# This notebook is organized as follows
# 
# 1. [Step 1](#initializenrpy): Import core NRPy+ modules and set NRPy+ parameters
# 1. [Step 2](#set_a_i): Set the vector $A_i$
# 1. [Step 3](#set_vi): Calculate $v^i$ from $B^i$ and $E_i$
# 1. [Step 4](#code_validation): Code Validation against `GiRaFFEfood_NRPy.GiRaFFEfood_NRPy` NRPy+ Module
# 1. [Step 5](#latex_pdf_output): Output this notebook to $\LaTeX$-formatted PDF file

# <a id='initializenrpy'></a>
# 
# # Step 1: Import core NRPy+ modules and set NRPy+ parameters \[Back to [top](#toc)\]
# $$\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"


# ##### <a id='set_a_i'></a>
# 
# # Step 2: Set the vector $A_i$ \[Back to [top](#toc)\]
# $$\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


# <a id='set_vi'></a>
# 
# # Step 3: Calculate $v^i$ from $B^i$ and $E_i$ \[Back to [top](#toc)\]
# $$\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)


# <a id='code_validation'></a>
# 
# # Step 4: Code Validation against `GiRaFFEfood_NRPy.GiRaFFEfood_NRPy` NRPy+ Module \[Back to [top](#toc)\]
# $$\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
# 1. this tutorial and 
# 2. the NRPy+ [`GiRaFFEfood_NRPy/GiRaFFEfood_NRPy_1D_tests.py`](../edit/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))


# <a id='latex_pdf_output'></a>
# 
# # Step 5: Output this notebook to $\LaTeX$-formatted PDF file \[Back to [top](#toc)\]
# $$\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](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")