Notebook

`GiRaFFE_NRPy`: A-to-B code¶

Author: Patrick Nelson¶

Notebook Status: Validated

Validation Notes: This code produces the expected magnetic fields for arbitrary vector potentials.

This module presents the functionality of GiRaFFE_NRPy_A2B.py.¶

Introduction:¶

This writes and documents the C code that GiRaFFE_NRPy uses to compute the magnetic fields from the vector potential. This is a relatively straightforward calculation, but requires care to be taken in the ghost zones.

We will need to compute $B^i$ everywhere in order to evolve $\tilde{S}_i$. However, $B^i = \epsilon^{ijk} \partial_j A_k$ requires derivatives of $A_i$, so getting $B^i$ in the ghostzones (and not just on the interior) will require some finesse. To do so, we will first compute the derivatives on the interior normally. Then, in the ghost zones, we will check if the point is on any face. If it is, we will use an appropriately-shifted stencil to avoid accessing points that are not on our grid. This will let us compute the derivative at the outermost gridpoints.

Table of Contents¶

$$\label{toc}$$

This notebook is organized as follows

Step 1: The Magnetic Field in terms of the Vector Potential
Step 2: Write out C code
1. Step 2.a: Compute the magnetic field at a point in a ghostzone
2. Step 2.b: Looping over all ghostzones
Step 3: Code Validation against original C code
Step 4: Output this notebook to $\LaTeX$-formatted PDF file

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)

import cmdline_helper as cmd
outdir = "GiRaFFE_NRPy/GiRaFFE_Ccode_validation/A2B/"
cmd.mkdir(outdir)

Step 1: The Magnetic Field in terms of the Vector Potential [Back to top]¶

$$\label{magnetic_field}$$

We start in the usual way - import the modules we need. We will also import the Levi-Civita symbol from indexedexp.py and use it to set the Levi-Civita tensor $\epsilon^{ijk} = [ijk]/\sqrt{\gamma}$.

In [2]:

# Step 1: The A-to-B driver
from outputC import outCfunction, lhrh # 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 indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support

thismodule = "GiRaFFE_NRPy_A2B"

CoordSystem     = "Cartesian"

# Set the finite-differencing order to 2
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", 2)

# Declare the three metric and compute the sqrt of its determinant.
gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01")
import GRHD.equations as gh
gh.compute_sqrtgammaDET(gammaDD)

# Import the Levi-Civita symbol and build the corresponding tensor.
# We already have a handy function to define the Levi-Civita symbol in indexedexp.py
LeviCivitaUUU = ixp.LeviCivitaTensorUUU_dim3_rank3(gh.sqrtgammaDET)

We will now create the structures we need: register the gridfunctions AD and BU, declare AD_dD for $A_{i,j}$ and zero BU for $B^i$. Then, we'll build the standard form of BU that we will use: $$B^i = \epsilon^{ijk} A_{k,j}.$$

In [3]:

AD = ixp.register_gridfunctions_for_single_rank1("EVOL","AD")
BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU")
AD_dD = ixp.declarerank2("AD_dD","nosym")
BU = ixp.zerorank1() # BU is already registered as a gridfunction, but we need to zero its values and declare it in this scope.
for i in range(3):
    for j in range(3):
        for k in range(3):
            BU[i] += LeviCivitaUUU[i][j][k] * AD_dD[k][j]

Step 2: Write out C code [Back to top]¶

$$\label{ccode}$$

Step 2.a: Compute the magnetic field at a point in a ghostzone [Back to top]¶

$$\label{ghostzones}$$

The rest of this file will consist of two functions. The first, compute_A2B_in_ghostzones(), will loop over the region specified by i0min, i0max, i1min, i1max, i2min, and i2max, passed to the function from the driver function. At each point, we consider the derivatives in each of the three directions that we will need to compute the curl; if we are on the edge of the grid, we will shift the stencil we use one point so that it no longer extends to outside our computational domain. Once we have all the derivatives we will need, we simply calculate the magnetic field as we otherwise would.

With the header files created to actually calculate $B^i$ from $A_k$, we will write the C code to direct how those are used. We will start by including needed header files. We will then define REAL as type double (so that we can easily change it everywhere if we want to use a different precision) and the IDX macros that we use to access specific points in the gridfunction arrays.

In [4]:

