GiRaFFE_NRPy C code library: Conservative-to-Primitive and Primitive-to-Conservative Solvers¶

Author: Patrick Nelson¶

This notebook writes and documents the C code that GiRaFFE_NRPy uses in order to update the Valencia 3-velocity at each timestep. It also computes corrections to the densitized Poynting flux in order to keep the physical quantities from violating the GRFFE constraints.¶

Notebook Status: Validated

Validation Notes: These functions have been validated to round-off precision against the coresponding functions in the original GiRaFFE.

NRPy+ Source Code for this module:¶

• GiRaFFE_NRPy/GiRaFFE_NRPy_C2P_P2C.py

Introduction:¶

These algorithms are adapted from the original GiRaFFE code (see arxiv:1704.00599v2), based on the description in arXiv:1310.3274v2. They have been fully NRPyfied and modified to use the Valencia 3-velocity instead of the drift velocity.

The algorithm to do this is as follows:

1. Apply fixes to ${\tilde S}_i$
   1. Enforce the orthogonality of ${\tilde S}_i$ and $B^i$
      1. ${\tilde S}_i \rightarrow {\tilde S}_i - ({\tilde S}_j {\tilde B}^j) {\tilde B}_i/{\tilde B}^2$
   2. Rescale ${\tilde S}_i$ to limit the Lorentz factor and enforce the velocity cap
      1. $f = \sqrt{(1-\gamma_{\max}^{-2}){\tilde B}^4/(16 \pi^2 \gamma {\tilde S}^2)}$
      2. ${\tilde S}_i \rightarrow {\tilde S}_i \min(1,f)$
   3. Recompute the velocities at the new timestep
      1. $v^i = 4 \pi \gamma^{ij} {\tilde S}_j \gamma^{-1/2} B^{-2}$
2. Enforce the Current Sheet prescription
   1. Zero the velocity normal to the sheet
      1. ${\tilde n}_i v^i = 0$
   2. Recompute the Poynting flux to be consistent.

Each of these steps can be toggled on/off by changing the following NRPy+ parameters, specified in the python module:

par.initialize_param(par.glb_param(type="bool", module=thismodule, parname="enforce_orthogonality_StildeD_BtildeU", defaultval=True))
par.initialize_param(par.glb_param(type="bool", module=thismodule, parname="enforce_speed_limit_StildeD", defaultval=True))
par.initialize_param(par.glb_param(type="bool", module=thismodule, parname="enforce_current_sheet_prescription", defaultval=True))

Table of Contents¶

This notebook is organized as follows

1. Step 1: The conservative-to-primitive solver
   1. Step 1.a: Enforce the orthogonality of $\tilde{S}_i$ and $B^i$
   2. Step 1.b: Rescale ${\tilde S}_i$ to limit the Lorentz factor and enforce the velocity cap
   3. Step 1.c: Recompute the velocities at the new timestep
   4. Step 1.d: Enforce the Current Sheet prescription
2. Step 2: The primitive-to-conservative solver
3. Step 3: Code Validation against
4. Step 4: Output this notebook to $\LaTeX$-formatted PDF file

Step 1: The conservative-to-primitive solver [Back to top]¶

We start with the Conservative-to-Primitive solver. This function is called after the vector potential and Poynting vector have been evolved at a timestep and updates the velocities. The algorithm will be as follows:

1. Apply fixes to ${\tilde S}_i$
   1. Enforce the orthogonality of ${\tilde S}_i$ and $B^i$
      1. ${\tilde S}_i \rightarrow {\tilde S}_i - ({\tilde S}_j {\tilde B}^j) {\tilde B}_i/{\tilde B}^2$
   2. Rescale ${\tilde S}_i$ to limit the Lorentz factor and enforce the velocity cap
      1. $f = \sqrt{(1-\gamma_{\max}^{-2}){\tilde B}^4/(16 \pi^2 \gamma {\tilde S}^2)}$
      2. ${\tilde S}_i \rightarrow {\tilde S}_i \min(1,f)$
   3. Recompute the velocities at the new timestep
      1. $v^i = 4 \pi \gamma^{ij} {\tilde S}_j \gamma^{-1/2} B^{-2}$
2. Enforce the Current Sheet prescription
   1. ${\tilde n}_i v^i = 0$

Step 1.a: Enforce the orthogonality of $\tilde{S}_i$ and $B^i$ [Back to top]¶

Now, we will enforce the orthogonality of the magnetic field and densitized poynting flux using Eq. 22 of arxiv:1704.00599v2: $${\tilde S}_i \rightarrow {\tilde S}_i - ({\tilde S}_j {\tilde B}^j) {\tilde B}_i/{\tilde B}^2$$ First, we compute the inner products ${\tilde S}_j {\tilde B}^j$ and ${\tilde B}^2 = \gamma_{ij} {\tilde B}^i {\tilde B}^j,$ where $\tilde{B}^i = B^i \sqrt{\gamma}$. Then, we subtract $({\tilde S}_j {\tilde B}^j) {\tilde B}_i/{\tilde B}^2$ from ${\tilde S}_i$. We thus guarantee that ${\tilde S}_i B^i=0$.

Having fixed ${\tilde S}_i$, we will also compute the related quantities ${\tilde S}^i = \gamma^{ij} {\tilde S}_j$ and ${\tilde S}^2 = {\tilde S}_i {\tilde S}^i$.

Step 1.b: Rescale ${\tilde S}_i$ to limit the Lorentz factor and enforce the velocity cap [Back to top]¶

