GiRaFFE_NRPy
: Solving the Induction Equation¶This notebook documents the function from the original GiRaFFE
that calculates the flux for $A_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), modified for the case of GRFFE.
Notebook Status: Validated
Validation Notes: This code has been validated by showing that it converges to the exact answer at the expected order
Our goal in this module is to write the code necessary to solve the induction equation $$ \partial_t A_i = \underbrace{\epsilon_{ijk} v^j B^k}_{\rm Flux\ terms} - \underbrace{\partial_i \left(\alpha \Phi - \beta^j A_j \right)}_{\rm Gauge\ terms}. $$ To properly handle the flux terms and avoiding problems with shocks, we cannot simply take a cross product of the velocity and magnetic field at the cell centers. Instead, we must solve the Riemann problem at the cell faces using the reconstructed values of the velocity and magnetic field on either side of the cell faces. The reconstruction is done using PPM (see here); in this module, we will assume that that step has already been done. Metric quantities are assumed to have been interpolated to cell faces, as is done in this tutorial.
Tóth's paper, Eqs. 30 and 31, are one of the first implementations of such a scheme. The original GiRaFFE used a 2D version of the algorithm from Del Zanna, et al. (2002); but since we are not using staggered grids, we can greatly simplify this algorithm with respect to the version used in the original GiRaFFE
. Instead, we will adapt the implementations of the algorithm used in Mewes, et al. (2020) and Giacomazzo, et al. (2011), Eqs. 3-11.
We first write the flux contribution to the induction equation RHS as $$ \partial_t A_i = -E_i, $$ where the electric field $E_i$ is given in ideal MHD (of which FFE is a subset) as $$ -E_i = \epsilon_{ijk} v^j B^k, $$ where $v^i$ is the drift velocity, $B^i$ is the magnetic field, and $\epsilon_{ijk} = \sqrt{\gamma} [ijk]$ is the Levi-Civita tensor. In Cartesian coordinates, \begin{align} -E_x &= [F^y(B^z)]_x = -[F^z(B^y)]_x \\ -E_y &= [F^z(B^x)]_y = -[F^x(B^z)]_y \\ -E_z &= [F^x(B^y)]_z = -[F^y(B^x)]_z, \\ \end{align} where $$ [F^i(B^j)]_k = \sqrt{\gamma} (v^i B^j - v^j B^i). $$ To compute the actual contribution to the RHS in some direction $i$, we average the above listed field as calculated on the $+j$, $-j$, $+k$, and $-k$ faces. That is, at some point $(i,j,k)$ on the grid, \begin{align} -E_x(x_i,y_j,z_k) &= \frac{1}{4} \left( [F_{\rm HLL}^y(B^z)]_{x(i,j+1/2,k)}+[F_{\rm HLL}^y(B^z)]_{x(i,j-1/2,k)}-[F_{\rm HLL}^z(B^y)]_{x(i,j,k+1/2)}-[F_{\rm HLL}^z(B^y)]_{x(i,j,k-1/2)} \right) \\ -E_y(x_i,y_j,z_k) &= \frac{1}{4} \left( [F_{\rm HLL}^z(B^x)]_{y(i,j,k+1/2)}+[F_{\rm HLL}^z(B^x)]_{y(i,j,k-1/2)}-[F_{\rm HLL}^x(B^z)]_{y(i+1/2,j,k)}-[F_{\rm HLL}^x(B^z)]_{y(i-1/2,j,k)} \right) \\ -E_z(x_i,y_j,z_k) &= \frac{1}{4} \left( [F_{\rm HLL}^x(B^y)]_{z(i+1/2,j,k)}+[F_{\rm HLL}^x(B^y)]_{z(i-1/2,j,k)}-[F_{\rm HLL}^y(B^x)]_{z(i,j+1/2,k)}-[F_{\rm HLL}^y(B^x)]_{z(i,j-1/2,k)} \right). \\ \end{align} Note the use of $F_{\rm HLL}$ here. This change signifies that the quantity output here is from the HLLE Riemann solver. Note also the indices on the fluxes. Values of $\pm 1/2$ indicate that these are computed on cell faces using the reconstructed values of $v^i$ and $B^i$ and the interpolated values of the metric gridfunctions. So, $$ F_{\rm HLL}^i(B^j) = \frac{c_{\rm min} F_{\rm R}^i(B^j) + c_{\rm max} F_{\rm L}^i(B^j) - c_{\rm min} c_{\rm max} (B_{\rm R}^j-B_{\rm L}^j)}{c_{\rm min} + c_{\rm max}}. $$
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 = \left. \left(-b \pm \sqrt{b^2-4ac}\right)\middle/ \left(2a\right) \right.$ 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.$$ An outline of a general finite-volume method is as follows, with the current step in bold:
We will assume in this notebook that the reconstructed velocities and magnetic fields are available on cell faces as input. We will also assume that the metric gridfunctions have been interpolated on the metric faces.
Solving the Riemann problem, then, consists of two substeps: First, we compute the flux through each face of the cell. Then, we add the average of these fluxes to the right-hand side of the evolution equation for the vector potential.
This notebook is organized as follows
# Step 0: Add NRPy's directory to the path
# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory
import os, sys # Standard Python modules for multiplatform OS-level functions
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
sys.path.append(nrpy_dir_path)
from outputC import outCfunction, outputC # NRPy+: Core C code output module
import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends
import NRPy_param_funcs as par # NRPy+: Parameter interface
import grid as gri # NRPy+: Functions having to do with numerical grids
import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
thismodule = "GiRaFFE_NRPy-Afield_flux"
import GRHD.equations as GRHD
# import GRFFE.equations as GRFFE
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}$.
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.
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}$$.
# 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,gammaUUii):
# Inputs: u0,vi,lapse,shift,gammadet,gupii
# Outputs: cplus,cminus
# a = 1/(alpha^2)
a = sp.sympify(1)/(lapse*lapse)
# b = 2 beta^i / alpha^2
b = sp.sympify(2) * shifti /(lapse*lapse)
# c = -g^{ii} + (beta^i)^2 / alpha^2
c = - gammaUUii + 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 - sp.sympify(4)*a*c
import Min_Max_and_Piecewise_Expressions as noif
detm = sp.sqrt(noif.max_noif(sp.sympify(0),detm))
global cplus,cminus
cplus = sp.Rational(1,2)*(-b/a + detm/a)
cminus = sp.Rational(1,2)*(-b/a - detm/a)
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. The only modifications we'll make are those necessary to eliminate calls to fabs(0)
in the C code. We use the same technique as above to replace the if
statements inherent to the MAX()
and MIN()
functions.
# We'll write this as a function, and call it within HLLE_solver, below.
def find_cmax_cmin(field_comp,gamma_faceDD,beta_faceU,alpha_face):
# Inputs: flux direction field_comp, 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
gamma_faceUU,unusedgammaDET = ixp.generic_matrix_inverter3x3(gamma_faceDD)
# Original needed for GRMHD
# find_cp_cm(alpha_face,beta_faceU[field_comp],gamma_faceUU[field_comp][field_comp])
# cpr = cplus
# cmr = cminus
# find_cp_cm(alpha_face,beta_faceU[field_comp],gamma_faceUU[field_comp][field_comp])
# cpl = cplus
# cml = cminus
find_cp_cm(alpha_face,beta_faceU[field_comp],gamma_faceUU[field_comp][field_comp])
cp = cplus
cm = cminus
# The following algorithms have been verified with random floats:
global cmax,cmin
# Now, we need to set cmax to the larger of cpr,cpl, and 0
import Min_Max_and_Piecewise_Expressions as noif
cmax = noif.max_noif(cp,sp.sympify(0))
# And then, set cmin to the smaller of cmr,cml, and 0
cmin = -noif.min_noif(cm,sp.sympify(0))
Here, we we calculate the flux and state vectors for the electric field. The flux vector is here given as $$ [F^i(B^j)]_k = \sqrt{\gamma} (v^i B^j - v^j B^i). $$ Here, $v^i$ is the drift velocity and $B^i$ is the magnetic field. This can be easily handled for an input flux direction $i$ with $$ [F^j(B^k)]_i = \epsilon_{ijk} v^j B^k, $$ where $\epsilon_{ijk} = \sqrt{\gamma} [ijk]$ and $[ijk]$ is the Levi-Civita symbol.
The state vector is simply the magnetic field $B^j$.
def calculate_flux_and_state_for_Induction(field_comp,flux_dirn, gammaDD,betaU,alpha,ValenciavU,BU):
# Define Levi-Civita symbol
def define_LeviCivitaSymbol_rank3(DIM=-1):
if DIM == -1:
DIM = par.parval_from_str("DIM")
LeviCivitaSymbol = ixp.zerorank3()
for i in range(DIM):
for j in range(DIM):
for k in range(DIM):
# From https://codegolf.stackexchange.com/questions/160359/levi-civita-symbol :
LeviCivitaSymbol[i][j][k] = (i - j) * (j - k) * (k - i) * sp.Rational(1,2)
return LeviCivitaSymbol
GRHD.compute_sqrtgammaDET(gammaDD)
# Here, we import the Levi-Civita tensor and compute the tensor with lower indices
LeviCivitaDDD = define_LeviCivitaSymbol_rank3()
for i in range(3):
for j in range(3):
for k in range(3):
LeviCivitaDDD[i][j][k] *= GRHD.sqrtgammaDET
global U,F
# Flux F = \epsilon_{ijk} v^j B^k
F = sp.sympify(0)
for j in range(3):
for k in range(3):
F += LeviCivitaDDD[field_comp][j][k] * (alpha*ValenciavU[j]-betaU[j]) * BU[k]
# U = B^i
U = BU[flux_dirn]
Now, we write a standard HLLE solver based on eq. 3.15 in the HLLE paper, $$ F^{\rm HLL} = \frac{c_{\rm min} F_{\rm R} + c_{\rm max} F_{\rm L} - c_{\rm min} c_{\rm max} (U_{\rm R}-U_{\rm L})}{c_{\rm min} + c_{\rm max}} $$
def HLLE_solver(cmax, cmin, Fr, Fl, Ur, Ul):
# This solves the Riemann problem for the flux of E_i in one direction
# F^HLL = (c_\min f_R + c_\max f_L - c_\min c_\max ( st_j_r - st_j_l )) / (c_\min + c_\max)
return (cmin*Fr + cmax*Fl - cmin*cmax*(Ur-Ul) )/(cmax + cmin)
Here, we will use the function we just wrote to calculate the flux through a face. We will pass the reconstructed Valencia 3-velocity and magnetic field on either side of an interface to this function (designated as the "left" and "right" sides) along with the value of the 3-metric, shift vector, and lapse function on the interface. The parameter flux_dirn
specifies which face through which we are calculating the flux. However, unlike when we used this method to calculate the flux term, the RHS of each component of $A_i$ does not depend on all three of the flux directions. Instead, the flux of one component of the $E_i$ field depends on flux through the faces in the other two directions. This will be handled when we generate the C function, as demonstrated in the example code after this next function.
Note that we allow the user to declare their own gridfunctions if they wish, and default to declaring basic symbols if they are not provided. The default names are chosen to imply interpolation of the metric gridfunctions and reconstruction of the primitives.
def calculate_E_i_flux(flux_dirn,alpha_face=None,gamma_faceDD=None,beta_faceU=None,\
Valenciav_rU=None,B_rU=None,Valenciav_lU=None,B_lU=None):
global E_fluxD
E_fluxD = ixp.zerorank1()
for field_comp in range(3):
find_cmax_cmin(field_comp,gamma_faceDD,beta_faceU,alpha_face)
calculate_flux_and_state_for_Induction(field_comp,flux_dirn, gamma_faceDD,beta_faceU,alpha_face,\
Valenciav_rU,B_rU)
Fr = F
Ur = U
calculate_flux_and_state_for_Induction(field_comp,flux_dirn, gamma_faceDD,beta_faceU,alpha_face,\
Valenciav_lU,B_lU)
Fl = F
Ul = U
E_fluxD[field_comp] += HLLE_solver(cmax, cmin, Fr, Fl, Ur, Ul)
Below, we will write some example code to use the above functions to generate C code for GiRaFFE_NRPy
. We need to write our own memory reads and writes because we need to add contributions from both faces in a given direction, which is expressed in the code as adding contributions from adjacent gridpoints to the RHS, which is not something FD_outputC
can handle. The .replace()
function calls adapt these reads and writes to the different directions. Note that, for reconstructions in a given direction, the fluxes are only added to the other two components, as can be seen in the equations we are implementing.
\begin{align}
-E_x(x_i,y_j,z_k) &= \frac{1}{4} \left( [F_{\rm HLL}^y(B^z)]_{x(i,j+1/2,k)}+[F_{\rm HLL}^y(B^z)]_{x(i,j-1/2,k)}-[F_{\rm HLL}^z(B^y)]_{x(i,j,k+1/2)}-[F_{\rm HLL}^z(B^y)]_{x(i,j,k-1/2)} \right) \\
-E_y(x_i,y_j,z_k) &= \frac{1}{4} \left( [F_{\rm HLL}^z(B^x)]_{y(i,j,k+1/2)}+[F_{\rm HLL}^z(B^x)]_{y(i,j,k-1/2)}-[F_{\rm HLL}^x(B^z)]_{y(i+1/2,j,k)}-[F_{\rm HLL}^x(B^z)]_{y(i-1/2,j,k)} \right) \\
-E_z(x_i,y_j,z_k) &= \frac{1}{4} \left( [F_{\rm HLL}^x(B^y)]_{z(i+1/2,j,k)}+[F_{\rm HLL}^x(B^y)]_{z(i-1/2,j,k)}-[F_{\rm HLL}^y(B^x)]_{z(i,j+1/2,k)}-[F_{\rm HLL}^y(B^x)]_{z(i,j-1/2,k)} \right). \\
\end{align}
From this, we can see that when, for instance, we reconstruct and interpolate in the $x$-direction, we must add only to the $y$- and $z$-components of the electric field.
Recall that when we reconstructed the velocity and magnetic field, we constructed to the $i-1/2$ face, so the data at $i+1/2$ is stored at $i+1$.
def generate_Afield_flux_function_files(out_dir,subdir,alpha_face,gamma_faceDD,beta_faceU,\
Valenciav_rU,B_rU,Valenciav_lU,B_lU,inputs_provided=True):
if not inputs_provided:
# declare all variables
alpha_face = sp.symbols(alpha_face)
beta_faceU = ixp.declarerank1("beta_faceU")
gamma_faceDD = ixp.declarerank2("gamma_faceDD","sym01")
Valenciav_rU = ixp.declarerank1("Valenciav_rU")
B_rU = ixp.declarerank1("B_rU")
Valenciav_lU = ixp.declarerank1("Valenciav_lU")
B_lU = ixp.declarerank1("B_lU")
Memory_Read = """const double alpha_face = auxevol_gfs[IDX4S(ALPHA_FACEGF, i0,i1,i2)];
const double gamma_faceDD00 = auxevol_gfs[IDX4S(GAMMA_FACEDD00GF, i0,i1,i2)];
const double gamma_faceDD01 = auxevol_gfs[IDX4S(GAMMA_FACEDD01GF, i0,i1,i2)];
const double gamma_faceDD02 = auxevol_gfs[IDX4S(GAMMA_FACEDD02GF, i0,i1,i2)];
const double gamma_faceDD11 = auxevol_gfs[IDX4S(GAMMA_FACEDD11GF, i0,i1,i2)];
const double gamma_faceDD12 = auxevol_gfs[IDX4S(GAMMA_FACEDD12GF, i0,i1,i2)];
const double gamma_faceDD22 = auxevol_gfs[IDX4S(GAMMA_FACEDD22GF, i0,i1,i2)];
const double beta_faceU0 = auxevol_gfs[IDX4S(BETA_FACEU0GF, i0,i1,i2)];
const double beta_faceU1 = auxevol_gfs[IDX4S(BETA_FACEU1GF, i0,i1,i2)];
const double beta_faceU2 = auxevol_gfs[IDX4S(BETA_FACEU2GF, i0,i1,i2)];
const double Valenciav_rU0 = auxevol_gfs[IDX4S(VALENCIAV_RU0GF, i0,i1,i2)];
const double Valenciav_rU1 = auxevol_gfs[IDX4S(VALENCIAV_RU1GF, i0,i1,i2)];
const double Valenciav_rU2 = auxevol_gfs[IDX4S(VALENCIAV_RU2GF, i0,i1,i2)];
const double B_rU0 = auxevol_gfs[IDX4S(B_RU0GF, i0,i1,i2)];
const double B_rU1 = auxevol_gfs[IDX4S(B_RU1GF, i0,i1,i2)];
const double B_rU2 = auxevol_gfs[IDX4S(B_RU2GF, i0,i1,i2)];
const double Valenciav_lU0 = auxevol_gfs[IDX4S(VALENCIAV_LU0GF, i0,i1,i2)];
const double Valenciav_lU1 = auxevol_gfs[IDX4S(VALENCIAV_LU1GF, i0,i1,i2)];
const double Valenciav_lU2 = auxevol_gfs[IDX4S(VALENCIAV_LU2GF, i0,i1,i2)];
const double B_lU0 = auxevol_gfs[IDX4S(B_LU0GF, i0,i1,i2)];
const double B_lU1 = auxevol_gfs[IDX4S(B_LU1GF, i0,i1,i2)];
const double B_lU2 = auxevol_gfs[IDX4S(B_LU2GF, i0,i1,i2)];
REAL A_rhsD0 = 0; REAL A_rhsD1 = 0; REAL A_rhsD2 = 0;
"""
Memory_Write = """rhs_gfs[IDX4S(AD0GF,i0,i1,i2)] += A_rhsD0;
rhs_gfs[IDX4S(AD1GF,i0,i1,i2)] += A_rhsD1;
rhs_gfs[IDX4S(AD2GF,i0,i1,i2)] += A_rhsD2;
"""
indices = ["i0","i1","i2"]
indicesp1 = ["i0+1","i1+1","i2+1"]
for flux_dirn in range(3):
calculate_E_i_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,\
Valenciav_rU,B_rU,Valenciav_lU,B_lU)
E_field_to_print = [\
sp.Rational(1,4)*E_fluxD[(flux_dirn+1)%3],
sp.Rational(1,4)*E_fluxD[(flux_dirn+2)%3],
]
E_field_names = [\
"A_rhsD"+str((flux_dirn+1)%3),
"A_rhsD"+str((flux_dirn+2)%3),
]
desc = "Calculate the electric flux on the left face in direction " + str(flux_dirn) + "."
name = "calculate_E_field_D" + str(flux_dirn) + "_right"
outCfunction(
outfile = os.path.join(out_dir,subdir,name+".h"), desc=desc, name=name,
params ="const paramstruct *params,const REAL *auxevol_gfs,REAL *rhs_gfs",
body = Memory_Read \
+outputC(E_field_to_print,E_field_names,"returnstring",params="outCverbose=False")\
+Memory_Write,
loopopts ="InteriorPoints",
rel_path_to_Cparams=os.path.join("../"))
desc = "Calculate the electric flux on the left face in direction " + str(flux_dirn) + "."
name = "calculate_E_field_D" + str(flux_dirn) + "_left"
outCfunction(
outfile = os.path.join(out_dir,subdir,name+".h"), desc=desc, name=name,
params ="const paramstruct *params,const REAL *auxevol_gfs,REAL *rhs_gfs",
body = Memory_Read.replace(indices[flux_dirn],indicesp1[flux_dirn]) \
+outputC(E_field_to_print,E_field_names,"returnstring",params="outCverbose=False")\
+Memory_Write,
loopopts ="InteriorPoints",
rel_path_to_Cparams=os.path.join("../"))
GiRaFFE_NRPy.Induction_Equation
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
Below are the gridfunction registrations we will need for testing. We will pass these to the above functions to self-validate the module that corresponds with this tutorial.
all_passed=True
def comp_func(expr1,expr2,basename,prefixname2="C2P_P2C."):
if str(expr1-expr2)!="0":
print(basename+" - "+prefixname2+basename+" = "+ str(expr1-expr2))
all_passed=False
def gfnm(basename,idx1,idx2=None,idx3=None):
if idx2 is None:
return basename+"["+str(idx1)+"]"
if idx3 is None:
return basename+"["+str(idx1)+"]["+str(idx2)+"]"
return basename+"["+str(idx1)+"]["+str(idx2)+"]["+str(idx3)+"]"
# 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 = gri.register_gridfunctions("AUXEVOL","alpha_face")
gamma_faceDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gamma_faceDD","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)
import GiRaFFE_NRPy.Afield_flux as Af
expr_list = []
exprcheck_list = []
namecheck_list = []
for flux_dirn in range(3):
calculate_E_i_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,\
Valenciav_rU,B_rU,Valenciav_lU,B_lU)
Af.calculate_E_i_flux(flux_dirn,alpha_face,gamma_faceDD,beta_faceU,\
Valenciav_rU,B_rU,Valenciav_lU,B_lU)
namecheck_list.extend([gfnm("E_fluxD",flux_dirn)])
exprcheck_list.extend([Af.E_fluxD[flux_dirn]])
expr_list.extend([E_fluxD[flux_dirn]])
for mom_comp in range(len(expr_list)):
comp_func(expr_list[mom_comp],exprcheck_list[mom_comp],namecheck_list[mom_comp])
import sys
if all_passed:
print("ALL TESTS PASSED!")
else:
print("ERROR: AT LEAST ONE TEST DID NOT PASS")
sys.exit(1)
ALL TESTS PASSED!
We will also check the output C code to make sure it matches what is produced by the python module.
import difflib
import sys
subdir = os.path.join("RHSs")
out_dir = os.path.join("GiRaFFE_standalone_Ccodes")
cmd.mkdir(out_dir)
cmd.mkdir(os.path.join(out_dir,subdir))
valdir = os.path.join("GiRaFFE_Ccodes_validation")
cmd.mkdir(valdir)
cmd.mkdir(os.path.join(valdir,subdir))
generate_Afield_flux_function_files(out_dir,subdir,alpha_face,gamma_faceDD,beta_faceU,\
Valenciav_rU,B_rU,Valenciav_lU,B_lU,inputs_provided=True)
Af.generate_Afield_flux_function_files(valdir,subdir,alpha_face,gamma_faceDD,beta_faceU,\
Valenciav_rU,B_rU,Valenciav_lU,B_lU,inputs_provided=True)
print("Printing difference between original C code and this code...")
# Open the files to compare
files = ["RHSs/calculate_E_field_D0_right.h",
"RHSs/calculate_E_field_D0_left.h",
"RHSs/calculate_E_field_D1_right.h",
"RHSs/calculate_E_field_D1_left.h",
"RHSs/calculate_E_field_D2_right.h",
"RHSs/calculate_E_field_D2_left.h"]
for file in files:
print("Checking file " + file)
with open(os.path.join(valdir,file)) as file1, open(os.path.join(out_dir,file)) as file2:
# Read the lines of each file
file1_lines = file1.readlines()
file2_lines = file2.readlines()
num_diffs = 0
for line in difflib.unified_diff(file1_lines, file2_lines, fromfile=os.path.join(valdir,file), tofile=os.path.join(out_dir,file)):
sys.stdout.writelines(line)
num_diffs = num_diffs + 1
if num_diffs == 0:
print("No difference. TEST PASSED!")
else:
print("ERROR: Disagreement found with .py file. See differences above.")
Output C function calculate_E_field_D0_right() to file GiRaFFE_standalone_Ccodes\RHSs\calculate_E_field_D0_right.h Output C function calculate_E_field_D0_left() to file GiRaFFE_standalone_Ccodes\RHSs\calculate_E_field_D0_left.h Output C function calculate_E_field_D1_right() to file GiRaFFE_standalone_Ccodes\RHSs\calculate_E_field_D1_right.h Output C function calculate_E_field_D1_left() to file GiRaFFE_standalone_Ccodes\RHSs\calculate_E_field_D1_left.h Output C function calculate_E_field_D2_right() to file GiRaFFE_standalone_Ccodes\RHSs\calculate_E_field_D2_right.h Output C function calculate_E_field_D2_left() to file GiRaFFE_standalone_Ccodes\RHSs\calculate_E_field_D2_left.h Output C function calculate_E_field_D0_right() to file GiRaFFE_Ccodes_validation\RHSs\calculate_E_field_D0_right.h Output C function calculate_E_field_D0_left() to file GiRaFFE_Ccodes_validation\RHSs\calculate_E_field_D0_left.h Output C function calculate_E_field_D1_right() to file GiRaFFE_Ccodes_validation\RHSs\calculate_E_field_D1_right.h Output C function calculate_E_field_D1_left() to file GiRaFFE_Ccodes_validation\RHSs\calculate_E_field_D1_left.h Output C function calculate_E_field_D2_right() to file GiRaFFE_Ccodes_validation\RHSs\calculate_E_field_D2_right.h Output C function calculate_E_field_D2_left() to file GiRaFFE_Ccodes_validation\RHSs\calculate_E_field_D2_left.h Printing difference between original C code and this code... Checking file RHSs/calculate_E_field_D0_right.h No difference. TEST PASSED! Checking file RHSs/calculate_E_field_D0_left.h No difference. TEST PASSED! Checking file RHSs/calculate_E_field_D1_right.h No difference. TEST PASSED! Checking file RHSs/calculate_E_field_D1_left.h No difference. TEST PASSED! Checking file RHSs/calculate_E_field_D2_right.h No difference. TEST PASSED! Checking file RHSs/calculate_E_field_D2_left.h No difference. TEST PASSED!
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-Induction_Equation.pdf (Note that clicking on this link may not work; you may need to open the PDF file through another means.)
import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-GiRaFFE_NRPy-Afield_flux",location_of_template_file=os.path.join(".."))
Notebook output to PDF is only supported on Linux systems, with pdflatex installed.