%%writefile $outdir/driver_AtoB.h
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#ifndef REAL
#define REAL double
#include "NGHOSTS.h" // A NRPy+-generated file, which is set based on FD_CENTDERIVS_ORDER.
#include "../CurviBoundaryConditions/gridfunction_defines.h"
#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )
#endif

Overwriting GiRaFFE_NRPy/GiRaFFE_Ccode_validation/A2B//driver_AtoB.h

Note the order of the if statements here; we first check if a point is not on an outermost-face because that is most likely. This will save time during execution by reducing the number of times the if statements are evaluated.

In [5]:

%%writefile $outdir/driver_AtoB.h
void compute_A2B_in_ghostzones(const paramstruct *restrict params,REAL *restrict in_gfs,REAL *restrict auxevol_gfs,
                                      const int i0min,const int i0max,
                                      const int i1min,const int i1max,
                                      const int i2min,const int i2max) {
#include "../set_Cparameters.h"
    for(int i2=i2min;i2<i2max;i2++) for(int i1=i1min;i1<i1max;i1++) for(int i0=i0min;i0<i0max;i0++) {
        REAL dx_Ay,dx_Az,dy_Ax,dy_Az,dz_Ax,dz_Ay;
        // Check to see if we're on the +x or -x face. If so, use a downwinded- or upwinded-stencil, respectively.
        // Otherwise, use a centered stencil.
        if (i0 > 0 && i0 < Nxx_plus_2NGHOSTS0-1) {
            dx_Ay = invdx0*(-1.0/2.0*in_gfs[IDX4S(AD1GF, i0-1,i1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD1GF, i0+1,i1,i2)]);
            dx_Az = invdx0*(-1.0/2.0*in_gfs[IDX4S(AD2GF, i0-1,i1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD2GF, i0+1,i1,i2)]);
        }
        else if (i0==0) {
            dx_Ay = invdx0*(-3.0/2.0*in_gfs[IDX4S(AD1GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD1GF, i0+1,i1,i2)] - 1.0/2.0*in_gfs[IDX4S(AD1GF, i0+2,i1,i2)]);
            dx_Az = invdx0*(-3.0/2.0*in_gfs[IDX4S(AD2GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD2GF, i0+1,i1,i2)] - 1.0/2.0*in_gfs[IDX4S(AD2GF, i0+2,i1,i2)]);
        }
        else {
            dx_Ay = invdx0*((3.0/2.0)*in_gfs[IDX4S(AD1GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD1GF, i0-1,i1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD1GF, i0-2,i1,i2)]);
            dx_Az = invdx0*((3.0/2.0)*in_gfs[IDX4S(AD2GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD2GF, i0-1,i1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD2GF, i0-2,i1,i2)]);
        }
        // As above, but in the y direction.
        if (i1 > 0 && i1 < Nxx_plus_2NGHOSTS1-1) {
            dy_Ax = invdx1*(-1.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1-1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1+1,i2)]);
            dy_Az = invdx1*(-1.0/2.0*in_gfs[IDX4S(AD2GF, i0,i1-1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD2GF, i0,i1+1,i2)]);
        }
        else if (i1==0) {
            dy_Ax = invdx1*(-3.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD0GF, i0,i1+1,i2)] - 1.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1+2,i2)]);
            dy_Az = invdx1*(-3.0/2.0*in_gfs[IDX4S(AD2GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD2GF, i0,i1+1,i2)] - 1.0/2.0*in_gfs[IDX4S(AD2GF, i0,i1+2,i2)]);
        }
        else {
            dy_Ax = invdx1*((3.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD0GF, i0,i1-1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1-2,i2)]);
            dy_Az = invdx1*((3.0/2.0)*in_gfs[IDX4S(AD2GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD2GF, i0,i1-1,i2)] + (1.0/2.0)*in_gfs[IDX4S(AD2GF, i0,i1-2,i2)]);
        }
        // As above, but in the z direction.
        if (i2 > 0 && i2 < Nxx_plus_2NGHOSTS2-1) {
            dz_Ax = invdx2*(-1.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1,i2-1)] + (1.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1,i2+1)]);
            dz_Ay = invdx2*(-1.0/2.0*in_gfs[IDX4S(AD1GF, i0,i1,i2-1)] + (1.0/2.0)*in_gfs[IDX4S(AD1GF, i0,i1,i2+1)]);
        }
        else if (i2==0) {
            dz_Ax = invdx2*(-3.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD0GF, i0,i1,i2+1)] - 1.0/2.0*in_gfs[IDX4S(AD0GF, i0,i1,i2+2)]);
            dz_Ay = invdx2*(-3.0/2.0*in_gfs[IDX4S(AD1GF, i0,i1,i2)] + 2*in_gfs[IDX4S(AD1GF, i0,i1,i2+1)] - 1.0/2.0*in_gfs[IDX4S(AD1GF, i0,i1,i2+2)]);
        }
        else {
            dz_Ax = invdx2*((3.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD0GF, i0,i1,i2-1)] + (1.0/2.0)*in_gfs[IDX4S(AD0GF, i0,i1,i2-2)]);
            dz_Ay = invdx2*((3.0/2.0)*in_gfs[IDX4S(AD1GF, i0,i1,i2)] - 2*in_gfs[IDX4S(AD1GF, i0,i1,i2-1)] + (1.0/2.0)*in_gfs[IDX4S(AD1GF, i0,i1,i2-2)]);
        }
        // Compute the magnetic field in the normal way, using the previously calculated derivatives.
        const double gammaDD00 = auxevol_gfs[IDX4S(GAMMADD00GF, i0,i1,i2)];
        const double gammaDD01 = auxevol_gfs[IDX4S(GAMMADD01GF, i0,i1,i2)];
        const double gammaDD02 = auxevol_gfs[IDX4S(GAMMADD02GF, i0,i1,i2)];
        const double gammaDD11 = auxevol_gfs[IDX4S(GAMMADD11GF, i0,i1,i2)];
        const double gammaDD12 = auxevol_gfs[IDX4S(GAMMADD12GF, i0,i1,i2)];
        const double gammaDD22 = auxevol_gfs[IDX4S(GAMMADD22GF, i0,i1,i2)];
        /*
        * NRPy+ Finite Difference Code Generation, Step 2 of 1: Evaluate SymPy expressions and write to main memory:
        */
        const double invsqrtg = 1.0/sqrt(gammaDD00*gammaDD11*gammaDD22
                                       - gammaDD00*gammaDD12*gammaDD12
                                       + 2*gammaDD01*gammaDD02*gammaDD12
                                       - gammaDD11*gammaDD02*gammaDD02
                                       - gammaDD22*gammaDD01*gammaDD01);
        auxevol_gfs[IDX4S(BU0GF, i0,i1,i2)] = (dy_Az-dz_Ay)*invsqrtg;
        auxevol_gfs[IDX4S(BU1GF, i0,i1,i2)] = (dz_Ax-dx_Az)*invsqrtg;
        auxevol_gfs[IDX4S(BU2GF, i0,i1,i2)] = (dx_Ay-dy_Ax)*invsqrtg;
    }
}