The next fix that we will apply limits the Lorentz factor using Eqs. 92 and 93 of arXiv:1310.3274v2. That is, we define the factor $f$ as $$f = \sqrt{(1-\Gamma_{\max}^{-2}){\tilde B}^4/(16 \pi^2 \gamma {\tilde S}^2)}.$$

If $f<1$, we rescale the components of ${\tilde S}_i$ by $f$. That is, we must set $${\tilde S}_i \rightarrow {\tilde S}_i \min(1,f).$$

Here, we will use a formulation of the min() function that does not use if: \begin{equation} \min(a,b) = \frac{1}{2} \left( a+b - \lvert a-b \rvert \right), \end{equation} or, in code,

min_noif(a,b) = sp.Rational(1,2)*(a+b-nrpyAbs(a-b))

Step 1.c: Recompute the velocities at the new timestep [Back to top]¶

Finally, we can calculate the velocities. In arxiv:1704.00599v2, Eq. 16 gives the drift velocity as $$v^i = 4 \pi \alpha \gamma^{ij} {\tilde S}_j \gamma^{-1/2} B^{-2} - \beta^i.$$ However, we wish to use the Valencia velocity instead. Since the Valencia velocity $\bar{v}^i = \frac{1}{\alpha} \left( v^i + \beta^i \right)$, we will code $$\bar{v}^i = 4 \pi \frac{\gamma^{ij} {\tilde S}_j}{\sqrt{\gamma} B^2}.$$

Step 1.d: Enforce the Current Sheet prescription [Back to top]¶

Now, we seek to handle any current sheets (a physically important phenomenon) that might form. This algorithm, given as Eq. 23 in arxiv:1704.00599v2, will preserve current sheets that form in the xy-plane by preventing our numerical scheme from dissipating them. After fixing the z-component of the velocity, we recompute the conservative variables $\tilde{S}_i$ to be consistent with the new velocities.

Thus, if we are within four gridpoints of $z=0$, we set the component of the velocity perpendicular to the current sheet to zero by $n_i v^i = 0$, where $n_i = \gamma_{ij} n^j$ is a unit normal to the current sheet and $n^j = \delta^{jz} = (0\ 0\ 1)$. For drift velocity, this means we just set $$ v^z = -\frac{\gamma_{xz} v^x + \gamma_{yz} v^y}{\gamma_{zz}}. $$ This reduces to $v^z = 0$ in flat space, as one would expect. We then express the Valencia velocity in terms of the drift velocity as $\bar{v}^i = \frac{1}{\alpha} \left( v^i + \beta^i \right)$.

Below is an experiment in coding this more abstractly. While it works, it's a bit harder to follow than the direct approach, which is what is coded above

# Set the Cartesian normal vector. This can be expanded later to arbitrary sheets and coordinate systems.
nU = ixp.zerorank1()
nU[2] = 1
# Lower the index, as usual:
nD = ixp.zerorank1()
for i in range(3):
    for j in range(3):
        nD[i] = gammaDD[i][j]*nU[j]

if par.parval_from_str("enforce_current_sheet_prescription"):
    # Calculate the drift velocity
    driftvU = ixp.declarerank1("driftvU")

    inner_product = sp.sympify(0)
    for i in range(3):
        inner_product += driftvU[i]*nD[i] # This is the portion of the drift velocity normal to the z plane
                                          # In flat space, this is just v^z
    # We'll use a sympy utility to solve for v^z. This should make it easier to generalize later
    newdriftvU2 = sp.solve(inner_product,driftvU[2]) # This outputs a list with a single element. 
                                                     # Take the 0th element so .subs() works right.
    newdriftvU2 = newdriftvU2[0] # In flat space this reduces to v^z=0
    for i in range(3):
        # Now, we substitute drift velocity in terms of our preferred Valencia velocity
        newdriftvU2 = newdriftvU2.subs(driftvU[i],alpha*ValenciavU[i]-betaU[i])

    # Now that we have the z component, it's time to substitute its Valencia form in.
    # Remember, we only do this if abs(z) < (k+0.01)*dz. Note that we add 0.01; this helps
    # avoid floating point errors and division by zero. This is the same as abs(z) - (k+0.01)*dz<0
    boundary = nrpyAbs(rfm.xx[2]) - (grid_points_from_z_plane+sp.Rational(1,100))*gri.dxx[2]
    ValenciavU[2] = min_normal0(boundary)*(newdriftvU2+betaU[2])/alpha \
                  + max_normal0(boundary)*ValenciavU[2]

Step 2: The primitive-to-conservative solver [Back to top]¶

This function is used to recompute the conservatives $\tilde{S}_i$ after the 3-velocity is changed as part of the current sheet prescription using Eq. 21 of arxiv:1704.00599v2. It implements the same equation used to compute the initial Poynting flux from the initial velocity: $$\tilde{S}_i = \gamma_{ij} \frac{\bar{v}^j \sqrt{\gamma}B^2}{4 \pi}$$ in terms of the Valencia 3-velocity. In the implementation here, we first calculate $B^2 = \gamma_{ij} B^i B^j$, then $v_i = \gamma_{ij} v^j$ before we calculate the equivalent expression $$\tilde{S}_i = \frac{\bar{v}_i \sqrt{\gamma}B^2}{4 \pi}.$$

Here, we will simply let the NRPy+ GRFFE module handle this part, since that is already validated.

Step 3: Code Validation [Back to top]¶

As a code validation check, we will verify agreement in the SymPy expressions between

1. this tutorial and
2. the NRPy+ GiRaFFE_NRPy/GiRaFFE_NRPy_C2P_P2C.py module.

We will test the C2P algorithms first.

Now, we will reset the variables and test the P2C function.

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