The Characteristic GRFFE speeds¶

Author: Patrick Nelson¶

This notebook documents the function from the original GiRaFFE that calculates the speeds of the characteristic GRFFE hydrodynamics waves.¶

Notebook Status: Validated

Validation Notes: This code has been validated to round-off level agreement with the corresponding code in the original GiRaFFE by means of its implementation in Tutorial-GiRaFFE_NRPy-Stilde-flux and Tutorial-GiRaFFE_NRPy-Afield_flux_handwritten

NRPy+ Source Code for this module:¶

• GiRaFFE_NRPy/GiRaFFE_NRPy_Characteristic_Speeds.py

Introduction¶

Here, 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. A complete derivation of these speeds is included below as an appendix; here, for brevity, we simply start with the implementation of the code in the original GiRaFFE. That code itself is adapted from the more general IllinoisGRMHD code; we will implement several more simplifications to help improve performance.

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

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.

Table of Contents¶

This notebook is organized as follows

1. Step 1: Preliminaries
2. Step 2: Apply GRFFE assumptions to simplify expressions
3. Step 3: Compute the characteristic speeds
   1. Step 3.a: Compute $c_+$ and $c_-$
   2. Step 3.b: Compute $c_\max$ and $c_\min$
4. Step 4: Code Validation against GiRaFFE_NRPy.GiRaFFE_NRPy_Characteristic_Speeds NRPy+ Module
5. Step 5: Complete Derivation of the Wave Speeds
6. Step 6: Output this notebook to $\LaTeX$-formatted PDF file

Step 1: Preliminaries [Back to top]¶

This first block of code imports the core NRPy+ functionality after first adding the main NRPy+ directory to the path.

Step 2: Apply GRFFE assumptions to simplify expressions [Back to top]¶

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 these, 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:

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,

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é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(\gamma^{ii} - (\beta^i)^2/\alpha^2\right) \\ &= -\gamma^{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}$$.

Step 3: Compute the characteristic speeds [Back to top]¶

Step 3.a: Compute $c_+$ and $c_-$ [Back to top]¶

We now have the expressions we need for $c_+$ and $c_-$: $$c_\pm = \frac{-b \pm \sqrt{b^2-4ac}}{2a},$$ where \begin{align} a &= 1/\alpha^2 \\ b &= 2 \beta^i / \alpha^2 \\ c &= -\gamma^{ii} + (\beta^i)^2/\alpha^2. \end{align}

Step 3.b: Compute $c_\max$ and $c_\min$ [Back to top]¶

In flat spacetime, where $\alpha=1$, $\beta^i=0$, and $\gamma^{ij} = \delta^{ij}$, $c_+ > 0$ and $c_- < 0$. For the HLLE solver, we will need both cmax and cmin to be positive; we also want to choose the speed that is larger in magnitude because overestimating the characteristic speeds will help damp unwanted oscillations. (However, in GRFFE, we only get one $c_+$ and one $c_-$, so we only need to fix the signs here.) Hence, the following function.

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

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

Here, as a code validation check, we verify agreement in the SymPy expressions for the GiRaFFE evolution equations and auxiliary quantities we intend to use between

1. this tutorial and
2. the NRPy+ GiRaFFE_NRPy.GiRaFFE_NRPy_Characteristic_Speeds module.

This first validation directly compares the sympy expressions. This is generally quicker and more reliable, but might overlook some complexities in implementing the C code.

Step 5: Complete Derivation of the Wave Speeds [Back to top]¶

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.. 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é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. We start from the following relations: \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 ] \\ \end{align} The last term of the above can be written as follow: $$ (g_{c b} + u_{c} u_{b}) k^c = (\delta^{\mu}_c + u_c u^{\mu} ) k^c $$

Then, \begin{align} 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 ] \\ &= (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 u^0)^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}) \\ 0 &= 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 (\gamma^{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}

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

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