Notebook

The Piecewise-Parabolic Method¶

Author: Patrick Nelson¶

Notebook Status: Validated

Validation Notes: This module is self-validated against its module and is also validated against the corresponding algorithm in the old GiRaFFE code in this tutorial.

This module presents the functionality of GiRaFFE_NRPy_PPM.py .¶

This notebook documents and ports the function from the original GiRaFFE that implements the reconstruction algorithm used by the piecewise-parabolic method (PPM) of Colella and Woodward (1984).

Introduction:¶

The differential equations that GiRaFFE evolves have two different terms that contribute to the time evolution of some quantity: the flux term and the source term. The flux terms will require us to use Riemann solvers to properly handle shocks in the system without Gibbs-like phenomena; such algorithms require data at cell interfaces. For a cell centered at point $(i,j,k)$ and reconstructed in $i$, these are at positions $(i \pm 1/2,j,k)$, and so forth for reconstructions in other directions (the direction in which the code does the reconstruction is specifiable as the input parameter flux_dirn). This reconstruction algorithm will compute the value of input quantities to the left and right side of the interface at $i-1/2$. It thus follows that the data at the interface $i+1/2$ can be found by accessing $i+1$ in the code. The algorithm will also apply slope-limiting and monotization algorithms that will serve to preserve stability near shocks.

This algorithm is not quite as accessible as the much simpler finite-difference methods; as such, this notebook is recommended as an introduction. It covers a simpler reconstruction scheme, and proved useful in preparing the documentation for this more complicated scheme.

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

The Reconstruction Step (This notebook)
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
2. Use that fit to define the state at the left and right interface of each cell
3. Apply a slope limiter to mitigate Gibbs phenomenon
Solving the Riemann Problem
1. Use the left and right reconstructed states to calculate the unique state at boundary
2. Use the unique state to estimate the derivative in the cell
Repeat the above for each conservative gridfunction in each direction

Table of Contents¶

$$\label{toc}$$

This notebook is organized as follows

Step 0: Preliminaries
Step 1: The reconstruction function
1. Step 1.a: Some definitions and declarations
2. Step 1.b: The function definition
3. Step 1.c: Interpolate the face values
4. Step 1.d: Monotonize the values within each cell
5. Step 1.e: Shift indices
Step 2: The slope limiter
Step 3: The monotonization algorithm
Step 4: Dependencies
Step 5: Code Validation against GiRaFFE_NRPy_PPM.py
Step 6: Example code in order to use this algorithm
Step 7: Output this notebook to $\LaTeX$-formatted PDF file

Step 0: Preliminaries [Back to top]¶

$$\label{prelim}$$

This first block of code just sets up a subdirectory within GiRaFFE_standalone_Ccodes/ to which we will write the C code.

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     # NRPy+: Multi-platform Python command-line interface
Ccodesdir = "GiRaFFE_standalone_Ccodes/PPM"
cmd.mkdir(os.path.join(Ccodesdir))

The structure gf_and_gz_struct is a C++ structure used to keep track of ghostzone information between routines. It contains a pointer and two arrays. It is specified by the following code:

// Keeping track of ghostzones between routines is a nightmare, so
//   we instead attach ghostzone info to each gridfunction and set
//   the ghostzone information correctly within each routine.
struct gf_and_gz_struct {
  REAL *gf;
  int gz_lo[4],gz_hi[4];
};

In NRPy+, we will have to interact with our arrays on a lower level than normal in order to get this pointer. NRPy+ stores gridfunctions all in one array, but values for a specific quantity (e.g., ValenciavU0) are still stored contiguously. That means that the pointer we are after is simply the one that points to the first point of that quantity, i.e. i0=0, i1=0, and i2=0. Recall the typical macro we use for memory access:

#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )

It thus follows that the pointer that needs to be stored in this structure is Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2*g, where g is, in this example, VALENCIAVU0GF.

Step 1: The reconstruction function [Back to top]¶

$$\label{reconstruction}$$

Step 1.a: Some definitions and declarations [Back to top]¶

$$\label{define}$$

This file contains the functions necessary for reconstruction. It is based on Colella & Woodward PPM in the case where pressure and density $P = \rho = 0$.

