


# The Riemann Solution on $\tilde{S}_i$ using HLLE
### Author: Patrick Nelson

This notebook documents the function from the original `GiRaFFE` that calculates the flux for $\tilde{S}_i$ according to the method of Harten, Lax, von Leer, and Einfeldt (HLLE), assuming that we have calculated the values of the velocity and magnetic field on the cell faces according to the piecewise-parabolic method (PPM) of [Colella and Woodward (1984)](https://crd.lbl.gov/assets/pubs_presos/AMCS/ANAG/A141984.pdf), modified for the case of GRFFE. 

**Notebook Status:** <font color=green><b> Validated against the original `GiRaFFE` </b></font>

**Validation Notes:** Everything should be written at this point, but we need to finish the rest of the algorithm before we can validate this.

The differential equations that `GiRaFFE` evolves are written in conservation form, and thus have two different terms that contribute to the time evolution of some quantity: the flux term and the source term. The PPM method is what the original `GiRaFFE` uses to handle the flux term; hopefully, using this instead of finite-differencing will fix some of the problems we've been having with `GiRaFFE_HO`.

In GRFFE, the evolution equation for the Poynting flux $\tilde{S}_i$ is given as 
$$
\boxed{\partial_t \tilde{S}_i + \underbrace{ \partial_j \left( \alpha \sqrt{\gamma} T^j_{{\rm EM} i} \right)}_{\rm Flux\ term} = \underbrace{\frac{1}{2} \alpha \sqrt{\gamma} T^{\mu \nu}_{\rm EM} \partial_i g_{\mu \nu}}_{\rm Source\ term}.}
$$
We can then see that, if we rewrite this, the right-hand side (RHS) describing the time evolution $\partial_t \tilde{S}_i$ consists of two terms: the flux term and the source term. The flux term in particular can be tricky, as it may be discontinuous due to shocks or other sharp features. This presents difficulties when we take finite-difference derivatives of that term, leading to the Gibbs phenomenon. So, we implement a different algorithm to take the derivative. 

It is the flux term $\alpha \sqrt{\gamma} T^j_{{\rm EM} i} = \alpha \sqrt{\gamma} g_{i \mu} T^{\mu i}_{\rm EM}$ itself that concerns us here, where $T^{\mu \nu}_{\rm EM} = b^2 u^\mu u^\nu + \frac{1}{2} b^2 g^{\mu \nu} - b^\mu b^\nu$; the following functions will compute this value so that we can easily take its derivative later. Having reconstructed the values of $v^i_{(n)}$ and $B^i$ on the cell faces, we will now compute the value of the flux of $\tilde{S}_i$ on each face. For each component of $\tilde{S}_i$ in each direction, we compute the flux as
$$
{\rm st\_j\_flux} = \frac{c_{\rm min} f_{\rm R} + c_{\rm max} f_{\rm L} - c_{\rm min} c_{\rm max} ({\rm st\_j\_r}-{\rm st\_j\_l})}{c_{\rm min} + c_{\rm max}},
$$
where 
$$
f = \alpha \sqrt{\gamma} \left( (\rho+b^2)(u^0 v^j) u_i + (P+\frac{1}{2}b^2) \delta_{ji} - b^j b_i \right)
$$
and
$$
{\rm st\_j} = \alpha \sqrt{\gamma} \left( (\rho+b^2) u^0 u_i - b^t b_i \right).
$$
Here, $i$ is direction in which we are computing the flux, and $j$ is the component of the momentum we are computing it for. Note that these two quantities are computed on both the left and right sides of the cell face. We will be able to draw heavily on the module [Tutorial-u0_smallb_Poynting-Cartesian.ipynb](Tutorial-u0_smallb_Poynting-Cartesian.ipynb) to compute $u^0$, $u^i$, and $b^\mu$, as well as the index-lowered forms of those vectors. Critically, these quantities depend on the Valencia 3-velocity $v^i_{(n)}$ and magnetic field $B^i$. We will not be using the normal gridfunctions for these, but rather the ones that we have previosly calculated on the left and right sides of the cell faces using the [Piecewise Parabolic Method](Tutorial-GiRaFFE_HO_Ccode_library-PPM.ipynb).

The speeds $c_\min$ and $c_\max$ are characteristic speeds that waves can travel through the plasma. In GRFFE, the expressions defining them reduce a function of only the metric quantities. $c_\min$ is the negative of the minimum amongst the speeds $c_-$ and $0$ and $c_\max$ is the maximum amongst the speeds $c_+$ and $0$. The speeds $c_\pm = \sqrt{b^2-4ac}$ must be calculated on both the left and right faces, where 
$$a = 1/\alpha^2,$$ 
$$b = 2 \beta^i / \alpha^2$$
and $$c = g^{ii} - (\beta^i)^2/\alpha^2.$$

Another point to consider is that since we are working on cell faces, not at the cell center, we can't use the normal metric values that we store. We will instead use the value of the metric interpolated onto the cell face, which in this tutorial, we will assume has been previously done. 

The algorithm for finite-volume methods in general is as follows: 

1. The Reconstruction Step - Piecewise Parabolic Method
    1. Within each cell, fit to a function that conserves the volume in that cell using information from the neighboring cells
        * For PPM, we will naturally use parabolas
    1. Use that fit to define the state at the left and right interface of each cell
    1. Apply a slope limiter to mitigate Gibbs phenomenon
1. Interpolate the value of the metric gridfunctions on the cell faces
1. **Solving the Riemann Problem - Harten, Lax, (This notebook, $\tilde{S}_i$ only)**
    1. **Use the left and right reconstructed states to calculate the unique state at boundary**
1. Use the unique state to estimate the derivative in the cell

<a id='toc'></a>

# Table of Contents
$$\label{toc}$$

This notebook is organized as follows

1. [Step 1](#prelim): Preliminaries
1. [Step 2](#s_i_flux): The $\tilde{S}_i$ function
    1. [Step 2.a](#define): The function declaration: inputs, outputs, and macros
    1. [Step 2.b](#speed_limit): Speed-limit the velocities
    1. [Step 2.c](#smallb): Computing $b^\mu$, the magnetic field in the comoving fluid frame
    1. [Step 2.d](#hydro_speed): GRFFE characteristic wave speeds
    1. [Step 2.e](#fluxes): Compute the HLLE fluxes
1. [Step 3](#code_validation): Code Validation against `GiRaFFE_HO.GiRaFFE_Higher_Order_v2` NRPy+ Module
1. [Step 4](#derive_speed): Complete Derivation of the Wave Speeds
1. [Step 5](#latex_pdf_output): Output this notebook to $\LaTeX$-formatted PDF file

<a id='prelim'></a>

# Step 1: Preliminaries \[Back to [top](#toc)\]
$$\label{prelim}$$

This first block of code just sets up a subdirectory within `GiRaFFE_standalone_Ccodes/` to which we will write the C code. We will also import the core NRPy+ functionality and register the needed gridfunctions. Doing so will let NRPy+ figure out where to read and write data from/to. 


```python
# 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)

from outputC import *            # NRPy+: Core C code output module
import finite_difference as fin  # NRPy+: Finite difference C code generation module
import NRPy_param_funcs as par   # NRPy+: Parameter interface
import grid as gri               # NRPy+: Functions having to do with numerical grids
import loop as lp                # NRPy+: Generate C code loops
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm   # NRPy+: Reference metric support
import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface
import shutil, os, sys           # Standard Python modules for multiplatform OS-level functions

outdir = "GiRaFFE_standalone_Ccodes/"
cmd.mkdir(outdir)
outdir += "PPM/"
cmd.mkdir(outdir)

par.set_parval_from_str("grid::DIM", 3)
DIM = par.parval_from_str("grid::DIM")

thismodule = "GiRaFFE_NRPy_Ccode_library-Stilde-flux"

```

We will also create the gridfunctions within NRPy+ so that it knows how to access them. Note that we declare different gridfunctions than normal, since we want to use the values of those gridfunctions interpolated on the cell faces.


```python
# These are the standard gridfunctions we've used before.
#ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU",DIM=3)
#gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01")
#betaU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","betaU")
#alpha = gri.register_gridfunctions("AUXEVOL",["alpha"])
#AD = ixp.register_gridfunctions_for_single_rank1("EVOL","AD",DIM=3)
#BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU",DIM=3)

# We will pass values of the gridfunction on the cell faces into the function. This requires us
# to declare them as C parameters in NRPy+. We will denote this with the _face infix/suffix.
alpha_face,gammadet_face = gri.register_gridfunctions("AUXEVOL",["alpha_face","gammadet_face"])
gamma_faceDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gamma_faceDD","sym01")
gamma_faceUU = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gamma_faceUU","sym01")
beta_faceU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","beta_faceU")

# We'll need some more gridfunctions, now, to represent the reconstructions of BU and ValenciavU
# on the right and left faces
Valenciav_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_rU",DIM=3)
B_rU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_rU",DIM=3)
Valenciav_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Valenciav_lU",DIM=3)
B_lU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","B_lU",DIM=3)

# ...and some more for the fluxes we calculate here. These three gridfunctions will each store
# the momentum flux of one component of StildeD in one direction; we'll be able to reuse them
# as we loop over each direction, reducing our memory costs.
Stilde_fluxD = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","Stilde_fluxD",DIM=3)
```

<a id='s_i_flux'></a>

# Step 2: The $\tilde{S}_i$ function \[Back to [top](#toc)\]
$$\label{s_i_flux}$$

<a id='define'></a>

## Step 2.a: The function declaration: inputs, outputs, and macros \[Back to [top](#toc)\]
$$\label{define}$$

First, we give the function definition. This function will take as parameters the integers `i0`, `i1`, and `i2`; these indicate a point on the grid. The function also needs the array of gridfunctions (which we use here for the reconstructed face values (see [the previous tutorial](Tutorial-GiRaFFE_HO_Ccode_library-PPM.ipynb) and the interpolated metric quantities on the cell interfaces). We also calculate quantities related to the conformal factor $\psi$ (including the square root of the determinant of the metric $\sqrt{\gamma}$) and the lapse $\alpha$; these lines will be handled in Python so the NRPy+ can optimize them. 


```python
# This product occurs very frequently, so we'll set it now. Conveniently,
# we store the metric determinant gamma directly. 
alpha_sqrt_gamma = alpha_face*sp.sqrt(gammadet_face)

# We'll also compute some powers of the conformal factor psi:
# Since psi^6 = sqrt(gammadet) by definition,
psi = sp.sqrt(gammadet_face)**(1.0/6.0)
psim4 = psi**(-4.0)
```

<a id='speed_limit'></a>

## Step 2.b: Speed-limit the velocities \[Back to [top](#toc)\]
$$\label{speed_limit}$$

Next, we compute the speed-limited velocities on the faces. This function is found in `inlined_functions.C`, however, we will replace it with our own that does the same thing (and is already compatible with the rest of our code). The original function also counts how often the speed had to be limited. We also introduce a variation of the previous speed limiter that allows the process to be incorporated directly into the final expressions that we will pass to `FD_outputC()`.

<a id='u0_no_if'></a>

### Step 2.b.i: Streamlining $u^0$ \[Back to [top](#toc)\]
$$\label{u0_no_if}$$

When we calculate $u^0$, we first check the Valencia 3-velocity and limit its speed if necessary. This is done as explained in the tutorial found [here](Tutorial-u0_smallb_Poynting-Cartesian.ipynb):

>Then our algorithm for computing $u^0$ is as follows:
>
>If
>$$R=\gamma_{ij}v^i_{(n)}v^j_{(n)}>1 - \frac{1}{\Gamma_{\rm max}^2},$$ 
then adjust the 3-velocity $v^i$ as follows:
>
>$$v^i_{(n)} = \sqrt{\frac{1 - \frac{1}{\Gamma_{\rm max}^2}}{R}}v^i_{(n)}.$$
>
>After this rescaling, we are then guaranteed that if $R$ is recomputed, it will be set to its ceiling value $R=R_{\rm max} = 1 - \frac{1}{\Gamma_{\rm max}^2}$.
>
>Then, regardless of whether the ceiling on $R$ was applied, $u^0$ can be safely computed via
>
>$$
u^0 = \frac{1}{\alpha \sqrt{1-R}}.
$$

However, we can once again eliminate the `if` statement used here to streamline our code generation. Note that when we violate the speed limit, we rescale by a factor $\sqrt{R_{\rm max}/R}$, where $R_{\rm max}=1 - \frac{1}{\Gamma_{\rm max}}$. It thus follows that, by letting $R_\max$ instead equal $\min(R,1 - \frac{1}{\Gamma_{\rm max}})$, we can *always* rescale by that same factor. In the case where we previously would not have rescaled, $\sqrt{R_{\rm max}/R}=1$, and we have not changed anything we didn't want to. 



```python
# This function rewrites the original version from Tutorial-u0_smallb_Poynting-Cartesian.ipynb
# without the if statement.
def compute_u0_noif(gammaDD,alpha,ValenciavU):
    # Inputs:  Metric gammaDD, lapse alpha, Valencia 3-velocity ValenciavU
    # Outputs: u^0, speed-limited ValenciavU
    
    # R = gamma_{ij} v^i v^j
    R = par.Cparameters("#define",thismodule,"TINYDOUBLE",1e-100)
    # We initialize R to TINYDOUBLE instead of 0 to avoid a 0/0 below,
    #    in the case that ValenciavU == 0 for all Valencia 3-velocities.
    # "Those tiny *doubles* make me warm all over
    #  with a feeling that I'm gonna love you till the end of time."
    #    - Adapted from Connie Francis' "Tiny Bubbles"
    for i in range(DIM):
        for j in range(DIM):
            R += gammaDD[i][j]*ValenciavU[i]*ValenciavU[j]

    # The default value isn't terribly important here, since we'll overwrite in the main C code
    GAMMA_SPEED_LIMIT = par.Cparameters("REAL",thismodule,"GAMMA_SPEED_LIMIT",10.0) # Default value based on
                                                                                    # IllinoisGRMHD.
                                                                                    # GiRaFFE default = 2000.0
    # Rmax = 1 - 1/Gamma_{\max}^2
    Rmax = 1 - 1/(GAMMA_SPEED_LIMIT*GAMMA_SPEED_LIMIT)
    # Now, we set Rmax = min(Rmax,R):
    # If Rmax>R, then Rmax = 0.5*(Rmax+R-Rmax+R) = R
    # If R>Rmax, then Rmax = 0.5*(Rmax+R+Rmax-R) = Rmax
    # If R==TINYDOUBLE, then Rmax = 0.5*(Rmax-Rmax) = 0, since, e.g., 10 +/- 1e-100 = 10 exactly in double precision    Rmax =  sp.Rational(1,2)*(Rmax+R-nrpyAbs(Rmax-R))
    Rmax =  sp.Rational(1,2)*(Rmax+R-nrpyAbs(Rmax-R))

    # With our rescaled Rmax, v^i = sqrt{Rmax/R} v^i
    rescaledValenciavU = ixp.zerorank1()
    for i in range(DIM):
        # If R == TINYDOUBLE, then Rmax=0,R=1e-100, and we get Rmax/R=0.
        #   If your velocities are of order 1e-100 and this is physically
        #   meaningful, there must be something wrong with your unit conversion.
        rescaledValenciavU[i] = ValenciavU[i]*sp.sqrt(Rmax/R)

    # We stick with the "rescaled" version since we redefine Rmax as detailed above
    # u^0 = 1/(alpha-sqrt(1-R_{\max}))
    rescaledu0 = 1/(alpha*sp.sqrt(1-Rmax))
    return rescaledValenciavU,rescaledu0

```

<a id='smallb'></a>

## Step 2.c: Computing $b^\mu$, the magnetic field in the comoving fluid frame \[Back to [top](#toc)\]
$$\label{smallb}$$

Now, we will compute the magnetic field in the comoving fluid frame $b^\mu$ , as defined in [Gammie's paper](https://arxiv.org/pdf/astro-ph/0301509.pdf).  We will also lower the index of the four velocities on the faces $u_\mu$.

Conveniently enough, the NRPy+ code that does this has already been written, and is found in `u0_smallb_Poynting__Cartesian.py`, which is documented [here](Tutorial-u0_smallb_Poynting-Cartesian.ipynb). We will call the function at least twice, once each for left and right side of the face. If we use the old version of the speed-limiting algorithm, we will need to call it a third time so that we generate a more generic version of `computeu0_Cfunction.h`.


```python
# We already did something like this in GiRaFFE_Higher_Order_v2; we'll 
# paste that code here, but using the face variables. This assumes that
# we've already interpolated the metric quantities on the faces.
import u0_smallb_Poynting__Cartesian.u0_smallb_Poynting__Cartesian as u0b

# Using the function coded above, we find the speed-limited v^i and u^0 on both 
# left and right faces.
Valenciav_rU,u4upperZero_r = compute_u0_noif(gamma_faceDD,alpha_face,Valenciav_rU)
Valenciav_lU,u4upperZero_l = compute_u0_noif(gamma_faceDD,alpha_face,Valenciav_lU)

# After importing the relevant module, we'll run the function contained therein 
# with the interpolated metric face values and the reconstructed velocities
# and magnetic fields. Then, we'll simply read in the expressions that that 
# function set, renaming them as convenient and necessary. In particular, we 
# add _r and _l tags and replace u0 with u4upperZero, which is necessary to avoid
# NRPy+ interpreting u0 as something it's not.
# LEFT
u0b.compute_u0_smallb_Poynting__Cartesian(gamma_faceDD,beta_faceU,alpha_face,Valenciav_lU,B_lU)
u_lD = ixp.zerorank1()
u_lU = ixp.zerorank1()

for i in range(DIM):
    u_lD[i] = u0b.uD[i].subs(u0b.u0,u4upperZero_l)
    u_lU[i] = u0b.uU[i].subs(u0b.u0,u4upperZero_l)

smallb4_lU = ixp.zerorank1(DIM=4)
smallb4_lD = ixp.zerorank1(DIM=4)
for mu in range(4):
    smallb4_lU[mu] = u0b.smallb4U[mu].subs(u0b.u0,u4upperZero_l)
    smallb4_lD[mu] = u0b.smallb4D[mu].subs(u0b.u0,u4upperZero_l)

smallb2_l = u0b.smallb2etk.subs(u0b.u0,u4upperZero_l)

# We must repeat this process on both faces.
# RIGHT
u0b.compute_u0_smallb_Poynting__Cartesian(gamma_faceDD,beta_faceU,alpha_face,Valenciav_rU,B_rU)
u_rD = ixp.zerorank1()
u_rU = ixp.zerorank1()

for i in range(DIM):
    u_rD[i] = u0b.uD[i].subs(u0b.u0,u4upperZero_r)
    u_rU[i] = u0b.uU[i].subs(u0b.u0,u4upperZero_r)

smallb4_rU = ixp.zerorank1(DIM=4)
smallb4_rD = ixp.zerorank1(DIM=4)
for mu in range(4):
    smallb4_rU[mu] = u0b.smallb4U[mu].subs(u0b.u0,u4upperZero_r)
    smallb4_rD[mu] = u0b.smallb4D[mu].subs(u0b.u0,u4upperZero_r)

smallb2_r = u0b.smallb2etk.subs(u0b.u0,u4upperZero_r)

```

<a id='hydro_speed'></a>

## Step 2.d: GRFFE characteristic wave speeds \[Back to [top](#toc)\]
$$\label{hydro_speed}$$

Next, we will find the speeds at which the hydrodynamics waves propagate. We start from the speed of light (since FFE deals with very diffuse plasmas), which is $c=1.0$ in our chosen units. We then find the speeds $c_+$ and $c_-$ on each face with the function `find_cp_cm`; then, we find minimum and maximum speeds possible from among those.



Below is the source code for `find_cp_cm`, edited to work with the NRPy+ version of GiRaFFE. One edit we need to make in particular is to the term `psim4*gupii` in the definition of `c`; that was written assuming the use of the conformal metric $\tilde{g}^{ii}$. Since we are not using that here, and are instead using the ADM metric, we should not multiply by $\psi^{-4}$.

```c
static inline void find_cp_cm(REAL &cplus,REAL &cminus,const REAL v02,const REAL u0,
                              const REAL vi,const REAL lapse,const REAL shifti,
                              const REAL gammadet,const REAL gupii) {
  const REAL u0_SQUARED=u0*u0;
  const REAL ONE_OVER_LAPSE_SQUARED = 1.0/(lapse*lapse);
  // sqrtgamma = psi6 -> psim4 = gammadet^(-1.0/3.0)
  const REAL psim4 = pow(gammadet,-1.0/3.0);
  //Find cplus, cminus:
  const REAL a = u0_SQUARED * (1.0-v02) + v02*ONE_OVER_LAPSE_SQUARED;
  const REAL b = 2.0* ( shifti*ONE_OVER_LAPSE_SQUARED * v02 - u0_SQUARED * vi * (1.0-v02) );
  const REAL c = u0_SQUARED*vi*vi * (1.0-v02) - v02 * ( gupii -
                                                               shifti*shifti*ONE_OVER_LAPSE_SQUARED);
  REAL detm = b*b - 4.0*a*c;
  //ORIGINAL LINE OF CODE:
  //if(detm < 0.0) detm = 0.0;
  //New line of code (without the if() statement) has the same effect:
  detm = sqrt(0.5*(detm + fabs(detm))); /* Based on very nice suggestion from Roland Haas */
  
  cplus = 0.5*(detm-b)/a;
  cminus = -0.5*(detm+b)/a;
  if (cplus < cminus) {
    const REAL cp = cminus;
    cminus = cplus;
    cplus = cp;
  }
}
```
Comments documenting this have been excised for brevity, but are reproduced in $\LaTeX$ [below](#derive_speed).

We could use this code directly, but there's substantial improvement we can make by changing the code into a NRPyfied form. Note the `if` statement; NRPy+ does not know how to handle them, so we must eliminate it if we want to leverage NRPy+'s full power. (Calls to `fabs()` are also cheaper than `if` statements.) This can be done if we rewrite this, taking inspiration from the other eliminated `if` statement documented in the above code block:
```c
  cp = 0.5*(detm-b)/a;
  cm = -0.5*(detm+b)/a;
  cplus  = 0.5*(cp+cm+fabs(cp-cm));
  cminus = 0.5*(cp+cm-fabs(cp-cm));
```
This can be simplified further, by substituting `cp` and `cm` into the below equations and eliminating terms as appropriate. First note that `cp+cm = -b/a` and that `cp-cm = detm/a`. Thus,
```c
  cplus  = 0.5*(-b/a + fabs(detm/a));
  cminus = 0.5*(-b/a - fabs(detm/a));
```
This fulfills the original purpose of the `if` statement in the original code because we have guaranteed that $c_+ \geq c_-$.

This leaves us with an expression that can be much more easily NRPyfied. So, we will rewrite the following in NRPy+, making only minimal changes to be proper Python. However, it turns out that we can make this even simpler. In GRFFE, $v_0^2$ is guaranteed to be exactly one. In GRMHD, this speed was calculated as $$v_{0}^{2} = v_{\rm A}^{2} + c_{\rm s}^{2}\left(1-v_{\rm A}^{2}\right),$$ where the Alfv&eacute;n speed $v_{\rm A}^{2}$ $$v_{\rm A}^{2} = \frac{b^{2}}{\rho_{b}h + b^{2}}.$$ So, we can see that when the density $\rho_b$ goes to zero, $v_{0}^{2} = v_{\rm A}^{2} = 1$. Then 
\begin{align}
a &= (u^0)^2 (1-v_0^2) + v_0^2/\alpha^2 \\
&= 1/\alpha^2 \\
b &= 2 \left(\beta^i v_0^2 / \alpha^2 - (u^0)^2 v^i (1-v_0^2)\right) \\
&= 2 \beta^i / \alpha^2 \\
c &= (u^0)^2 (v^i)^2 (1-v_0^2) - v_0^2 \left(g^{ii} - (\beta^i)^2/\alpha^2\right) \\
&= g^{ii} - (\beta^i)^2/\alpha^2,
\end{align}
are simplifications that should save us some time; we can see that $a \geq 0$ is guaranteed. Note that we also force `detm` to be positive. Thus, `detm/a` is guaranteed to be positive itself, rendering the calls to `nrpyAbs()` superfluous. Furthermore, we eliminate any dependence on the Valencia 3-velocity and the time compoenent of the four-velocity, $u^0$. This leaves us free to solve the quadratic in the familiar way: $$c_\pm = \frac{-b \pm \sqrt{b^2-4ac}}{2a}$$.


```python
# We'll write this as a function so that we can calculate the expressions on-demand for any choice of i
def find_cp_cm(lapse,shifti,gupii):
    # Inputs:  u0,vi,lapse,shift,gammadet,gupii
    # Outputs: cplus,cminus 
    
    # a = 1/(alpha^2)
    a = 1/(lapse*lapse)
    # b = 2 beta^i / alpha^2
    b = 2 * shifti /(lapse*lapse)
    # c = -g^{ii} + (beta^i)^2 / alpha^2
    c = - gupii + shifti*shifti/(lapse*lapse)
    
    # Now, we are free to solve the quadratic equation as usual. We take care to avoid passing a
    # negative value to the sqrt function.
    detm = b*b - 4*a*c
    detm = sp.sqrt(sp.Rational(1,2)*(detm + nrpyAbs(detm)))
    cplus  = sp.Rational(1,2)*(-b/a + detm/a)
    cminus = sp.Rational(1,2)*(-b/a - detm/a)
    return cplus,cminus

```

We will now write a function in NRPy+ similar to the one used in the old `GiRaFFE`, allowing us to generate the expressions with less need to copy-and-paste code; the key difference is that this one will be in Python, and generate optimized C code integrated into the rest of the operations. Notice that since we eliminated the dependence on velocities, none of the input quantities are different on either side of the face. So, this function won't really do much besides guarantee that `cmax` and `cmin` are positive and negative, respectively, but we'll leave the machinery here since it is likely to be a useful guide to somebody who wants to something similar. We use the same technique as above to replace the `if` statements inherent to the `MAX()` and `MIN()` functions.


```python
# We'll write this as a function, and call it within HLLE_solver, below.
def find_cmax_cmin(flux_dirn,gamma_faceUU,beta_faceU,alpha_face,gammadet_face):
    # Inputs:  flux direction flux_dirn, Inverse metric gamma_faceUU, shift beta_faceU,
    #          lapse alpha_face, metric determinant gammadet_face
    # Outputs: maximum and minimum characteristic speeds cmax and cmin
    # First, we need to find the characteristic speeds on each face
    cpr,cmr = find_cp_cm(alpha_face,beta_faceU[flux_dirn],gamma_faceUU[flux_dirn][flux_dirn])
    cpl,cml = find_cp_cm(alpha_face,beta_faceU[flux_dirn],gamma_faceUU[flux_dirn][flux_dirn])
    
    # The following algorithms have been verified with random floats:
    
    # Now, we need to set cmax to the larger of cpr,cpl, and 0
    cmax = sp.Rational(1,2)*(cpr+cpl+nrpyAbs(cpr-cpl))
    cmax = sp.Rational(1,2)*(cmax+nrpyAbs(cmax))
    
    # And then, set cmin to the smaller of cmr,cml, and 0
    cmin =  sp.Rational(1,2)*(cmr+cml-nrpyAbs(cmr-cml))
    cmin = -sp.Rational(1,2)*(cmin-nrpyAbs(cmin))
    return cmax,cmin

```

<a id='fluxes'></a>

## Step 2.e: Compute the HLLE fluxes \[Back to [top](#toc)\]
$$\label{fluxes}$$

Finally, we can compute the flux in each direction. This momentum flux in the $m$ direction is defined as $\alpha \sqrt{\gamma} T^m_{\ \ j}$, based on the input `flux_dirn`. We have already defined $\alpha \sqrt{\gamma}$, so all we need to do is calculate $T^m_{\ \ j}$, where $T^{\mu \nu}_{\rm EM} = b^2 u^\mu u^\nu + \frac{1}{2} b^2 g^{\mu \nu} - b^\mu b^\nu$. In doing this index-lowering operation, recall that $g^{\mu \nu} g_{\nu \alpha} = \delta^\mu_\alpha$. We will do so in accordance with the method published by [Harten, Lax, and von Leer](https://epubs.siam.org/doi/pdf/10.1137/1025002)  and [Einfeldt](https://epubs.siam.org/doi/10.1137/0725021) (hereafter HLLE) to solve the Riemann problem. So, we define $f(u) = T^m_{\ \ j}$ on each face as 
$$
f = \alpha \sqrt{\gamma} \left( (\rho+b^2)(u^0 v^m) u_j + (P+\frac{1}{2}b^2) \delta^m_j - b^m b_j \right);
$$
Because $\rho = P = 0$ in GRFFE and $u^0 v^m = u^m$ in general (since $v^m$ is the drift velocity here), this simplifies to 
$$
f = \alpha \sqrt{\gamma} \left( b^2 u^m u_j + \frac{1}{2}b^2 \delta^m_j - b^m b_j \right).
$$
We use $j$ to correspond to the component of the flux we are calculating; that is, $j=1$ corresponds to $x$, and so forth (however, remember that in a NRPy+ 3-vector, the numbers will still run from 0 to 2). $\delta^i_j$ is the standard Kronecker delta. We also define `st_x_r` and `st_x_l` (called $u_{\rm R}$ and $u_{\rm L}$ in HLL; we'll not do that, since we're using $u$ for four-velocity):
$$
{\rm st\_j} = \alpha \sqrt{\gamma} \left( (\rho+b^2) u^0 u_j - b^0 b_j \right),
$$
which, in GRFFE, simplifies to 
$$
{\rm st\_j} = \alpha \sqrt{\gamma} \left( b^2 u^0 u_j - b^0 b_j \right),
$$

and combine based on eq. 3.15 in the HLLE paper,
$$
{\rm st\_j\_flux} = \frac{c_{\rm min} f_{\rm R} + c_{\rm max} f_{\rm L} - c_{\rm min} c_{\rm max} ({\rm st\_x\_r}-{\rm st\_x\_l})}{c_{\rm min} + c_{\rm max}}
$$

We'll write the HLLE step as a function so that we can loop over `flux_dirn` and `mom_comp` and write each version needed as we need it.


```python
# Note the Kronecker delta symbol in the above; let's create it real quick:
# We'll zero-offset everything here, unlike in the original GiRaFFE.
kronecker_delta = ixp.zerorank2()
for i in range(DIM): 
    kronecker_delta[i][i] = 1

# We'll rewrite this assuming that we've passed the entire reconstructed
# gridfunctions. You could also do this with only one point, but then you'd 
# need to declare everything as a Cparam in NRPy+

def HLLE_solver(flux_dirn, mom_comp): 
    # This solves the Riemann problem for the mom_comp component of the momentum
    # flux StildeD in the flux_dirn direction.
    
    # There's a caveat here to consider with the indices here: remember the general rule
    # that a latin index must be incremented if it's in a 4-vector. 
    # f = alpha sqrt{gamma} ( b^2 u^m u_j + (1/2) b^2 delta^m_j - b^m b_j )
    Fr = alpha_sqrt_gamma * (smallb2_r * u_rU[flux_dirn] * u_rD[mom_comp] 
                             + smallb2_r * kronecker_delta[flux_dirn][mom_comp] / 2.0
                             - smallb4_rU[flux_dirn+1] * smallb4_rD[mom_comp+1])
    
    # f = alpha sqrt{gamma} ( b^2 u^m u_j + (1/2) b^2 delta^m_j - b^m b_j )
    Fl = alpha_sqrt_gamma * (smallb2_l * u_lU[flux_dirn] * u_lD[mom_comp] 
                             + smallb2_l * kronecker_delta[flux_dirn][mom_comp] / 2.0
                             - smallb4_lU[flux_dirn+1] * smallb4_lD[mom_comp+1])
    
    # st_j = alpha sqrt{gamma} ( b^2 u^0 u_j - b^0 b_j )
    st_j_r = alpha_sqrt_gamma * (smallb2_r*u4upperZero_r*u_rD[mom_comp] - smallb4_rU[0]*smallb4_rD[mom_comp+1])
    st_j_l = alpha_sqrt_gamma * (smallb2_l*u4upperZero_l*u_lD[mom_comp] - smallb4_lU[0]*smallb4_lD[mom_comp+1])
    
    cmaxL,cminL = find_cmax_cmin(flux_dirn,gamma_faceUU,beta_faceU,alpha_face,gammadet_face)
    
    # st_j_flux = (c_\min f_R + c_\max f_L - c_\min c_\max ( st_j_r - st_j_l )) / (c_\min + c_\max)
    return (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(st_j_r-st_j_l) )/(cmaxL + cminL)

```

Finally, now that we have the equations written, we can write out the files. After writing a brief post string (to close the function), we will loop over `flux_dirn`, creating one C file for each flux direction. In each file, we will include the math to calculate each momentum-flux component `mom_comp` in that direction.

We have written the function `HLLE_solver()` so that we can easily compute the flux as specified by those two indices.


```python
def calculate_Stilde_flux(flux_dirn,Stilde_fluxD):
    for mom_comp in range(DIM):
        Stilde_fluxD[mom_comp] = HLLE_solver(flux_dirn,mom_comp)
    return Stilde_fluxD

```

<a id='code_validation'></a>

# Step 3:  Code Validation against `GiRaFFE_HO.GiRaFFE_Higher_Order_v2` NRPy+ Module \[Back to [top](#toc)\]
$$\label{code_validation}$$


Here, as a code validation check, we verify agreement in the SymPy expressions for the $\texttt{GiRaFFE}$ evolution equations and auxiliary quantities we intend to use between
1. this tutorial and 
2. the NRPy+ [GiRaFFE_HO.GiRaFFE_NRPy_Ccode_library_Stilde_flux](../edit/GiRaFFE_HO/GiRaFFE_NRPy_Ccode_library_Stilde_flux.py) module.



```python
# Reset the list of gridfunctions, as registering a gridfunction
#   twice will spawn an error.
gri.glb_gridfcs_list = []

import GiRaFFE_HO.Stilde_flux as stflux

print("Consistency check between GiRaFFE_Higher_Order tutorial and NRPy+ module: ALL SHOULD BE ZERO.")

for i in range(DIM):
    print("Checking the flux in direction "+str(i))
    Stilde_fluxD = calculate_Stilde_flux(i,Stilde_fluxD)
    gri.glb_gridfcs_list = []
    stflux.calculate_Stilde_flux(i,False)
    for j in range(DIM):
        print("Stilde_fluxD["+str(j)+"] - stflux.Stilde_fluxD["+str(j)+"] = "\
              + str(Stilde_fluxD[j] - stflux.Stilde_fluxD[j]))

# Stilde_fluxD = calculate_Stilde_flux(i,Stilde_fluxD)
# print("u4upperZero_r - stflux.u4upperZero_r = "\
#               + str(u4upperZero_r - stflux.u4upperZero_r))

# for j in range(DIM):
#     print("Valenciav_rU["+str(j)+"] - stflux.Valenciav_rU["+str(j)+"] = "\
#           + str(Valenciav_rU[j] - stflux.Valenciav_rU[j]))
```

    Consistency check between GiRaFFE_Higher_Order tutorial and NRPy+ module: ALL SHOULD BE ZERO.
    Checking the flux in direction 0
    Stilde_fluxD[0] - stflux.Stilde_fluxD[0] = 0
    Stilde_fluxD[1] - stflux.Stilde_fluxD[1] = 0
    Stilde_fluxD[2] - stflux.Stilde_fluxD[2] = 0
    Checking the flux in direction 1
    Stilde_fluxD[0] - stflux.Stilde_fluxD[0] = 0
    Stilde_fluxD[1] - stflux.Stilde_fluxD[1] = 0
    Stilde_fluxD[2] - stflux.Stilde_fluxD[2] = 0
    Checking the flux in direction 2
    Stilde_fluxD[0] - stflux.Stilde_fluxD[0] = 0
    Stilde_fluxD[1] - stflux.Stilde_fluxD[1] = 0
    Stilde_fluxD[2] - stflux.Stilde_fluxD[2] = 0


<a id='derive_speed'></a>

# Step 4: Complete Derivation of the Wave Speeds \[Back to [top](#toc)\]
$$\label{derive_speed}$$

This computes phase speeds in the direction given by flux_dirn. Note that we replace the full dispersion relation with a simpler one, which overestimates the maximum speeds by a factor of ~2. See full discussion around Eqs. 49 and 50 in [Duez, et al.](http://arxiv.org/pdf/astro-ph/0503420.pdf). In summary, we solve the dispersion relation (in, e.g., the $x$-direction) with a wave vector of $k_\mu = (-\omega,k_x,0,0)$. So, we solve the approximate dispersion relation $\omega_{\rm cm}^2 = [v_A^2 + c_s^2 (1-v_A^2)]k_{\rm cm}^2$ for the wave speed $\omega/k_x$, where the sound speed $c_s = \sqrt{\Gamma P/(h \rho_0)}$, the Alfv&eacute;n speed $v_A = 1$ (in GRFFE), $\omega_{\rm cm} = -k_\mu k^\mu$ is the frequency in the comoving frame, $k_{\rm cm}^2 = K_\mu K^\mu$ is the wavenumber squared in the comoving frame, and $K_\mu = (g_{\mu\nu} + u_\mu u_\nu)k^\nu$ is the part of the wave vector normal to the four-velocity $u^\mu$. See below for a complete derivation.

What follows is a complete derivation of the quadratic we solve.
\begin{align}
w_{\rm cm} &= (-k_0 u^0 - k_x u^x) \\
k_{\rm cm}^2 &= K_{\mu} K^{\mu}, \\
K_{\mu} K^{\mu} &= (g_{\mu a} + u_{\mu} u_a) k^a  g^{\mu b} [ (g_{c b} + u_c u_b) k^c ] \\
\rightarrow g^{\mu b} (g_{c b} + u_{c} u_{b}) k^c &= (\delta^{\mu}_c + u_c u^{\mu} ) k^c \\
                 &= (g_{\mu a} + u_{\mu} u_a) k^a  (\delta^{\mu}_c + u_c u^{\mu} ) k^c \\
                 &=[(g_{\mu a} + u_{\mu} u_a) \delta^{\mu}_c + (g_{\mu a} + u_{\mu} u_a) u_c u^{\mu} ] k^c k^a \\
                 &=[(g_{c a} + u_c u_a) + (u_c u_a -  u_a u_c] k^c k^a \\
                 &=(g_{c a} + u_c u_a) k^c k^a \\
                 &= k_a k^a + u^c u^a k_c k_a \\
k^a = g^{\mu a} k_{\mu} &= g^{0 a} k_0 + g^{x a} k_x \\
k_a k^a &= k_0 g^{0 0} k_0 + k_x k_0 g^{0 x} + g^{x 0} k_0 k_x + g^{x x} k_x k_x \\
         &= g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2 \\
u^c u^a k_c k_a &= (u^0 k_0 + u^x k_x) (u^0 k_0 + u^x k_x) = (u^0 k_0)^2 + 2 u^x k_x u^0 k_0 + (u^x k_x)^2 \\
(k_0 u0)^2  + 2 k_x u^x k_0 u^0 + (k_x u^x)^2 &= v_0^2 [ (u^0 k_0)^2 + 2 u^x k_x u^0 k_0 + (u^x k_x)^2 + g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2] \\
(1-v_0^2) (u^0 k_0 + u^x k_x)^2 &= v_0^2 (g^{00} (k_0)^2 + 2 g^{x0} k_0 k_x + g^{xx} (k_x)^2) \\
(1-v_0^2) (u^0 k_0/k_x + u^x)^2 &= v_0^2 (g^{00} (k_0/k_x)^2 + 2 g^{x0} k_0/k_x + g^{xx}) \\
(1-v_0^2) (u^0 X + u^x)^2 &= v_0^2 (g^{00} X^2 + 2 g^{x0} X + g^{xx}) \\
(1-v_0^2) ((u^0)^2 X^2 + 2 u^x (u^0) X + (u^x)^2) &= v_0^2 (g^{00} X^2 + 2 g^{x0} X + g^{xx}) \\
&X^2 ( (1-v_0^2) (u^0)^2 - v_0^2 g^{00}) + X (2 u^x u^0 (1-v_0^2) - 2 v_0^2 g^{x0}) + (1-v_0^2) (u^x)^2 - v_0^2 g^{xx} \\
a &= (1-v_0^2) (u^0)^2 - v_0^2 g^{00} = (1-v_0^2) (u^0)^2 + v_0^2/\alpha^2 \leftarrow {\rm VERIFIED} \\
b &= 2 u^x u^0 (1-v_0^2) - 2 v_0^2 \beta^x/\alpha^2 \leftarrow {\rm VERIFIED,\ } X\rightarrow -X, {\rm because\ } X = -w/k_1, {\rm \ and\ we\ are\ solving\ for} -X. \\
c &= (1-v_0^2) (u^x)^2 - v_0^2 (g^{xx}\psi^{-4} - (\beta^x/\alpha)^2) \leftarrow {\rm VERIFIED} \\
v_0^2 &= v_A^2 + c_s^2 (1 - v_A^2) \\
\end{align}



<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-GiRaFFE_HO_Ccode_library-Stilde-flux.pdf](Tutorial-GiRaFFE_HO_Ccode_library-Stilde-flux.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)


```python
!jupyter nbconvert --to latex --template latex_nrpy_style.tplx --log-level='WARN' Tutorial-GiRaFFE_NRPy_Ccode_library-Stilde-flux.ipynb
!pdflatex -interaction=batchmode Tutorial-GiRaFFE_NRPy_Ccode_library-Stilde-flux.tex
!pdflatex -interaction=batchmode Tutorial-GiRaFFE_NRPy_Ccode_library-Stilde-flux.tex
!pdflatex -interaction=batchmode Tutorial-GiRaFFE_NRPy_Ccode_library-Stilde-flux.tex
!rm -f Tut*.out Tut*.aux Tut*.log
```

    [NbConvertApp] Converting notebook Tutorial-GiRaFFE_NRPy_Ccode_library-Stilde-flux.ipynb to latex
    [NbConvertApp] Writing 19720600 bytes to Tutorial-GiRaFFE_NRPy_Ccode_library-Stilde-flux.tex
    This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)
     restricted \write18 enabled.
    entering extended mode
    This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)
     restricted \write18 enabled.
    entering extended mode
    This is pdfTeX, Version 3.14159265-2.6-1.40.18 (TeX Live 2017/Debian) (preloaded format=pdflatex)
     restricted \write18 enabled.
    entering extended mode


d43f1143257b6f20ef8b51cc563e9704beb88178
ipynb
Jupyter Notebook
in_progress/Tutorial-GiRaFFE_NRPy_Ccode_library-Stilde-flux.ipynb
KAClough/nrpytutorial
2cc3b22cb1092bc10890237dd8ee3b6881c36b52
BSD-2-Clause
2019-12-23T05:31:25.000Z
2019-12-23T05:31:25.000Z
in_progress/Tutorial-GiRaFFE_NRPy_Ccode_library-Stilde-flux.ipynb
KAClough/nrpytutorial
2cc3b22cb1092bc10890237dd8ee3b6881c36b52
BSD-2-Clause
in_progress/Tutorial-GiRaFFE_NRPy_Ccode_library-Stilde-flux.ipynb
KAClough/nrpytutorial
2cc3b22cb1092bc10890237dd8ee3b6881c36b52
BSD-2-Clause
Qwen/Qwen-72B
1. YES
2. YES