Overwriting GiRaFFE_NRPy/GiRaFFE_Ccode_validation/A2B//driver_AtoB.h

Step 2.b: Multiple-order algorithm for interior [Back to top]¶

$$\label{driver}$$

To help avoid Gibbs' phenomenon, we will implement an algorithm to lower the finite-differencing order near shocks. This will be done by calculating the magnetic field at multiple orders and using the lower order if the different orders disagree. What we write here will take the place of a NRPy+-generated 'body', as the parameter of the NRPy+ function outCfunction(). It will act on each point and each component individually. This will be much easier with some auxiliary functions, which we will write first. One will be a basic relative error function, another will compute a component of the magnetic field based on the input vector potential, and the last will return the best-choice order.

In [6]:

%%writefile -a $outdir/driver_AtoB.h
REAL relative_error(REAL a, REAL b) {
    if((a+b)!=0.0) {
        return 2.0*fabs(a-b)/fabs(a+b);
    }
    else {
        return 0.0;
    }
}

#define M2 0
#define M1 1
#define P0 2
#define P1 3
#define P2 4
#define CN4 0
#define CN2 1
#define UP2 2
#define DN2 3
#define UP1 4
#define DN1 5
void compute_Bx_pointwise(REAL *Bx, const REAL invdy, const REAL *Ay, const REAL invdz, const REAL *Az) {
    REAL dz_Ay,dy_Az;
    dz_Ay = invdz*((Ay[P1]-Ay[M1])*2.0/3.0 - (Ay[P2]-Ay[M2])/12.0);
    dy_Az = invdy*((Az[P1]-Az[M1])*2.0/3.0 - (Az[P2]-Az[M2])/12.0);
    Bx[CN4] = dy_Az - dz_Ay;

    dz_Ay = invdz*(Ay[P1]-Ay[M1])/2.0;
    dy_Az = invdy*(Az[P1]-Az[M1])/2.0;
    Bx[CN2] = dy_Az - dz_Ay;

    dz_Ay = invdz*(-1.5*Ay[P0]+2.0*Ay[P1]-0.5*Ay[P2]);
    dy_Az = invdy*(-1.5*Az[P0]+2.0*Az[P1]-0.5*Az[P2]);
    Bx[UP2] = dy_Az - dz_Ay;

    dz_Ay = invdz*(1.5*Ay[P0]-2.0*Ay[M1]+0.5*Ay[M2]);
    dy_Az = invdy*(1.5*Az[P0]-2.0*Az[M1]+0.5*Az[M2]);
    Bx[DN2] = dy_Az - dz_Ay;

    dz_Ay = invdz*(Ay[P1]-Ay[P0]);
    dy_Az = invdy*(Az[P1]-Az[P0]);
    Bx[UP1] = dy_Az - dz_Ay;

    dz_Ay = invdz*(Ay[P0]-Ay[M1]);
    dy_Az = invdy*(Az[P0]-Az[M1]);
    Bx[DN1] = dy_Az - dz_Ay;
}