We start by defining the values of MINUS2...PLUS2 as $\{0, \ldots ,4\}$ for the sake of convenience later on, when we will read data from positions $i-2$ and $i+2$ to store in an array; we also define MAXNUMINDICES as 5 so we can easily loop over the above. We include loop_defines_reconstruction_NRPy.h for some macros that will allow us to conveniently write common loops that we will use and give the function prototypes for our slope limiter, slope_limit_NRPy(), and our monotization algorithm, monotonize_NRPy().

In [2]:

%%writefile $Ccodesdir/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c
/*****************************************
 * PPM Reconstruction Interface.
 * Zachariah B. Etienne (2013)
 *
 * This version of PPM implements the standard
 * Colella & Woodward PPM, but in the GRFFE
 * limit, where P=rho=0. Thus, e.g., ftilde=0.
 *****************************************/

#define MINUS2 0
#define MINUS1 1
#define PLUS0  2
#define PLUS1  3
#define PLUS2  4
#define MAXNUMINDICES 5
//    ^^^^^^^^^^^^^ Be _sure_ to define MAXNUMINDICES appropriately!

#ifndef MIN
#define MIN(a,b) ( ((a) < (b)) ? (a) : (b) )
#endif
#ifndef MIN
#define MAX(a,b) ( ((a) > (b)) ? (a) : (b) )
#endif
#ifndef SQR
#define SQR(x) ((x) * (x))
#endif
// FIXME: Should make this zero-offset for NRPy+ standards. Probably a wrapper function for compatibility with a minimum of other changes?
const int kronecker_delta[4][3] = { { 0,0,0 },
                                    { 1,0,0 },
                                    { 0,1,0 },
                                    { 0,0,1 } };

// You'll find the #define's for LOOP_DEFINE and SET_INDEX_ARRAYS_NRPY inside:
#include "loop_defines_reconstruction_NRPy.h"

static inline REAL slope_limit_NRPy(const REAL *dU,const REAL *dUp1);
static inline void monotonize_NRPy(const REAL *U,REAL *Ur,REAL *Ul);

Overwriting GiRaFFE_standalone_Ccodes/PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

Step 1.b: The function definition [Back to top]¶

$$\label{func}$$

Here, we start the function definition for the main function for our reconstruction, reconstruct_set_of_prims_PPM_GRFFE_NRPy(). Among its parameters are the arrays that define the grid (that will need to be replaced with NRPy+ equivalents), a flux direction, the integer array specifying which primitives to reconstruct (as well as the number of primitives to reconstruct), the input structure in_prims, the output structures out_prims_r and out_prims_l, and a temporary array (this will be used to help switch variable names).

We then check the number of ghostzones and error out if there are too few - this method requires three. Note the for loop here; it continues through the next two cells as well, looping over each primitive we will reconstruct in the chosen direction.

In [3]:

%%writefile -a $Ccodesdir/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

void reconstruct_set_of_prims_PPM_GRFFE_NRPy(const paramstruct *params,REAL *auxevol_gfs,const int flux_dirn,
                                               const int num_prims_to_reconstruct,const int *which_prims_to_reconstruct,
                                               const gf_and_gz_struct *in_prims,gf_and_gz_struct *out_prims_r,
                                               gf_and_gz_struct *out_prims_l,REAL *temporary) {
#include "../set_Cparameters.h"
  const int Nxx_plus_2NGHOSTS[3] = {Nxx_plus_2NGHOSTS0,Nxx_plus_2NGHOSTS1,Nxx_plus_2NGHOSTS2};

  REAL U[NUM_RECONSTRUCT_GFS][MAXNUMINDICES],dU[NUM_RECONSTRUCT_GFS][MAXNUMINDICES],slope_lim_dU[NUM_RECONSTRUCT_GFS][MAXNUMINDICES],
    Ur[NUM_RECONSTRUCT_GFS][MAXNUMINDICES],Ul[NUM_RECONSTRUCT_GFS][MAXNUMINDICES];
  int ijkgz_lo_hi[4][2];

  for(int ww=0;ww<num_prims_to_reconstruct;ww++) {
    const int whichvar=which_prims_to_reconstruct[ww];

    if(in_prims[whichvar].gz_lo[flux_dirn]!=0 || in_prims[whichvar].gz_hi[flux_dirn]!=0) {
      printf("TOO MANY GZ'S! WHICHVAR=%d: %d %d %d : %d %d %d DIRECTION %d",whichvar,
		  in_prims[whichvar].gz_lo[1],in_prims[whichvar].gz_lo[2],in_prims[whichvar].gz_lo[3],
		  in_prims[whichvar].gz_hi[1],in_prims[whichvar].gz_hi[2],in_prims[whichvar].gz_hi[3],flux_dirn);
      exit(0);
    }

Appending to GiRaFFE_standalone_Ccodes/PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

Step 1.c: Interpolate the face values [Back to top]¶

$$\label{face}$$

In Loop 1, we will interpolate the face values at the left and right interfaces, Ur and Ul, respectively. This is done on a point-by-point basis as defined by the LOOP_DEFINE.

After reading in the relevant values from memory, we calculate the simple dU: \begin{align} dU_{-1} &= U_{-1} - U_{-2} \\ dU_{+0} &= U_{+0} - U_{-1} \\ dU_{+1} &= U_{+1} - U_{+0} \\ dU_{+2} &= U_{+2} - U_{+1}. \\ \end{align} From that, we compute the slope-limited slope_lim_dU, or $\nabla U$ (see below). Then, we compute the face values using eq. A1 from arxiv:astro-ph/050342, adapted from 1.9 in Colella and Woodward (1984): \begin{align} U_r &= \frac{1}{2} \left( U_{+1} + U_{+0} \right) + \frac{1}{6} \left( \nabla U_{+0} - \nabla U_{+1} \right) \\ U_l &= \frac{1}{2} \left( U_{+0} + U_{-1} \right) + \frac{1}{6} \left( \nabla U_{-1} - \nabla U_{+0} \right). \\ \end{align} (Note, however, that we use the standard coefficient $1/6$ instead of $1/8$.) Finally, we write the values to memory in the output structures.

In [4]:

%%writefile -a $Ccodesdir/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

    // *** LOOP 1: Interpolate to Ur and Ul, which are face values ***
    //  You will find that Ur depends on U at MINUS1,PLUS0, PLUS1,PLUS2, and
    //                     Ul depends on U at MINUS2,MINUS1,PLUS0,PLUS1.
    //  However, we define the below loop from MINUS2 to PLUS2. Why not split
    //     this up and get additional points? Maybe we should. In GRMHD, the
    //     reason is that later on, Ur and Ul depend on ftilde, which is
    //     defined from MINUS2 to PLUS2, so we would lose those points anyway.
    //     But in GRFFE, ftilde is set to zero, so there may be a potential
    //     for boosting performance here.
    LOOP_DEFINE(2,2,  Nxx_plus_2NGHOSTS,flux_dirn,  ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {
      SET_INDEX_ARRAYS_NRPY(-2,2,flux_dirn);
      /* *** LOOP 1a: READ INPUT *** */
      // Read in a primitive at all gridpoints between m = MINUS2 & PLUS2, where m's direction is given by flux_dirn. Store to U.
      for(int ii=MINUS2;ii<=PLUS2;ii++) U[whichvar][ii] = in_prims[whichvar].gf[index_arr[flux_dirn][ii]];

      /* *** LOOP 1b: DO COMPUTATION *** */
      /* First, compute simple dU = U(i) - U(i-1), where direction of i
       *         is given by flux_dirn, and U is a primitive variable:
       *         {vx,vy,vz,Bx,By,Bz}. */
      // Note that for Ur and Ul at i, we must compute dU(i-1),dU(i),dU(i+1),
      //         and dU(i+2)
      dU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2];
      dU[whichvar][PLUS0]  = U[whichvar][PLUS0] - U[whichvar][MINUS1];
      dU[whichvar][PLUS1]  = U[whichvar][PLUS1] - U[whichvar][PLUS0];
      dU[whichvar][PLUS2]  = U[whichvar][PLUS2] - U[whichvar][PLUS1];

      // Then, compute slope-limited dU, using MC slope limiter:
      slope_lim_dU[whichvar][MINUS1]=slope_limit_NRPy(&dU[whichvar][MINUS1],&dU[whichvar][PLUS0]);
      slope_lim_dU[whichvar][PLUS0] =slope_limit_NRPy(&dU[whichvar][PLUS0], &dU[whichvar][PLUS1]);
      slope_lim_dU[whichvar][PLUS1] =slope_limit_NRPy(&dU[whichvar][PLUS1], &dU[whichvar][PLUS2]);

      // Finally, compute face values Ur and Ul based on the PPM prescription
      //   (Eq. A1 in http://arxiv.org/pdf/astro-ph/0503420.pdf, but using standard 1/6=(1.0/6.0) coefficient)
      // Ur[PLUS0] represents U(i+1/2)
      // We applied a simplification to the following line: Ur=U+0.5*(U(i+1)-U) + ... = 0.5*(U(i+1)+U) + ...
      Ur[whichvar][PLUS0] = 0.5*(U[whichvar][PLUS1] + U[whichvar][PLUS0] ) + (1.0/6.0)*(slope_lim_dU[whichvar][PLUS0]  - slope_lim_dU[whichvar][PLUS1]);
      // Ul[PLUS0] represents U(i-1/2)
      // We applied a simplification to the following line: Ul=U(i-1)+0.5*(U-U(i-1)) + ... = 0.5*(U+U(i-1)) + ...
      Ul[whichvar][PLUS0] = 0.5*(U[whichvar][PLUS0] + U[whichvar][MINUS1]) + (1.0/6.0)*(slope_lim_dU[whichvar][MINUS1] - slope_lim_dU[whichvar][PLUS0]);

      /* *** LOOP 1c: WRITE OUTPUT *** */
      // Store right face values to {vxr,vyr,vzr,Bxr,Byr,Bzr},
      //    and left face values to {vxl,vyl,vzl,Bxl,Byl,Bzl}
      out_prims_r[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ur[whichvar][PLUS0];
      out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ul[whichvar][PLUS0];
    }

Appending to GiRaFFE_standalone_Ccodes/PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

Step 1.d: Monotonize the values within each cell [Back to top]¶

$$\label{monotonize}$$

We skip Loop 2 in GRFFE because it modifies variables that are 0 here by definition; then, we flatten the data in Loop 3 (but since we flatten based on ftilde_gf, which is 0 in GRFFE, we again don't really do anything). Also in Loop 3, we call the monotonize_NRPy() function on the face values. This function adjusts the face values to ensure that the data is monotonic within each cell to avoid the Gibbs phenomenon.

In [5]:

%%writefile -a $Ccodesdir/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

    // *** LOOP 2 (REMOVED): STEEPEN RHOB. RHOB DOES NOT EXIST IN GRFFE EQUATIONS ***
  }

  // *** LOOP 3: FLATTEN BASED ON FTILDE AND MONOTONIZE ***
  for(int ww=0;ww<num_prims_to_reconstruct;ww++) {
    const int whichvar=which_prims_to_reconstruct[ww];
    // ftilde() depends on P(MINUS2,MINUS1,PLUS1,PLUS2), THUS IS SET TO ZERO IN GRFFE
    LOOP_DEFINE(2,2,  Nxx_plus_2NGHOSTS,flux_dirn,  ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {
      SET_INDEX_ARRAYS_NRPY(0,0,flux_dirn);

      U[whichvar][PLUS0]  = in_prims[whichvar].gf[index_arr[flux_dirn][PLUS0]];
      Ur[whichvar][PLUS0] = out_prims_r[whichvar].gf[index_arr[flux_dirn][PLUS0]];
      Ul[whichvar][PLUS0] = out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]];

      // ftilde_gf was computed in the function compute_ftilde_gf(), called before this routine
      //REAL ftilde = ftilde_gf[index_arr[flux_dirn][PLUS0]];
      // ...and then flatten (local operation)
      Ur[whichvar][PLUS0]   = Ur[whichvar][PLUS0];
      Ul[whichvar][PLUS0]   = Ul[whichvar][PLUS0];

      // Then monotonize
      monotonize_NRPy(&U[whichvar][PLUS0],&Ur[whichvar][PLUS0],&Ul[whichvar][PLUS0]);

      out_prims_r[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ur[whichvar][PLUS0];
      out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ul[whichvar][PLUS0];
    }
    // Note: ftilde=0 in GRFFE. Ur depends on ftilde, which depends on points of U between MINUS2 and PLUS2
    out_prims_r[whichvar].gz_lo[flux_dirn]+=2;
    out_prims_r[whichvar].gz_hi[flux_dirn]+=2;
    // Note: ftilde=0 in GRFFE. Ul depends on ftilde, which depends on points of U between MINUS2 and PLUS2
    out_prims_l[whichvar].gz_lo[flux_dirn]+=2;
    out_prims_l[whichvar].gz_hi[flux_dirn]+=2;
  }

Appending to GiRaFFE_standalone_Ccodes/PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

Step 1.e: Shift indices [Back to top]¶

$$\label{shift}$$

In Loop 4, we will shift the indices of Ur and Ul. So far, we have been concerned with the behavior of the data within a single cell. In that context, it makes sense to call the value of data at the left end of the cell Ul and the data at the right end of the cell Ur. However, going forward, we will be concerned about the behavior of the data at the interface between cells. In this context, it sense to call the value of data on the left of the interface (which is at the right end of the cell!) Ul and the data on the right of the interface Ur. So, using the array temporary, we switch the two names while shifting Ur appropriately.

In [6]:

%%writefile -a $Ccodesdir/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

// *** LOOP 4: SHIFT Ur AND Ul ***
  /* Currently face values are set so that
   *      a) Ur(i) represents U(i+1/2), and
   *      b) Ul(i) represents U(i-1/2)
   *    Here, we shift so that the indices are consistent:
   *      a) U(i-1/2+epsilon) = oldUl(i)   = newUr(i)
   *      b) U(i-1/2-epsilon) = oldUr(i-1) = newUl(i)
   *    Note that this step is not strictly necessary if you keep
   *      track of indices when computing the flux. */
  for(int ww=0;ww<num_prims_to_reconstruct;ww++) {
    const int whichvar=which_prims_to_reconstruct[ww];
    LOOP_DEFINE(3,2,  Nxx_plus_2NGHOSTS,flux_dirn,  ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {
      SET_INDEX_ARRAYS_NRPY(-1,0,flux_dirn);
      temporary[index_arr[flux_dirn][PLUS0]] = out_prims_r[whichvar].gf[index_arr[flux_dirn][MINUS1]];
    }

    LOOP_DEFINE(3,2,  Nxx_plus_2NGHOSTS,flux_dirn,  ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {
      SET_INDEX_ARRAYS_NRPY(0,0,flux_dirn);
      // Then shift so that Ur represents the gridpoint at i-1/2+epsilon,
      //                and Ul represents the gridpoint at i-1/2-epsilon.
      // Ur(i-1/2) = Ul(i-1/2)     = U(i-1/2+epsilon)
      // Ul(i-1/2) = Ur(i+1/2 - 1) = U(i-1/2-epsilon)
      out_prims_r[whichvar].gf[index_arr[flux_dirn][PLUS0]] = out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]];
      out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = temporary[index_arr[flux_dirn][PLUS0]];
    }
    // Ul was just shifted, so we lost another ghostzone.
    out_prims_l[whichvar].gz_lo[flux_dirn]+=1;
    out_prims_l[whichvar].gz_hi[flux_dirn]+=0;
    // As for Ur, we didn't need to get rid of another ghostzone,
    //    but we did ... seems wasteful!
    out_prims_r[whichvar].gz_lo[flux_dirn]+=1;
    out_prims_r[whichvar].gz_hi[flux_dirn]+=0;

  }
}

Appending to GiRaFFE_standalone_Ccodes/PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

Step 2: The slope limiter [Back to top]¶

$$\label{slope_limit}$$

The first function here implements the Monotonized Central (MC) reconstruction slope limiter: $$ MC(a,b) = \left \{ \begin{array}{ll} 0 & {\rm if} ab \leq 0 \\ {\rm sign}(a) \min(2|a|,2|b|, |a+b|/2) & {\rm otherwise.} \end{array} \right. $$

This is adapted from eq. 1.8 of Colella and Woodward (1984).

In [7]:

%%writefile -a $Ccodesdir/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

// Set SLOPE_LIMITER_COEFF = 2.0 for MC, 1 for minmod
#define SLOPE_LIMITER_COEFF 2.0

//Eq. 60 in JOURNAL OF COMPUTATIONAL PHYSICS 123, 1-14 (1996)
//   [note the factor of 2 missing in the |a_{j+1} - a_{j}| term].
//   Recall that dU = U_{i} - U_{i-1}.
static inline REAL slope_limit_NRPy(const REAL *dU,const REAL *dUp1) {
  if((*dU)*(*dUp1) > 0.0) {
    //delta_m_U=0.5 * [ (u_(i+1)-u_i) + (u_i-u_(i-1)) ] = (u_(i+1) - u_(i-1))/2  <-- first derivative, second-order; this should happen most of the time (smooth flows)
    const REAL delta_m_U = 0.5*((*dU) + (*dUp1));
    // EXPLANATION OF BELOW LINE OF CODE.
    // In short, sign_delta_a_j = sign(delta_m_U) = (0.0 < delta_m_U) - (delta_m_U < 0.0).
    //    If delta_m_U>0, then (0.0 < delta_m_U)==1, and (delta_m_U < 0.0)==0, so sign_delta_a_j=+1
    //    If delta_m_U<0, then (0.0 < delta_m_U)==0, and (delta_m_U < 0.0)==1, so sign_delta_a_j=-1
    //    If delta_m_U==0,then (0.0 < delta_m_U)==0, and (delta_m_U < 0.0)==0, so sign_delta_a_j=0
    const int sign_delta_m_U = (0.0 < delta_m_U) - (delta_m_U < 0.0);
    //Decide whether to use 2nd order derivative or first-order derivative, limiting slope.
    return sign_delta_m_U*MIN(fabs(delta_m_U),MIN(SLOPE_LIMITER_COEFF*fabs(*dUp1),SLOPE_LIMITER_COEFF*fabs(*dU)));
  }
  return 0.0;
}

Appending to GiRaFFE_standalone_Ccodes/PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

Step 3: The monotonization algorithm [Back to top]¶

$$\label{monotonize_def}$$

The next function monotonizes the slopes using the algorithm from Colella and Woodward (1984), eq. 1.10. We want the slope to be monotonic in a cell in order to reduce the impact of the Gibbs phenomenon. So, we consider three values in the cell: the cell average, U; on the left interface of the cell, Ul; and on the right interface of the cell, Ur. The goal of the algorithm is to ensure monotonicity; so, it first checks to see if the cell contains a local extremum. If it does, we make the interpolation function a constant.We must then also consider the case where U is "close" to Ur or Ul, and an interpolating polynomial between them would not be monotonic over the cell. So, the basic algorithm is as follows:

Compute dU = Ur - Ul
Compute mU = 0.5*(Ur+Ul).
If the cell has an extremum:
1. Set Ur = U
2. Set Ul = U
If U is too close to Ul
1. Move Ul farther away
If U is too close to Ur
1. Move Ur farther away

More rigorous definitions of "Too Close" and "Farther Away" are derived from parabolas with vertices on the interfaces, as can be seen in the code below:

In [8]:

%%writefile -a $Ccodesdir/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

static inline void monotonize_NRPy(const REAL *U,REAL *Ur,REAL *Ul) {
  const REAL dU = *Ur - *Ul;
  const REAL mU = 0.5*(*Ur+*Ul);

  if ( ((*Ur)-(*U))*((*U)-(*Ul)) <= 0.0) {
    *Ur = *U;
    *Ul = *U;
    return;
  }
  if ( dU*((*U)-mU) > (1.0/6.0)*SQR(dU)) {
    *Ul = 3.0*(*U) - 2.0*(*Ur);
    return;
  }
  if ( dU*((*U)-mU) < -(1.0/6.0)*SQR(dU)) {
    *Ur = 3.0*(*U) - 2.0*(*Ul);
    return;
  }
}

Appending to GiRaFFE_standalone_Ccodes/PPM/reconstruct_set_of_prims_PPM_GRFFE_NRPy.c

Step 4: Dependencies [Back to top]¶

$$\label{loops}$$

This above functions include calls to several functions not given there. We will include those files below. These set the indices over which we will need to loop to gather the data needed to perform the reconstruction at some point i,j,k, as well as start such a loop.

In [9]:

%%writefile $Ccodesdir/loop_defines_reconstruction_NRPy.h
#ifndef loop_defines_reconstruction_NRPy_H_
#define loop_defines_reconstruction_NRPy_H_

#define LOOP_DEFINE(gz_shift_lo,gz_shift_hi,  ext,flux_dirn,  ijkgz_lo_hi,gz_lo,gz_hi) \
  for(int rr=1;rr<=3;rr++) {                                            \
    ijkgz_lo_hi[rr][0]=          gz_lo[rr];                             \
    ijkgz_lo_hi[rr][1]=ext[rr-1]-gz_hi[rr];                             \
  }                                                                     \
  ijkgz_lo_hi[flux_dirn][0] += gz_shift_lo;                             \
  ijkgz_lo_hi[flux_dirn][1] -= gz_shift_hi;                             \
  /* The following line is valid C99 */                                 \
  _Pragma("omp parallel for private(U,dU,slope_lim_dU,Ur,Ul)")          \
  for(int k=ijkgz_lo_hi[3][0];k<ijkgz_lo_hi[3][1];k++)                  \
    for(int j=ijkgz_lo_hi[2][0];j<ijkgz_lo_hi[2][1];j++)                \
      for(int i=ijkgz_lo_hi[1][0];i<ijkgz_lo_hi[1][1];i++)

// This define only sets indices.
// FIXME: benchmark with and without the if() statement.
// FIXME: try without index_arr being defined in all directions.
#define SET_INDEX_ARRAYS_NRPY(IMIN,IMAX,flux_dirn)                           \
  const int max_shift=(MAXNUMINDICES/2);                                \
  /* DEBUGGING ONLY:  if(IMIN<-max_shift || IMAX>max_shift) CCTK_VError(VERR_DEF_PARAMS,"FIX MAXNUMINDICES!"); */ \
  int index_arr[4][MAXNUMINDICES];                                      \
  for(int idx=IMIN;idx<=IMAX;idx++) {                                   \
    index_arr[flux_dirn][idx+max_shift]=                                \
      IDX3S(                                                            \
            i+idx*kronecker_delta[flux_dirn][0],                        \
            j+idx*kronecker_delta[flux_dirn][1],                        \
            k+idx*kronecker_delta[flux_dirn][2]);                       \
  }

#define SET_INDEX_ARRAYS_NRPY_3DBLOCK(IJKLOHI)                               \
  const int max_shift=(MAXNUMINDICES/2);                                \
  int index_arr_3DB[MAXNUMINDICES][MAXNUMINDICES][MAXNUMINDICES];       \
  for(int idx_k=IJKLOHI[4];idx_k<=IJKLOHI[5];idx_k++) for(int idx_j=IJKLOHI[2];idx_j<=IJKLOHI[3];idx_j++) for(int idx_i=IJKLOHI[0];idx_i<=IJKLOHI[1];idx_i++) { \
        index_arr_3DB[idx_k+max_shift][idx_j+max_shift][idx_i+max_shift]=IDX3S(i+idx_i,j+idx_j,k+idx_k); \
      }

#endif /* loop_defines_reconstruction_NRPy_H_ */

Overwriting GiRaFFE_standalone_Ccodes/PPM/loop_defines_reconstruction_NRPy.h

Step 5: Code Validation against `GiRaFFE_NRPy_PPM.py` [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 generated by GiRaFFE_NRPy_PPM.py

In [10]:

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

import GiRaFFE_NRPy.GiRaFFE_NRPy_PPM as PPM
PPM.GiRaFFE_NRPy_PPM(valdir)

import difflib
import sys

print("Printing difference between original C code and this code...")
# Open the files to compare
files = ["reconstruct_set_of_prims_PPM_GRFFE_NRPy.c",
         "loop_defines_reconstruction_NRPy.h"]

for file in files:
    print("Checking file " + file)
    with open(os.path.join(valdir,file)) as file1, open(os.path.join(Ccodesdir,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(Ccodesdir+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.")
            sys.exit(1)

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

Step 6: Example code in order to use this algorithm [Back to top]¶

$$\label{example}$$

Recall that the gridfunction and ghostzone structures (here rewritten in proper C) require a pointer gf

// Keeping track of ghostzones between routines is a nightmare, so
//   we instead attach ghostzone info to each gridfunction and set
//   the ghostzone information correctly within each routine.
typedef struct __gf_and_gz_struct__ {
  REAL *gf;
  int gz_lo[4],gz_hi[4];
} gf_and_gz_struct;

that can be found in NRPy+ as

#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )

So, the prerequisites that go into calling the PPM function will be as follows.

gf_and_gz_struct in_prims[NUM_RECONSTRUCT_GFS], out_prims_r[NUM_RECONSTRUCT_GFS], out_prims_l[NUM_RECONSTRUCT_GFS];
  int which_prims_to_reconstruct[NUM_RECONSTRUCT_GFS],num_prims_to_reconstruct;

  int ww=0;
  in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*VALENCIAVU0GF; 
    out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU0GF; 
    out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU0GF; 
    ww++;
  in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*VALENCIAVU1GF; 
    out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU1GF; 
    out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU1GF; 
    ww++;
  in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*VALENCIAVU2GF; 
    out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_RU2GF; 
    out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*VALENCIAV_LU2GF; 
    ww++;
  in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*BU0GF; 
    out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU0GF; 
    out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU0GF; 
    ww++;
  in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*BU1GF; 
    out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU1GF; 
    out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU1GF; 
    ww++;
  in_prims[ww].gf      = auxevol_gfs + Nxxp2NG012*BU2GF; 
    out_prims_r[ww].gf = auxevol_gfs + Nxxp2NG012*B_RU2GF; 
    out_prims_l[ww].gf = auxevol_gfs + Nxxp2NG012*B_LU2GF; 
    ww++;

  // Prims are defined AT ALL GRIDPOINTS, so we set the # of ghostzones to zero:
  for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { in_prims[i].gz_lo[j]=0; in_prims[i].gz_hi[j]=0; }
  // Left/right variables are not yet defined, yet we set the # of gz's to zero by default:
  for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { out_prims_r[i].gz_lo[j]=0; out_prims_r[i].gz_hi[j]=0; }
  for(int i=0;i<NUM_RECONSTRUCT_GFS;i++) for(int j=1;j<=3;j++) { out_prims_l[i].gz_lo[j]=0; out_prims_l[i].gz_hi[j]=0; }

  ww=0;
  which_prims_to_reconstruct[ww]=VALENCIAVU0; ww++;
  which_prims_to_reconstruct[ww]=VALENCIAVU1; ww++;
  which_prims_to_reconstruct[ww]=VALENCIAVU2; ww++;
  which_prims_to_reconstruct[ww]=BU0;         ww++;
  which_prims_to_reconstruct[ww]=BU1;         ww++;
  which_prims_to_reconstruct[ww]=BU2;         ww++;
  num_prims_to_reconstruct=ww;
  // This function is housed in the file: "reconstruct_set_of_prims_PPM_GRFFE_NRPy.c"
  reconstruct_set_of_prims_PPM_GRFFE_NRPy(params, auxevol_gfs, flux_dirn, num_prims_to_reconstruct,                                                          which_prims_to_reconstruct, in_prims, out_prims_r, out_prims_l, temporary);

The array temporary should point to unused gridfunction for the sake of efficiency; in the actual code, AevolParen works quite well. For the unit test, with a smaller grid, a dedicated array works fine:

REAL temporary[Nxxp2NG012];

This supposes that the following constants have been set:

const int VX=0,VY=1,VZ=2,BX=3,BY=4,BZ=5;
const int NUM_RECONSTRUCT_GFS = 6;
const int Nxxp2NG012 = Nxx_plus_2NGHOSTS0*Nxx_plus_2NGHOSTS1*Nxx_plus_2NGHOSTS2;

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

In [11]:

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

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