#define TOLERANCE_A2B 1.0e-4
REAL find_accepted_Bx_order(REAL *Bx) {
    REAL accepted_val = Bx[CN4];
    REAL Rel_error_o2_vs_o4 = relative_error(Bx[CN2],Bx[CN4]);
    REAL Rel_error_oCN2_vs_oDN2 = relative_error(Bx[CN2],Bx[DN2]);
    REAL Rel_error_oCN2_vs_oUP2 = relative_error(Bx[CN2],Bx[UP2]);

    if(Rel_error_o2_vs_o4 > TOLERANCE_A2B) {
        accepted_val = Bx[CN2];
        if(Rel_error_o2_vs_o4 > Rel_error_oCN2_vs_oDN2 || Rel_error_o2_vs_o4 > Rel_error_oCN2_vs_oUP2) {
            // Should we use AND or OR in if statement?
            if(relative_error(Bx[UP2],Bx[UP1]) < relative_error(Bx[DN2],Bx[DN1])) {
                accepted_val = Bx[UP2];
            }
            else {
                accepted_val = Bx[DN2];
            }
        }
    }
    return accepted_val;
}

Appending to GiRaFFE_NRPy/GiRaFFE_Ccode_validation/A2B//driver_AtoB.h

The algorithm will start by reading in all the quantities it will need and storing them as temporary variables. This includes all three components of the vector potential in two different directions each, reflecting the derivatives we will need to take. Five points are needed in each direction in order to fill the templates for fourth- and second-order centered-differencing and first-order forward- and backward-differencing. We calculate the inverse of the square root of the metric determinant and then compute all three components of the magnetic field at all orders. Then, if the relative differences between orders are found to be high, the order is dropped. The final result is written to memory.

*This is currently not used, it does not reliably improve shock handling*

In [7]:

# order_lowering_body = """REAL AD0_1[5],AD0_2[5],AD1_2[5],AD1_0[5],AD2_0[5],AD2_1[5];
# const double gammaDD00 = auxevol_gfs[IDX4S(GAMMADD00GF, i0,i1,i2)];
# const double gammaDD01 = auxevol_gfs[IDX4S(GAMMADD01GF, i0,i1,i2)];
# const double gammaDD02 = auxevol_gfs[IDX4S(GAMMADD02GF, i0,i1,i2)];
# const double gammaDD11 = auxevol_gfs[IDX4S(GAMMADD11GF, i0,i1,i2)];
# const double gammaDD12 = auxevol_gfs[IDX4S(GAMMADD12GF, i0,i1,i2)];
# const double gammaDD22 = auxevol_gfs[IDX4S(GAMMADD22GF, i0,i1,i2)];
# AD0_2[M2] = in_gfs[IDX4S(AD0GF, i0,i1,i2-2)];
# AD0_2[M1] = in_gfs[IDX4S(AD0GF, i0,i1,i2-1)];
# AD0_1[M2] = in_gfs[IDX4S(AD0GF, i0,i1-2,i2)];
# AD0_1[M1] = in_gfs[IDX4S(AD0GF, i0,i1-1,i2)];
# AD0_1[P0] = AD0_2[P0] = in_gfs[IDX4S(AD0GF, i0,i1,i2)];
# AD0_1[P1] = in_gfs[IDX4S(AD0GF, i0,i1+1,i2)];
# AD0_1[P2] = in_gfs[IDX4S(AD0GF, i0,i1+2,i2)];
# AD0_2[P1] = in_gfs[IDX4S(AD0GF, i0,i1,i2+1)];
# AD0_2[P2] = in_gfs[IDX4S(AD0GF, i0,i1,i2+2)];
# AD1_2[M2] = in_gfs[IDX4S(AD1GF, i0,i1,i2-2)];
# AD1_2[M1] = in_gfs[IDX4S(AD1GF, i0,i1,i2-1)];
# AD1_0[M2] = in_gfs[IDX4S(AD1GF, i0-2,i1,i2)];
# AD1_0[M1] = in_gfs[IDX4S(AD1GF, i0-1,i1,i2)];
# AD1_2[P0] = AD1_0[P0] = in_gfs[IDX4S(AD1GF, i0,i1,i2)];
# AD1_0[P1] = in_gfs[IDX4S(AD1GF, i0+1,i1,i2)];
# AD1_0[P2] = in_gfs[IDX4S(AD1GF, i0+2,i1,i2)];
# AD1_2[P1] = in_gfs[IDX4S(AD1GF, i0,i1,i2+1)];
# AD1_2[P2] = in_gfs[IDX4S(AD1GF, i0,i1,i2+2)];
# AD2_1[M2] = in_gfs[IDX4S(AD2GF, i0,i1-2,i2)];
# AD2_1[M1] = in_gfs[IDX4S(AD2GF, i0,i1-1,i2)];
# AD2_0[M2] = in_gfs[IDX4S(AD2GF, i0-2,i1,i2)];
# AD2_0[M1] = in_gfs[IDX4S(AD2GF, i0-1,i1,i2)];
# AD2_0[P0] = AD2_1[P0] = in_gfs[IDX4S(AD2GF, i0,i1,i2)];
# AD2_0[P1] = in_gfs[IDX4S(AD2GF, i0+1,i1,i2)];
# AD2_0[P2] = in_gfs[IDX4S(AD2GF, i0+2,i1,i2)];
# AD2_1[P1] = in_gfs[IDX4S(AD2GF, i0,i1+1,i2)];
# AD2_1[P2] = in_gfs[IDX4S(AD2GF, i0,i1+2,i2)];
# const double invsqrtg = 1.0/sqrt(gammaDD00*gammaDD11*gammaDD22
#                                - gammaDD00*gammaDD12*gammaDD12
#                                + 2*gammaDD01*gammaDD02*gammaDD12
#                                - gammaDD11*gammaDD02*gammaDD02
#                                - gammaDD22*gammaDD01*gammaDD01);

# REAL BU0[4],BU1[4],BU2[4];
# compute_Bx_pointwise(BU0,invdx2,AD1_2,invdx1,AD2_1);
# compute_Bx_pointwise(BU1,invdx0,AD2_0,invdx2,AD0_2);
# compute_Bx_pointwise(BU2,invdx1,AD0_1,invdx0,AD1_0);

# auxevol_gfs[IDX4S(BU0GF, i0,i1,i2)] = find_accepted_Bx_order(BU0)*invsqrtg;
# auxevol_gfs[IDX4S(BU1GF, i0,i1,i2)] = find_accepted_Bx_order(BU1)*invsqrtg;
# auxevol_gfs[IDX4S(BU2GF, i0,i1,i2)] = find_accepted_Bx_order(BU2)*invsqrtg;
# """

Step 2.c: Looping over all ghostzones [Back to top]¶

$$\label{driver}$$

This function is responsible for driving the logic needed to compute $B^i$ from $A_i$ whenever it is needed in the main code (typically done along with applying BCs after a timestep in the main code). We first loop over the interior and calculate $B^i$ there at the FD order currently set in NRPy+. We then set imin and imax such that they define the interior of the grid.

With the Interior handled, we call compute_A2B_in_ghostzones several times. Note that imin and imax are passed such that each call loops over one face in one ghostzone (e.g., the innermost ghostzone on the $+x$ face). On the first pass through the loop, we calculate $B^i$ in the innermost ghostzone. We also decrement imin and increment imax on each pass through the loop. This serves two purposes: the next time through the loop, compute_A2B_in_ghostzones() will act on the next-outer ghost zone. on the next-outer ghost zone, and the size of each face will increase, eventually allowing us to act on each point.

In [8]:

desc="Compute the magnetic field from the vector potential everywhere, including ghostzones"
#     body     = order_lowering_body,
name="driver_A_to_B"
driver_Ccode = outCfunction(
    outfile  = "returnstring", desc=desc, name=name,
    params   = "const paramstruct *restrict params,REAL *restrict in_gfs,REAL *restrict auxevol_gfs",
    body     = fin.FD_outputC("returnstring",[lhrh(lhs=gri.gfaccess("out_gfs","BU0"),rhs=BU[0]),\
                                              lhrh(lhs=gri.gfaccess("out_gfs","BU1"),rhs=BU[1]),\
                                              lhrh(lhs=gri.gfaccess("out_gfs","BU2"),rhs=BU[2])]),
    postloop = """
    int imin[3] = { NGHOSTS_A2B, NGHOSTS_A2B, NGHOSTS_A2B };
    int imax[3] = { NGHOSTS+Nxx0, NGHOSTS+Nxx1, NGHOSTS+Nxx2 };
    // Now, we loop over the ghostzones to calculate the magnetic field there.
    for(int which_gz = 0; which_gz < NGHOSTS_A2B; which_gz++) {
        // After updating each face, adjust imin[] and imax[]
        //   to reflect the newly-updated face extents.
        compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imin[0]-1,imin[0], imin[1],imax[1], imin[2],imax[2]); imin[0]--;
        compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imax[0],imax[0]+1, imin[1],imax[1], imin[2],imax[2]); imax[0]++;

        compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imin[0],imax[0], imin[1]-1,imin[1], imin[2],imax[2]); imin[1]--;
        compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imin[0],imax[0], imax[1],imax[1]+1, imin[2],imax[2]); imax[1]++;

        compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imin[0],imax[0], imin[1],imax[1], imin[2]-1,imin[2]); imin[2]--;
        compute_A2B_in_ghostzones(params,in_gfs,auxevol_gfs,imin[0],imax[0], imin[1],imax[1], imax[2],imax[2]+1); imax[2]++;
    }
""",
    loopopts="InteriorPoints",
    rel_path_to_Cparams=os.path.join("../")).replace("= NGHOSTS","= NGHOSTS_A2B").replace("NGHOSTS+Nxx0","Nxx_plus_2NGHOSTS0-NGHOSTS_A2B").replace("NGHOSTS+Nxx1","Nxx_plus_2NGHOSTS1-NGHOSTS_A2B").replace("NGHOSTS+Nxx2","Nxx_plus_2NGHOSTS2-NGHOSTS_A2B")

with open(os.path.join(outdir,"driver_AtoB.h"),"a") as file:
    file.write(driver_Ccode)

Step 3: Code Validation against original C code [Back to top]¶

$$\label{code_validation}$$

To validate the code in this tutorial we check for agreement between the files

that were written in this tutorial and
those that are stored in GiRaFFE_NRPy/GiRaFFE_Ccode_library or generated by GiRaFFE_NRPy_A2B.py

In [9]:

# Define the directory that we wish to validate against:
valdir = "GiRaFFE_NRPy/GiRaFFE_Ccode_library/A2B/"

import GiRaFFE_NRPy.GiRaFFE_NRPy_A2B as A2B
A2B.GiRaFFE_NRPy_A2B(valdir,gammaDD,AD,BU)

import difflib
import sys

print("Printing difference between original C code and this code...")
# Open the files to compare
file = "driver_AtoB.h"

print("Checking file " + file)
with open(os.path.join(valdir,file)) as file1, open(os.path.join(outdir,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(outdir+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.")

Printing difference between original C code and this code...
Checking file driver_AtoB.h
No difference. TEST PASSED!

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

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

In [10]:

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-GiRaFFE_NRPy-A2B",location_of_template_file=os.path.join(".."))

Created Tutorial-GiRaFFE_NRPy-A2B.tex, and compiled LaTeX file to PDF file
    Tutorial-GiRaFFE_NRPy-A2B.pdf

GiRaFFE_NRPy: A-to-B code¶