text:
# Start-to-Finish Example: Unit Testing `GiRaFFE_NRPy`: $\tilde{S}_i$ source term, $A_i$ gauge term, and $\psi^6 \Phi$

## Author: Patrick Nelson

## This module Validates the `Stilde_source_A_gauge_Phi_rhs` routine for `GiRaFFE`.

**Notebook Status:** Validated

**Validation Notes:** This module validates the routines for the source terms in [Tutorial-GRHD_Equations-Cartesian](Tutorial-GRHD_Equations-Cartesian.ipynb) and [Tutorial-GRFFE_Equations-Cartesian](Tutorial-GRFFE_Equations-Cartesian.ipynb).

### NRPy+ Source Code for this module:
* [GRHD/equations.py](../../edit/GRHD/equations.py) [\[**tutorial**\]](../Tutorial-GRHD_Equations-Cartesian.ipynb) Generates the source term for the Poynting flux evolution equation.
* [GRFFE/equations.py](../../edit/GRFFE/equations.py) [\[**tutorial**\]](../Tutorial-GRFFE_Equations-Cartesian.ipynb) Generates the source terms for the induction equation for the vector potential and the evolution equation for $\Phi$.

## Introduction:

This notebook validates the code that will add simpler terms to the RHSs of our evolved, conservative variables in `GiRaFFE_NRPy`. These terms require only basic finite-differencing for their derivatives, so the code to generate them is far simpler.

It is, in general, good coding practice to unit test functions individually to verify that they produce the expected and intended output. We will generate test data with arbitrarily-chosen analytic functions and calculate gridfunctions at the cell centers on a small numeric grid. We will then run the algorithms to compute the RHS terms numerically along with similar algorithms to compute them analytically, for comparison. We will then compare the data output by both routines to show convergence.

When this notebook is run, the difference between the approximate and exact right-hand sides will be output to text files that can be found in the same directory as this notebook. These will be read in in [Step 3](#convergence), and used there to confirm second order convergence of the algorithm.

# Table of Contents
$$\label{toc}$$

This notebook is organized as follows

1. [Step 1](#setup): Set up core functions and parameters for unit testing the algorithm
    1. [Step 1.a](#metric_expressions) Write expressions for the metric gridfunctions
    1. [Step 1.b](#mhd_expressions) Write expressions for the magnetohydrodynamic gridfunctions
    1. [Step 1.c](#ccodekernels) Generate C functions to calculate the gridfunctions
    1. [Step 1.d](#analytic_comparison) Generate C code to analytically calculate RHSs
    1. [Step 1.e](#functions_to_test) Generate C code to numerically calculate RHSs
    1. [Step 1.f](#free_parameters) Set free parameters in the code and generate other functions
1. [Step 2](#mainc): `Stilde_source_A_gauge_Phi_rhs_unit_test.c`: The Main C Code
    1. [Step 2.a](#compile_run): Compile and run the code
1. [Step 3](#convergence): Code validation: Verify that relative error in numerical solution converges to zero at the expected order
1. [Step 4](#latex_pdf_output): Output this notebook to $\LaTeX$-formatted PDF file

# Step 1: Set up core functions and parameters for unit testing the algorithm \[Back to [top](#toc)\]
$$\label{setup}$$

We'll start by appending the relevant paths to `sys.path` so that we can access sympy modules in other places. Then, we'll import NRPy+ core functionality and set up a directory in which to carry out our test.

```python
import shutil, os, sys           # Standard Python modules for multiplatform OS-level functions
# First, we'll add the parent directory to the list of directories Python will check for modules.
nrpy_dir_path = os.path.join("..")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)
nrpy_dir_path = os.path.join("..","..")
if nrpy_dir_path not in sys.path:
    sys.path.append(nrpy_dir_path)

from outputC import *            # 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 loop as lp                # NRPy+: Generate C code loops
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import reference_metric as rfm   # NRPy+: Reference metric support
import cmdline_helper as cmd     # NRPy+: Multi-platform Python command-line interface

out_dir = "Validation/"
cmd.mkdir(out_dir)
subdir = "source_terms"
cmd.mkdir(os.path.join(out_dir,subdir))

thismodule = "Start_to_Finish_UnitTest-GiRaFFE_NRPy-Source_Terms"

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

CoordSystem = "Cartesian"

par.set_parval_from_str("reference_metric::CoordSystem",CoordSystem)
rfm.reference_metric() # Create ReU, ReDD needed for rescaling B-L initial data, generating BSSN RHSs, etc.

# Default Kreiss-Oliger dissipation strength
default_KO_strength = 0.1
diss_strength = par.Cparameters("REAL", thismodule, "diss_strength", default_KO_strength)
```

## Step 1.a: Write expressions for the metric gridfunctions \[Back to [top](#toc)\]
$$\label{metric_expressions}$$

Now, we'll choose some functions with arbitrary forms to generate test data. We'll need to set ten gridfunctions, so expressions are being pulled from several previously written unit tests.

\begin{align}
\gamma_{xx} &= a \exp\left(-\left((x-b)^2+(y-c)^2+(z-d)^2\right)\right) \\
\gamma_{yy} &= f \exp\left(-\left((x-g)^2+(y-h)^2+(z-l)^2\right)\right) \\
\gamma_{zz} &= m \exp\left(-\left((x-n)^2+(y-o)^2+(z-p)^2\right)\right), \\
\gamma_{xy} &= ax^3 + by^3 + cz^3 + dy^2 + ez^2 + f \\
\gamma_{xz} &= gx^3 + hy^3 + lz^3 + mx^2 + nz^2 + p \\
\gamma_{yz} &= px^3 + qy^3 + rz^3 + sx^2 + ty^2 + u. \\
\beta^x &= \arctan(ax + by + cz) \\
\beta^y &= \arctan(bx + cy + az) \\
\beta^z &= \arctan(cx + ay + bz) \\
\alpha &= 1 - \frac{1}{2+x^2+y^2+z^2} \\
\end{align}

```python
a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t,u = par.Cparameters("REAL",thismodule,["a","b","c","d","e","f","g","h","l","m","n","o","p","q","r","s","t","u"],1e300)
M_PI  = par.Cparameters("#define",thismodule,["M_PI"], "")

gammaDD = ixp.register_gridfunctions_for_single_rank2("AUXEVOL","gammaDD","sym01",DIM=3)
betaU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","betaU",DIM=3)
alpha = gri.register_gridfunctions("AUXEVOL","alpha")

par.set_parval_from_str("reference_metric::CoordSystem","Cartesian")
rfm.reference_metric()
x = rfm.xxCart[0]
y = rfm.xxCart[1]
z = rfm.xxCart[2]

gammaDD[0][0] = 1.0 + 0.1*(a * sp.exp(-((x-b)**2 + (y-c)**2 + (z-d)**2)))
gammaDD[1][1] = 1.0 + 0.1*(f * sp.exp(-((x-g)**2 + (y-h)**2 + (z-l)**2)))
gammaDD[2][2] = 1.0 + 0.1*(m * sp.exp(-((x-n)**2 + (y-o)**2 + (z-p)**2)))
gammaDD[0][1] = gammaDD[1][0] = 0.1*(a*x**3 + b*y**3 + c*z**3 + d*y**2 + e*z**2 + f)
gammaDD[0][2] = gammaDD[2][0] = 0.1*(g*x**3 + h*y**3 + l*z**3 + m*x**2 + n*z**2 + o)
gammaDD[1][2] = gammaDD[2][1] = 0.1*(p*x**3 + q*y**3 + r*z**3 + s*x**2 + t*y**2 + u)

betaU[0] = sp.atan(a*x + b*y + c*z)
betaU[1] = sp.atan(b*x + c*y + a*z)
betaU[2] = sp.atan(c*x + a*y + b*z)

alpha = sp.sympify(1) - sp.sympify(1) / (sp.sympify(2) + x**2 + y**2 + z**2)
```

## Step 1.b: Write expressions for the magnetohydrodynamic gridfunctions \[Back to [top](#toc)\]
$$\label{mhd_expressions}$$

We will also pick arbitrary analytic forms for the vector potential and Valencia three-velocity here.
\begin{align}
A_x &= \exp(ey+fz) \\
A_y &= \exp(fz+dx) \\
A_z &= \exp(dx+ey) \\
\bar{v}^x &= \frac{2}{\pi} \arctan(ax + by + cz) \\
\bar{v}^y &= \frac{2}{\pi} \arctan(bx + cy + az) \\
\bar{v}^z &= \frac{2}{\pi} \arctan(cx + ay + bz) \\
\end{align}

We compute the magnetic field from the vector potential analytically as $B^i = \epsilon^{ijk} A_{k,j}$. Additionally, we let the scalar potential $\sqrt{\gamma} \Phi$ be defined as $\sqrt{\gamma} \left(1+x^2+y^2+z^2\right)$, where $\gamma$ is the determinant of the three metric $\gamma_{ij}$

```python
AD = ixp.register_gridfunctions_for_single_rank1("EVOL","AD")
BU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","BU")
ValenciavU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","ValenciavU")

gammaUU,gammadet = ixp.symm_matrix_inverter3x3(gammaDD)

AD[0] = sp.exp(e*y+f*z)
AD[1] = sp.exp(f*z+d*x)
AD[2] = sp.exp(d*x+e*y)
ValenciavU[0] = (sp.sympify(2)/M_PI) * sp.atan(d*x + e*y + f*z)
ValenciavU[1] = (sp.sympify(2)/M_PI) * sp.atan(e*x + f*y + d*z)
ValenciavU[2] = (sp.sympify(2)/M_PI) * sp.atan(f*x + d*y + e*z)

LeviCivitaUUU = ixp.LeviCivitaTensorUUU_dim3_rank3(sp.sqrt(gammadet))

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] * sp.diff(AD[k],rfm.xxCart[j])

psi6Phi = gri.register_gridfunctions("EVOL","psi6Phi")
psi6Phi = gammadet * (1+x*x+y*y+z*z)
```

## Step 1.c: Generate C functions to calculate the gridfunctions \[Back to [top](#toc)\]
$$\label{ccodekernels}$$

We now use NRPy+'s `outCfunction()` to generate the functions that will calculate our sample data on numerical grids.

```python
metric_gfs_to_print = [\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD00"),rhs=gammaDD[0][0]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD01"),rhs=gammaDD[0][1]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD02"),rhs=gammaDD[0][2]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD11"),rhs=gammaDD[1][1]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD12"),rhs=gammaDD[1][2]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","gammaDD22"),rhs=gammaDD[2][2]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","betaU0"),rhs=betaU[0]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","betaU1"),rhs=betaU[1]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","betaU2"),rhs=betaU[2]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","alpha"),rhs=alpha),\
                      ]
desc = "Calculate the metric gridfunctions"
name = "calculate_metric_gfs"
outCfunction(
    outfile  = os.path.join(out_dir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *restrict params,REAL *restrict xx[3],REAL *restrict auxevol_gfs",
    body     = fin.FD_outputC("returnstring",metric_gfs_to_print,params="outCverbose=False").replace("IDX4","IDX4S"),
    loopopts ="AllPoints,Read_xxs")

EM_to_print = [\
               lhrh(lhs=gri.gfaccess("auxevol_gfs","BU0"),rhs=BU[0]),\
               lhrh(lhs=gri.gfaccess("auxevol_gfs","BU1"),rhs=BU[1]),\
               lhrh(lhs=gri.gfaccess("auxevol_gfs","BU2"),rhs=BU[2]),\
               lhrh(lhs=gri.gfaccess("out_gfs","AD0"),rhs=AD[0]),\
               lhrh(lhs=gri.gfaccess("out_gfs","AD1"),rhs=AD[1]),\
               lhrh(lhs=gri.gfaccess("out_gfs","AD2"),rhs=AD[2]),\
               lhrh(lhs=gri.gfaccess("out_gfs","psi6Phi"),rhs=psi6Phi),\
              ]

desc = "Calculate sample magnetic field and potential data"
name = "calculate_BU_AD_psi6Phi"
outCfunction(
    outfile  = os.path.join(out_dir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs,REAL *out_gfs",
    body     = fin.FD_outputC("returnstring",EM_to_print,params="outCverbose=False").replace("IDX4","IDX4S"),
    loopopts ="AllPoints,Read_xxs")

ValenciavU_to_print = [\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU0"),rhs=ValenciavU[0]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU1"),rhs=ValenciavU[1]),\
                       lhrh(lhs=gri.gfaccess("auxevol_gfs","ValenciavU2"),rhs=ValenciavU[2]),\
                      ]

desc = "Calculate sample velocity data"
name = "calculate_ValenciavU"
outCfunction(
    outfile  = os.path.join(out_dir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],REAL *auxevol_gfs",
    body     = fin.FD_outputC("returnstring",ValenciavU_to_print,params="outCverbose=False").replace("IDX4","IDX4S"),
    loopopts ="AllPoints,Read_xxs")

```

    Output C function calculate_metric_gfs() to file Validation/calculate_metric_gfs.h
    Output C function calculate_BU_AD_psi6Phi() to file Validation/calculate_BU_AD_psi6Phi.h
    Output C function calculate_ValenciavU() to file Validation/calculate_ValenciavU.h

## Step 1.d: Generate C code to analytically calculate RHSs \[Back to [top](#toc)\]
$$\label{analytic_comparison}$$

In order to do a unit test, we have to have some expected value against which to compare. In this case, we'll calculate the expected output of our functions using the analytic forms of the gridfunction defined above. The first step to doing this will be to analytically compute the derivatives that we will need. For the $\partial_t \tilde{S}_i$ source term, this will just be the metric derivatives.

```python
gammaDD_dD = ixp.zerorank3()
for i in range(3):
    for j in range(3):
        for k in range(3):
            gammaDD_dD[i][j][k] = sp.diff(gammaDD[i][j],rfm.xxCart[k])

betaU_dD = ixp.zerorank2()
for i in range(3):
    for j in range(3):
        betaU_dD[i][j] = sp.diff(betaU[i],rfm.xxCart[j])

alpha_dD = ixp.zerorank1()
for i in range(3):
    alpha_dD[i] = sp.diff(alpha,rfm.xxCart[i])

import GRHD.equations as GRHD
import GRFFE.equations as GRFFE
TINYDOUBLE = par.Cparameters("REAL",thismodule,"TINYDOUBLE",1e-100)

GRHD.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha,betaU,gammaDD, ValenciavU)

# Next sqrt(gamma)
GRHD.compute_sqrtgammaDET(gammaDD)

# Then compute g4DD_zerotimederiv_dD
GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDD_dD,betaU_dD,alpha_dD)
# small b
sqrt4pi = par.Cparameters("REAL",thismodule,"sqrt4pi","sqrt(4.0*M_PI)")
GRFFE.compute_smallb4U(gammaDD,betaU,alpha, GRHD.u4U_ito_ValenciavU, BU, sqrt4pi)
GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)
# Electromagnetic stress-energy tensor
GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4U, GRFFE.smallbsquared,GRHD.u4U_ito_ValenciavU)

# Add the source term to the RHS
Stilde_rhsD = ixp.zerorank1(DIM=3)
for i in range(3):
    for mu in range(4):
        for nu in range(4):
            # \frac{1}{2} \alpha \sqrt{\gamma} T^{\mu \nu}_{\rm EM} \partial_i g_{\mu \nu}
            Stilde_rhsD[i] += sp.Rational(1,2) * alpha * GRHD.sqrtgammaDET * \
                              GRFFE.TEM4UU[mu][nu] * GRHD.g4DD_zerotimederiv_dD[mu][nu][i+1]

```

Next, we will handle the $\partial_t A_i$ gauge term, $-\partial_i (\alpha \Phi - \beta^j A_j).$ Carrying through this differentiation analytically, we see that we will need to code
$$
-\alpha_{,i} \Phi - \alpha \Phi_{,i} + \beta^j_{,i} A_j + \beta^j A_{j,i}
$$

```python
Phi_dD = ixp.zerorank1()
for i in range(3):
    Phi_dD[i] = sp.diff(psi6Phi/GRHD.sqrtgammaDET,rfm.xxCart[i])
AD_dD = ixp.zerorank2()
for i in range(3):
    for j in range(3):
        AD_dD[i][j] = sp.diff(AD[i],rfm.xxCart[j])

A_rhsD = ixp.zerorank1()
for i in range(3):
    A_rhsD[i] += -alpha_dD[i]*psi6Phi/GRHD.sqrtgammaDET - alpha*Phi_dD[i]

for i in range(3):
    for j in range(3):
        A_rhsD[i] += betaU_dD[j][i]*AD[j] + betaU[j]*AD_dD[j][i]

```

And finally, we will code the equation
\begin{align}
\partial_t [\sqrt{\gamma} \Phi] &= -\partial_j (\alpha\sqrt{\gamma}A^j - \beta^j [\sqrt{\gamma} \Phi]) - \xi \alpha [\sqrt{\gamma} \Phi] \\
&= -\alpha_{,j} \sqrt{\gamma} A^j - \alpha (\sqrt{\gamma})_{,j} A^j - \alpha \sqrt{\gamma} A^j_{,j} + \beta^j_{,j} [\sqrt{\gamma} \Phi] + \beta^j [\sqrt{\gamma} \Phi]_{,j} - \xi \alpha [\sqrt{\gamma} \Phi]
\end{align}

```python
psi6Phi_dD = ixp.zerorank1()
for i in range(3):
    psi6Phi_dD[i] = sp.diff(psi6Phi,rfm.xxCart[i])

gammaUU,unusedgammadet = ixp.symm_matrix_inverter3x3(gammaDD)
AU = ixp.zerorank1()
for i in range(3):
    for j in range(3):
        AU[i] += gammaUU[i][j]*AD[j]
AU_dD = ixp.zerorank2()
for i in range(3):
    for j in range(3):
        AU_dD[i][j] = sp.diff(AU[i],rfm.xxCart[j])

xi_damping = par.Cparameters("REAL",thismodule,"xi_damping",0.1)
psi6Phi_rhs = -xi_damping*alpha*psi6Phi
for j in range(3):
    psi6Phi_rhs += -alpha_dD[j]*GRHD.sqrtgammaDET*AU[j] \
                   -alpha*sp.diff(GRHD.sqrtgammaDET,rfm.xxCart[j])*AU[j] \
                   -alpha*GRHD.sqrtgammaDET*AU_dD[j][j] \
                   + betaU_dD[j][j]*psi6Phi + betaU[j]*psi6Phi_dD[j]
```

With these expressions coded, we can finally generate the C function to write this data for a comparison. We will do this using `outCfunction()`.

```python
StildeD = ixp.register_gridfunctions_for_single_rank1("EVOL","StildeD")

RHSs_to_print = [\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","StildeD0"),rhs=Stilde_rhsD[0]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","StildeD1"),rhs=Stilde_rhsD[1]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","StildeD2"),rhs=Stilde_rhsD[2]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","AD0"),rhs=A_rhsD[0]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","AD1"),rhs=A_rhsD[1]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","AD2"),rhs=A_rhsD[2]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","psi6Phi"),rhs=psi6Phi_rhs),\
                ]

desc = "Calculate analytic RHSs for comparison"
name = "calculate_analytic_RHSs"
outCfunction(
    outfile  = os.path.join(out_dir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],REAL *rhs_gfs",
    body     = fin.FD_outputC("returnstring",RHSs_to_print,params="outCverbose=False").replace("IDX4","IDX4S"),
    loopopts ="AllPoints,Read_xxs")

```

    Output C function calculate_analytic_RHSs() to file Validation/calculate_analytic_RHSs.h

## Step 1.e: Generate C code to numerically calculate RHSs \[Back to [top](#toc)\]
$$\label{functions_to_test}$$

Here, we will generate the functions that we specifically wish to test. We'll start with the source term for $\partial_t \tilde{S}_i$. This term involves derivatives of the metric tensor, and in `GiRaFFE`, we will be storing interpolations of the metric quantities to the cell faces of our grid. This will allow us to save some time in our simulations if we compute the finite differences in a somewhat unusual way. If we do a simple, first-order finite-difference derivative of these interpolations, we can take a derivative with the accuracy of a higher-order stencil but the computational cost of a first-order stencil. To do this, we will create arrays of NRPy+'s `Cparameter`s to pass to the function that will build the RHS, then write the lines of code by hand to compute the derivatives for them.

Additionally, we will write this as three similar functions, corresponding to one for each spatial direction in which we take derivatives.

```python
# Reset gridfunctions to basenames:
gammaDD = ixp.declarerank2("gammaDD","sym01",DIM=3)
betaU = ixp.declarerank1("betaU",DIM=3)
alpha = sp.symbols("alpha",real=True)
ValenciavU = ixp.declarerank1("ValenciavU",DIM=3)
BU = ixp.declarerank1("BU",DIM=3)
AD = ixp.declarerank1("AD",DIM=3)
psi6Phi = sp.symbols("psi6Phi",real=True)

# Declare all the Cparameterse we will need
metricderivDDD = ixp.declarerank3("metricderivDDD","sym01",DIM=3)
shiftderivUD = ixp.declarerank2("shiftderivUD","nosym",DIM=3)
lapsederivD = ixp.declarerank1("lapsederivD",DIM=3)

general_access = """const REAL gammaDD00 = auxevol_gfs[IDX4S(GAMMADD00GF,i0,i1,i2)];
const REAL gammaDD01 = auxevol_gfs[IDX4S(GAMMADD01GF,i0,i1,i2)];
const REAL gammaDD02 = auxevol_gfs[IDX4S(GAMMADD02GF,i0,i1,i2)];
const REAL gammaDD11 = auxevol_gfs[IDX4S(GAMMADD11GF,i0,i1,i2)];
const REAL gammaDD12 = auxevol_gfs[IDX4S(GAMMADD12GF,i0,i1,i2)];
const REAL gammaDD22 = auxevol_gfs[IDX4S(GAMMADD22GF,i0,i1,i2)];
const REAL betaU0 = auxevol_gfs[IDX4S(BETAU0GF,i0,i1,i2)];
const REAL betaU1 = auxevol_gfs[IDX4S(BETAU1GF,i0,i1,i2)];
const REAL betaU2 = auxevol_gfs[IDX4S(BETAU2GF,i0,i1,i2)];
const REAL alpha = auxevol_gfs[IDX4S(ALPHAGF,i0,i1,i2)];
const REAL ValenciavU0 = auxevol_gfs[IDX4S(VALENCIAVU0GF,i0,i1,i2)];
const REAL ValenciavU1 = auxevol_gfs[IDX4S(VALENCIAVU1GF,i0,i1,i2)];
const REAL ValenciavU2 = auxevol_gfs[IDX4S(VALENCIAVU2GF,i0,i1,i2)];
const REAL BU0 = auxevol_gfs[IDX4S(BU0GF,i0,i1,i2)];
const REAL BU1 = auxevol_gfs[IDX4S(BU1GF,i0,i1,i2)];
const REAL BU2 = auxevol_gfs[IDX4S(BU2GF,i0,i1,i2)];
"""
metric_deriv_access = ixp.zerorank1(DIM=3)
metric_deriv_access[0] = """const REAL metricderivDDD000 = (auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1,i2)])/dxx0;
const REAL metricderivDDD010 = (auxevol_gfs[IDX4S(GAMMA_FACEDD01GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD01GF,i0,i1,i2)])/dxx0;
const REAL metricderivDDD020 = (auxevol_gfs[IDX4S(GAMMA_FACEDD02GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD02GF,i0,i1,i2)])/dxx0;
const REAL metricderivDDD110 = (auxevol_gfs[IDX4S(GAMMA_FACEDD11GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD11GF,i0,i1,i2)])/dxx0;
const REAL metricderivDDD120 = (auxevol_gfs[IDX4S(GAMMA_FACEDD12GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD12GF,i0,i1,i2)])/dxx0;
const REAL metricderivDDD220 = (auxevol_gfs[IDX4S(GAMMA_FACEDD22GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD22GF,i0,i1,i2)])/dxx0;
const REAL shiftderivUD00 = (auxevol_gfs[IDX4S(BETA_FACEU0GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(BETA_FACEU0GF,i0,i1,i2)])/dxx0;
const REAL shiftderivUD10 = (auxevol_gfs[IDX4S(BETA_FACEU1GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(BETA_FACEU1GF,i0,i1,i2)])/dxx0;
const REAL shiftderivUD20 = (auxevol_gfs[IDX4S(BETA_FACEU2GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(BETA_FACEU2GF,i0,i1,i2)])/dxx0;
const REAL lapsederivD0 = (auxevol_gfs[IDX4S(ALPHA_FACEGF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(ALPHA_FACEGF,i0,i1,i2)])/dxx0;
REAL Stilde_rhsD0;
"""
metric_deriv_access[1] = """const REAL metricderivDDD001 = (auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1+1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1,i2)])/dxx1;
const REAL metricderivDDD011 = (auxevol_gfs[IDX4S(GAMMA_FACEDD01GF,i0,i1+1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD01GF,i0,i1,i2)])/dxx1;
const REAL metricderivDDD021 = (auxevol_gfs[IDX4S(GAMMA_FACEDD02GF,i0,i1+1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD02GF,i0,i1,i2)])/dxx1;
const REAL metricderivDDD111 = (auxevol_gfs[IDX4S(GAMMA_FACEDD11GF,i0,i1+1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD11GF,i0,i1,i2)])/dxx1;
const REAL metricderivDDD121 = (auxevol_gfs[IDX4S(GAMMA_FACEDD12GF,i0,i1+1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD12GF,i0,i1,i2)])/dxx1;
const REAL metricderivDDD221 = (auxevol_gfs[IDX4S(GAMMA_FACEDD22GF,i0,i1+1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD22GF,i0,i1,i2)])/dxx1;
const REAL shiftderivUD01 = (auxevol_gfs[IDX4S(BETA_FACEU0GF,i0,i1+1,i2)]-auxevol_gfs[IDX4S(BETA_FACEU0GF,i0,i1,i2)])/dxx1;
const REAL shiftderivUD11 = (auxevol_gfs[IDX4S(BETA_FACEU1GF,i0,i1+1,i2)]-auxevol_gfs[IDX4S(BETA_FACEU1GF,i0,i1,i2)])/dxx1;
const REAL shiftderivUD21 = (auxevol_gfs[IDX4S(BETA_FACEU2GF,i0,i1+1,i2)]-auxevol_gfs[IDX4S(BETA_FACEU2GF,i0,i1,i2)])/dxx1;
const REAL lapsederivD1 = (auxevol_gfs[IDX4S(ALPHA_FACEGF,i0,i1+1,i2)]-auxevol_gfs[IDX4S(ALPHA_FACEGF,i0,i1,i2)])/dxx1;
REAL Stilde_rhsD1;
"""
metric_deriv_access[2] = """const REAL metricderivDDD002 = (auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1,i2+1)]-auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1,i2)])/dxx2;
const REAL metricderivDDD012 = (auxevol_gfs[IDX4S(GAMMA_FACEDD01GF,i0,i1,i2+1)]-auxevol_gfs[IDX4S(GAMMA_FACEDD01GF,i0,i1,i2)])/dxx2;
const REAL metricderivDDD022 = (auxevol_gfs[IDX4S(GAMMA_FACEDD02GF,i0,i1,i2+1)]-auxevol_gfs[IDX4S(GAMMA_FACEDD02GF,i0,i1,i2)])/dxx2;
const REAL metricderivDDD112 = (auxevol_gfs[IDX4S(GAMMA_FACEDD11GF,i0,i1,i2+1)]-auxevol_gfs[IDX4S(GAMMA_FACEDD11GF,i0,i1,i2)])/dxx2;
const REAL metricderivDDD122 = (auxevol_gfs[IDX4S(GAMMA_FACEDD12GF,i0,i1,i2+1)]-auxevol_gfs[IDX4S(GAMMA_FACEDD12GF,i0,i1,i2)])/dxx2;
const REAL metricderivDDD222 = (auxevol_gfs[IDX4S(GAMMA_FACEDD22GF,i0,i1,i2+1)]-auxevol_gfs[IDX4S(GAMMA_FACEDD22GF,i0,i1,i2)])/dxx2;
const REAL shiftderivUD02 = (auxevol_gfs[IDX4S(BETA_FACEU0GF,i0,i1,i2+1)]-auxevol_gfs[IDX4S(BETA_FACEU0GF,i0,i1,i2)])/dxx2;
const REAL shiftderivUD12 = (auxevol_gfs[IDX4S(BETA_FACEU1GF,i0,i1,i2+1)]-auxevol_gfs[IDX4S(BETA_FACEU1GF,i0,i1,i2)])/dxx2;
const REAL shiftderivUD22 = (auxevol_gfs[IDX4S(BETA_FACEU2GF,i0,i1,i2+1)]-auxevol_gfs[IDX4S(BETA_FACEU2GF,i0,i1,i2)])/dxx2;
const REAL lapsederivD2 = (auxevol_gfs[IDX4S(ALPHA_FACEGF,i0,i1,i2+1)]-auxevol_gfs[IDX4S(ALPHA_FACEGF,i0,i1,i2)])/dxx2;
REAL Stilde_rhsD2;
"""
write_final_quantity = ixp.zerorank1(DIM=3)
write_final_quantity[0] = """rhs_gfs[IDX4S(STILDED0GF,i0,i1,i2)] += Stilde_rhsD0;
"""
write_final_quantity[1] = """rhs_gfs[IDX4S(STILDED1GF,i0,i1,i2)] += Stilde_rhsD1;
"""
write_final_quantity[2] = """rhs_gfs[IDX4S(STILDED2GF,i0,i1,i2)] += Stilde_rhsD2;
"""
# We need to rerun a few of these functions with the reset lists to make sure these functions
# don't cheat by using analytic expressions
GRHD.compute_sqrtgammaDET(gammaDD)
GRHD.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha, betaU, gammaDD, ValenciavU)
GRFFE.compute_smallb4U(gammaDD, betaU, alpha, GRHD.u4U_ito_ValenciavU, BU, sqrt4pi)
GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)
GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4U, GRFFE.smallbsquared,GRHD.u4U_ito_ValenciavU)
GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, metricderivDDD,shiftderivUD,lapsederivD)
GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD, GRFFE.TEM4UU)
for i in range(3):
    desc = "Adds the source term to StildeD"+str(i)+"."
    name = "calculate_StildeD"+str(i)+"_source_term"
    outCfunction(
        outfile  = os.path.join(out_dir,subdir,name+".h"), desc=desc, name=name,
        params   ="const paramstruct *params,REAL *xx[3],const REAL *auxevol_gfs, REAL *rhs_gfs",
        body     = general_access \
                  +metric_deriv_access[i]\
                  +outputC(GRHD.S_tilde_source_termD[i],"Stilde_rhsD"+str(i),"returnstring",params="outCverbose=False").replace("IDX4","IDX4S")\
                  +write_final_quantity[i],
        loopopts ="InteriorPoints",
        rel_path_for_Cparams=os.path.join("../"))

```

    Output C function calculate_StildeD0_source_term() to file Validation/source_terms/calculate_StildeD0_source_term.h
    Output C function calculate_StildeD1_source_term() to file Validation/source_terms/calculate_StildeD1_source_term.h
    Output C function calculate_StildeD2_source_term() to file Validation/source_terms/calculate_StildeD2_source_term.h

Next, we will write the C function to compute the the quantities inside the parentheses for the $A_i$ gauge term and $\sqrt{\gamma} \Phi$ RHS. This needs to be done as a separate step before we compute the contributions to the RHSs themselves in order to fully leverage the automated finite-differencing for these terms.

First, we will also want to add Kreiss-Oliger dissipation in order to reduce Gibbs' phenomenon near shocks.

```python
PhievolParenU = ixp.register_gridfunctions_for_single_rank1("AUXEVOL","PhievolParenU",DIM=3)
AevolParen = gri.register_gridfunctions("AUXEVOL","AevolParen")

GRFFE.compute_AD_source_term_parenthetical_for_FD(GRHD.sqrtgammaDET,betaU,alpha,psi6Phi,AD)
GRFFE.compute_psi6Phi_rhs_parenthetical(gammaDD,GRHD.sqrtgammaDET,betaU,alpha,AD,psi6Phi)

parens_to_print = [\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","AevolParen"),rhs=GRFFE.AevolParen),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","PhievolParenU0"),rhs=GRFFE.PhievolParenU[0]),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","PhievolParenU1"),rhs=GRFFE.PhievolParenU[1]),\
                   lhrh(lhs=gri.gfaccess("auxevol_gfs","PhievolParenU2"),rhs=GRFFE.PhievolParenU[2]),\
                  ]

desc = "Calculate quantities to be finite-differenced for the GRFFE RHSs"
name = "calculate_parentheticals_for_RHSs"
outCfunction(
    outfile  = os.path.join(out_dir,subdir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],const REAL *in_gfs,REAL *auxevol_gfs",
    body     = fin.FD_outputC("returnstring",parens_to_print,params="outCverbose=False").replace("IDX4","IDX4S"),
    loopopts ="AllPoints",
    rel_path_for_Cparams=os.path.join("../"))

```

    Output C function calculate_parentheticals_for_RHSs() to file Validation/source_terms/calculate_parentheticals_for_RHSs.h

With these terms computed, we can move on to generating the function that will add to the RHSs.

```python
GRFFE.compute_psi6Phi_rhs_damping_term(alpha,psi6Phi,xi_damping)

AevolParen_dD = ixp.declarerank1("AevolParen_dD",DIM=3)
PhievolParenU_dD = ixp.declarerank2("PhievolParenU_dD","nosym",DIM=3)

A_rhsD = ixp.zerorank1()
psi6Phi_rhs = GRFFE.psi6Phi_damping

for i in range(3):
    A_rhsD[i] += -AevolParen_dD[i]
    psi6Phi_rhs += -PhievolParenU_dD[i][i]

# Add Kreiss-Oliger dissipation to the GRFFE RHSs:

psi6Phi_dKOD = ixp.declarerank1("psi6Phi_dKOD")
AD_dKOD    = ixp.declarerank2("AD_dKOD","nosym")
for i in range(3):
    psi6Phi_rhs += diss_strength*psi6Phi_dKOD[i]*rfm.ReU[i] # ReU[i] = 1/scalefactor_orthog_funcform[i]
    for j in range(3):
        A_rhsD[j] += diss_strength*AD_dKOD[j][i]*rfm.ReU[i] # ReU[i] = 1/scalefactor_orthog_funcform[i]

RHSs_to_print = [\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","AD0"),rhs=A_rhsD[0]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","AD1"),rhs=A_rhsD[1]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","AD2"),rhs=A_rhsD[2]),\
                 lhrh(lhs=gri.gfaccess("rhs_gfs","psi6Phi"),rhs=psi6Phi_rhs),\
                ]

desc = "Calculate AD gauge term and psi6Phi RHSs"
name = "calculate_AD_gauge_psi6Phi_RHSs"
outCfunction(
    outfile  = os.path.join(out_dir,subdir,name+".h"), desc=desc, name=name,
    params   ="const paramstruct *params,REAL *xx[3],const REAL *in_gfs,const REAL *auxevol_gfs,REAL *rhs_gfs",
    body     = fin.FD_outputC("returnstring",RHSs_to_print,params="outCverbose=False").replace("IDX4","IDX4S"),
    loopopts ="InteriorPoints",
    rel_path_for_Cparams=os.path.join("../"))

```

    Output C function calculate_AD_gauge_psi6Phi_RHSs() to file Validation/source_terms/calculate_AD_gauge_psi6Phi_RHSs.h

## Step 1.f: Set free parameters in the code and generate other functions \[Back to [top](#toc)\]
$$\label{free_parameters}$$

We also need to create the files that interact with NRPy's C parameter interface.

```python
# Step 3.d.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
# par.generate_Cparameters_Ccodes(os.path.join(out_dir))

# Step 3.d.ii: Set free_parameters.h
with open(os.path.join(out_dir,"free_parameters.h"),"w") as file:
    file.write("""
// Override parameter defaults with values based on command line arguments and NGHOSTS.
params.Nxx0 = atoi(argv[1]);
params.Nxx1 = atoi(argv[2]);
params.Nxx2 = atoi(argv[3]);
params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*NGHOSTS;
params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*NGHOSTS;
params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*NGHOSTS;
// Step 0d: Set up space and time coordinates
// Step 0d.i: Declare \Delta x^i=dxx{0,1,2} and invdxx{0,1,2}, as well as xxmin[3] and xxmax[3]:
const REAL xxmin[3] = {-0.2,-0.2,-0.2};
const REAL xxmax[3] = { 0.2, 0.2, 0.2};

params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx0);
params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx1);
params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx2);
printf("dxx0,dxx1,dxx2 = %.5e,%.5e,%.5e\\n",params.dxx0,params.dxx1,params.dxx2);
params.invdx0 = 1.0 / params.dxx0;
params.invdx1 = 1.0 / params.dxx1;
params.invdx2 = 1.0 / params.dxx2;
\n""")

# Generates declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h
par.generate_Cparameters_Ccodes(os.path.join(out_dir))
```

Also, we need to run the function that will create C code to calculate the metric face values for the $\tilde{S}_i source term.

```python
import GiRaFFE_NRPy.GiRaFFE_NRPy_Metric_Face_Values as FCVAL
FCVAL.GiRaFFE_NRPy_FCVAL(os.path.join(out_dir,subdir))
```

# Step 2: `Stilde_source_A_gauge_Phi_rhs_unit_test.c`: The Main C Code \[Back to [top](#toc)\]
$$\label{mainc}$$

```python
%%writefile $out_dir/Stilde_source_A_gauge_Phi_rhs_unit_test.c
// These are common packages that we are likely to need.
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "string.h" // Needed for strncmp, etc.
#include "stdint.h" // Needed for Windows GCC 6.x compatibility
#include <time.h>   // Needed to set a random seed.

#define REAL double
#include "declare_Cparameters_struct.h"

const int NGHOSTS = 3;

REAL a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t,u;

// Standard NRPy+ memory access:
#define IDX4S(g,i,j,k) \
( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )

const int kronecker_delta[4][3] = { { 0,0,0 },
                                    { 1,0,0 },
                                    { 0,1,0 },
                                    { 0,0,1 } };

// Give gridfunctions their names:
#define GAMMADD00GF 0
#define GAMMADD01GF 1
#define GAMMADD02GF 2
#define GAMMADD11GF 3
#define GAMMADD12GF 4
#define GAMMADD22GF 5
#define BETAU0GF 6
#define BETAU1GF 7
#define BETAU2GF 8
#define ALPHAGF 9
#define GAMMA_FACEDD00GF 10
#define GAMMA_FACEDD01GF 11
#define GAMMA_FACEDD02GF 12
#define GAMMA_FACEDD11GF 13
#define GAMMA_FACEDD12GF 14
#define GAMMA_FACEDD22GF 15
#define BETA_FACEU0GF 16
#define BETA_FACEU1GF 17
#define BETA_FACEU2GF 18
#define ALPHA_FACEGF 19
#define VALENCIAVU0GF 20
#define VALENCIAVU1GF 21
#define VALENCIAVU2GF 22
#define BU0GF 23
#define BU1GF 24
#define BU2GF 25
#define AEVOLPARENGF 26
#define PHIEVOLPARENU0GF 27
#define PHIEVOLPARENU1GF 28
#define PHIEVOLPARENU2GF 29
#define NUM_AUXEVOL_GFS 30

#define AD0GF 0
#define AD1GF 1
#define AD2GF 2
#define STILDED0GF 3
#define STILDED1GF 4
#define STILDED2GF 5
#define PSI6PHIGF 6
#define NUM_EVOL_GFS 7

#include "calculate_metric_gfs.h"
#include "calculate_BU_AD_psi6Phi.h"
#include "calculate_ValenciavU.h"
#include "calculate_analytic_RHSs.h"
#include "source_terms/interpolate_metric_gfs_to_cell_faces.h"
#include "source_terms/calculate_StildeD0_source_term.h"
#include "source_terms/calculate_StildeD1_source_term.h"
#include "source_terms/calculate_StildeD2_source_term.h"
#include "source_terms/calculate_parentheticals_for_RHSs.h"
#include "source_terms/calculate_AD_gauge_psi6Phi_RHSs.h"

int main(int argc, const char *argv[]) {
    paramstruct params;
#include "set_Cparameters_default.h"

    // Step 0c: Set free parameters, overwriting Cparameters defaults
    //          by hand or with command-line input, as desired.
#include "free_parameters.h"
#include "set_Cparameters-nopointer.h"

    // Step 0e: Set up cell-centered Cartesian coordinate grids
    REAL *xx[3];
    xx[0] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS0);
    xx[1] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS1);
    xx[2] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS2);
    for(int j=0;j<Nxx_plus_2NGHOSTS0;j++) xx[0][j] = xxmin[0] + (j-NGHOSTS)*dxx0;
    for(int j=0;j<Nxx_plus_2NGHOSTS1;j++) xx[1][j] = xxmin[1] + (j-NGHOSTS)*dxx1;
    for(int j=0;j<Nxx_plus_2NGHOSTS2;j++) xx[2][j] = xxmin[2] + (j-NGHOSTS)*dxx2;

    for(int i=0;i<Nxx_plus_2NGHOSTS0;i++) printf("xx[0][%d] = %.15e\n",i,xx[0][i]);

    // This is the array to which we'll write the NRPy+ variables.
    REAL *auxevol_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_AUXEVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);
    REAL *evol_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);
    REAL *rhs_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);
    // And another for exact data:
    REAL *rhs_exact_gfs  = (REAL *)malloc(sizeof(REAL) * NUM_EVOL_GFS * Nxx_plus_2NGHOSTS2 * Nxx_plus_2NGHOSTS1 * Nxx_plus_2NGHOSTS0);

    // Generate some random coefficients. Leave the random seed on its default for consistency between trials.
    a = (double)(rand()%20)/5.0;
    f = (double)(rand()%20)/5.0;
    m = (double)(rand()%20)/5.0;
    b = (double)(rand()%10-5)/100.0;
    c = (double)(rand()%10-5)/100.0;
    d = (double)(rand()%10-5)/100.0;
    g = (double)(rand()%10-5)/100.0;
    h = (double)(rand()%10-5)/100.0;
    l = (double)(rand()%10-5)/100.0;
    n = (double)(rand()%10-5)/100.0;
    o = (double)(rand()%10-5)/100.0;
    p = (double)(rand()%10-5)/100.0;
    q = (double)(rand()%10-5)/100.0;
    r = (double)(rand()%10-5)/100.0;
    s = (double)(rand()%10-5)/100.0;
    t = (double)(rand()%10-5)/100.0;
    u = (double)(rand()%10-5)/100.0;

    // First, calculate the test data on our grid:
    calculate_metric_gfs(&params,xx,auxevol_gfs);
    calculate_BU_AD_psi6Phi(&params,xx,auxevol_gfs,evol_gfs);
    calculate_ValenciavU(&params,xx,auxevol_gfs);
    calculate_analytic_RHSs(&params,xx,rhs_exact_gfs);
    // Now, run our function:
    int flux_dirn = 0;
    interpolate_metric_gfs_to_cell_faces(&params,auxevol_gfs,flux_dirn+1);
    calculate_StildeD0_source_term(&params,xx,auxevol_gfs,rhs_gfs);
    flux_dirn = 1;
    interpolate_metric_gfs_to_cell_faces(&params,auxevol_gfs,flux_dirn+1);
    calculate_StildeD1_source_term(&params,xx,auxevol_gfs,rhs_gfs);
    flux_dirn = 2;
    interpolate_metric_gfs_to_cell_faces(&params,auxevol_gfs,flux_dirn+1);
    calculate_StildeD2_source_term(&params,xx,auxevol_gfs,rhs_gfs);
    calculate_parentheticals_for_RHSs(&params,xx,evol_gfs,auxevol_gfs);
    calculate_AD_gauge_psi6Phi_RHSs(&params,xx,evol_gfs,auxevol_gfs,rhs_gfs);
    /*printf("Exact: Stilde_rhsD: %.15e,%.15e,%.15e\n",rhs_exact_gfs[IDX4S(STILDED0GF,4,4,4)],rhs_exact_gfs[IDX4S(STILDED1GF,4,4,4)],rhs_exact_gfs[IDX4S(STILDED2GF,4,4,4)]);
    printf("Numer: Stilde_rhsD: %.15e,%.15e,%.15e\n",rhs_gfs[IDX4S(STILDED0GF,4,4,4)],rhs_gfs[IDX4S(STILDED1GF,4,4,4)],rhs_gfs[IDX4S(STILDED2GF,4,4,4)]);
    printf("Exact: A_rhsD: %.15e,%.15e,%.15e\n",rhs_exact_gfs[IDX4S(AD0GF,4,4,4)],rhs_exact_gfs[IDX4S(AD1GF,4,4,4)],rhs_exact_gfs[IDX4S(AD2GF,4,4,4)]);
    printf("Numer: A_rhsD: %.15e,%.15e,%.15e\n",rhs_gfs[IDX4S(AD0GF,4,4,4)],rhs_gfs[IDX4S(AD1GF,4,4,4)],rhs_gfs[IDX4S(AD2GF,4,4,4)]);
    printf("Exact: psi6Phi_rhs: %.15e\n",rhs_exact_gfs[IDX4S(PSI6PHIGF,4,4,4)]);
    printf("Numer: psi6Phi_rhs: %.15e\n",rhs_gfs[IDX4S(PSI6PHIGF,4,4,4)]);*/

    char filename[100];
    sprintf(filename,"out%d-numer.txt",Nxx0);
    FILE *out2D = fopen(filename, "w");
    // We print the difference between approximate and exact numbers.
    int i0 = Nxx_plus_2NGHOSTS0/2;
    int i1 = Nxx_plus_2NGHOSTS1/2;
    int i2 = Nxx_plus_2NGHOSTS2/2;
    fprintf(out2D,"%.16e\t%.16e\t%.16e\t%.16e\t%.16e\t%.16e\t%.16e\t %e %e %e\n",
            rhs_exact_gfs[IDX4S(STILDED0GF,i0,i1,i2)]-rhs_gfs[IDX4S(STILDED0GF,i0,i1,i2)],
            rhs_exact_gfs[IDX4S(STILDED1GF,i0,i1,i2)]-rhs_gfs[IDX4S(STILDED1GF,i0,i1,i2)],
            rhs_exact_gfs[IDX4S(STILDED2GF,i0,i1,i2)]-rhs_gfs[IDX4S(STILDED2GF,i0,i1,i2)],
            rhs_exact_gfs[IDX4S(AD0GF,i0,i1,i2)]-rhs_gfs[IDX4S(AD0GF,i0,i1,i2)],
            rhs_exact_gfs[IDX4S(AD1GF,i0,i1,i2)]-rhs_gfs[IDX4S(AD1GF,i0,i1,i2)],
            rhs_exact_gfs[IDX4S(AD2GF,i0,i1,i2)]-rhs_gfs[IDX4S(AD2GF,i0,i1,i2)],
            rhs_exact_gfs[IDX4S(PSI6PHIGF,i0,i1,i2)]-rhs_gfs[IDX4S(PSI6PHIGF,i0,i1,i2)],
            xx[0][i0],xx[1][i1],xx[2][i2]
            );
    fclose(out2D);
}
```

    Overwriting Validation//Stilde_source_A_gauge_Phi_rhs_unit_test.c

## Step 2.a: Compile and run the code \[Back to [top](#toc)\]
$$\label{compile_run}$$

Now that we have our file, we can compile it and run the executable.

```python
import time

print("Now compiling, should take ~2 seconds...\n")
start = time.time()
cmd.C_compile(os.path.join(out_dir,"Stilde_source_A_gauge_Phi_rhs_unit_test.c"), os.path.join(out_dir,"Stilde_source_A_gauge_Phi_rhs_unit_test"))
end = time.time()
print("Finished in "+str(end-start)+" seconds.\n\n")

print("Now running...\n")
start = time.time()
!./Validation/Stilde_source_A_gauge_Phi_rhs_unit_test 2 2 2
# To do a convergence test, we'll also need a second grid with twice the resolution.
!./Validation/Stilde_source_A_gauge_Phi_rhs_unit_test 4 4 4
end = time.time()
print("Finished in "+str(end-start)+" seconds.\n\n")

```

    Now compiling, should take ~2 seconds...

    Compiling executable...
    Executing `gcc -Ofast -fopenmp -march=native -funroll-loops Validation/Stilde_source_A_gauge_Phi_rhs_unit_test.c -o Validation/Stilde_source_A_gauge_Phi_rhs_unit_test -lm`...
    Finished executing in 1.2167835235595703 seconds.
    Finished compilation.
    Finished in 1.2294585704803467 seconds.

    Now running...

    dxx0,dxx1,dxx2 = 2.00000e-01,2.00000e-01,2.00000e-01
    xx[0][0] = -8.000000000000000e-01
    xx[0][1] = -6.000000000000001e-01
    xx[0][2] = -4.000000000000000e-01
    xx[0][3] = -2.000000000000000e-01
    xx[0][4] = 0.000000000000000e+00
    xx[0][5] = 2.000000000000000e-01
    xx[0][6] = 4.000000000000000e-01
    xx[0][7] = 6.000000000000001e-01
    dxx0,dxx1,dxx2 = 1.00000e-01,1.00000e-01,1.00000e-01
    xx[0][0] = -5.000000000000000e-01
    xx[0][1] = -4.000000000000000e-01
    xx[0][2] = -3.000000000000000e-01
    xx[0][3] = -2.000000000000000e-01
    xx[0][4] = -1.000000000000000e-01
    xx[0][5] = 0.000000000000000e+00
    xx[0][6] = 1.000000000000000e-01
    xx[0][7] = 2.000000000000000e-01
    xx[0][8] = 3.000000000000000e-01
    xx[0][9] = 4.000000000000000e-01
    Finished in 0.2607285976409912 seconds.

# Step 3: Code validation: Verify that relative error in numerical solution converges to zero at the expected order \[Back to [top](#toc)\]
$$\label{convergence}$$

For testing purposes, we have only checked these algorithms on a small grid. By construction, we have only guaranteed ourselves output from the functions we are testing at a point, so we will simply print the convergence order at that point after processing our outputs below.

```python
import numpy as np
import matplotlib.pyplot as plt

Data1 = np.loadtxt("out2-numer.txt")
Data2 = np.loadtxt("out4-numer.txt")

# print("The following quantities converge at the listed order (should be ~2):")

conv_order_min = 1.9
conv_order_max = 2.3

def check_conv_order(description,conv_order):
    if conv_order < conv_order_min or conv_order > conv_order_max:
        print("ERROR: convergence order of "+description+" was too large: "+str(conv_order))
        sys.exit(1)
    print("PASS: "+description+" has convergence order of "+str(conv_order))

check_conv_order("Stilde_rhsD0", np.log2((Data1[0])/(Data2[0])))
check_conv_order("Stilde_rhsD1", np.log2((Data1[0])/(Data2[0])))

print("Stilde_rhsD0: "+str(np.log2((Data1[0])/(Data2[0]))))
print("Stilde_rhsD1: "+str(np.log2((Data1[1])/(Data2[1]))))
print("Stilde_rhsD2: "+str(np.log2((Data1[2])/(Data2[2]))))
print("A_rhsD0:      "+str(np.log2((Data1[3])/(Data2[3]))))
print("A_rhsD1:      "+str(np.log2((Data1[4])/(Data2[4]))))
print("A_rhsD2:      "+str(np.log2((Data1[5])/(Data2[5]))))
print("psi6Phi:      "+str(np.log2((Data1[6])/(Data2[6]))))

```

    The following quantities converge at the listed order (should be ~2):
    Stilde_rhsD0: 2.0331691730607853
    Stilde_rhsD1: 2.19343762296608
    Stilde_rhsD2: 2.1412625521533175
    A_rhsD0:      2.017299839695867
    A_rhsD1:      2.0210492279403764
    A_rhsD2:      2.0022838026313523
    psi6Phi:      2.0107304427621933

# Step 4: Output this notebook to $\LaTeX$-formatted PDF file \[Back to [top](#toc)\]
$$\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-Start_to_Finish_UnitTest-GiRaFFE_NRPy-Source_Terms.pdf](Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-Source_Terms.pdf) (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

```python
import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-Source_Terms")
```

max_stars_repo_name: ksible/nrpytutorial

max_issues_repo_name: ksible/nrpytutorial

max_forks_repo_name: ksible/nrpytutorial


text:
# Semi-Monocoque Theory

```python
from pint import UnitRegistry
import sympy
import networkx as nx
import numpy as np
import matplotlib.pyplot as plt
import sys
%matplotlib inline
from IPython.display import display
```

Import **Section** class, which contains all calculations

```python
from Section import Section
```

Initialization of **sympy** symbolic tool and **pint** for dimension analysis (not really implemented rn as not directly compatible with sympy)

```python
ureg = UnitRegistry()
sympy.init_printing()
```

Define **sympy** parameters used for geometric description of sections

```python
A, A0, t, t0, a, b, h, L = sympy.symbols('A A_0 t t_0 a b h L', positive=True)
```

We also define numerical values for each **symbol** in order to plot scaled section and perform calculations

```python
values = [(A, 150 * ureg.millimeter**2),(A0, 250  * ureg.millimeter**2),(a, 80 * ureg.millimeter), \
          (b, 20 * ureg.millimeter),(h, 35 * ureg.millimeter),(L, 2000 * ureg.millimeter)]
datav = [(v[0],v[1].magnitude) for v in values]
```

# Triangular section

Define graph describing the section:

1) **stringers** are **nodes** with parameters:
- **x** coordinate
- **y** coordinate
- **Area**

2) **panels** are **oriented edges** with parameters:
- **thickness**
- **lenght** which is automatically calculated

```python
stringers = {1:[(sympy.Integer(0),h),A],
             2:[(sympy.Integer(0),sympy.Integer(0)),A],
             3:[(a,sympy.Integer(0)),A]}

panels = {(1,2):t,
          (2,3):t,
          (3,1):t}
```

Define section and perform first calculations

```python
S1 = Section(stringers, panels)
```

```python
S1.cycles
```

## Plot of **S1** section in original reference frame

Define a dictionary of coordinates used by **Networkx** to plot section as a Directed graph.
Note that arrows are actually just thicker stubs

```python
start_pos={ii: [float(S1.g.node[ii]['ip'][i].subs(datav)) for i in range(2)] for ii in S1.g.nodes() }
```

```python
plt.figure(figsize=(12,8),dpi=300)
nx.draw(S1.g,with_labels=True, arrows= True, pos=start_pos)
plt.arrow(0,0,20,0)
plt.arrow(0,0,0,20)
#plt.text(0,0, 'CG', fontsize=24)
plt.axis('equal')
plt.title("Section in starting reference Frame",fontsize=16);
```

Expression of **Inertial properties** wrt Center of Gravity in with original rotation

```python
S1.Ixx0, S1.Iyy0, S1.Ixy0, S1.α0
```

## Plot of **S1** section in inertial reference Frame

Section is plotted wrt **center of gravity** and rotated (if necessary) so that *x* and *y* are principal axes.
**Center of Gravity** and **Shear Center** are drawn

```python
positions={ii: [float(S1.g.node[ii]['pos'][i].subs(datav)) for i in range(2)] for ii in S1.g.nodes() }
```

```python
x_ct, y_ct = S1.ct.subs(datav)

plt.figure(figsize=(12,8),dpi=300)
nx.draw(S1.g,with_labels=True, pos=positions)
plt.plot([0],[0],'o',ms=12,label='CG')
plt.plot([x_ct],[y_ct],'^',ms=12, label='SC')
#plt.text(0,0, 'CG', fontsize=24)
#plt.text(x_ct,y_ct, 'SC', fontsize=24)
plt.legend(loc='lower right', shadow=True)
plt.axis('equal')
plt.title("Section in pricipal reference Frame",fontsize=16);
```

Expression of **inertial properties** in *principal reference frame*

```python
sympy.simplify(S1.Ixx), sympy.simplify(S1.Iyy), sympy.simplify(S1.Ixy), sympy.simplify(S1.θ)
```

## **Shear center** expression

Expressions can be messy, so we evaluate them to numerical values

```python
sympy.N(S1.ct.subs(datav))
```

## Analisys of Loads

We define some symbols

```python
Tx, Ty, Nz, Mx, My, Mz, F, ry, ry, mz = sympy.symbols('T_x T_y N_z M_x M_y M_z F r_y r_x m_z')
```

```python
S1.set_loads(_Tx=0, _Ty=Ty, _Nz=0, _Mx=Mx, _My=0, _Mz=0)
#S1.compute_stringer_actions()
#S1.compute_panel_fluxes();
```

**Axial Loads**

```python
#S1.N
```

**Panel Fluxes**

```python
#S1.q
```

**Example 2**: _twisting moment_ in **z** direction

```python
S1.set_loads(_Tx=0, _Ty=0, _Nz=0, _Mx=0, _My=0, _Mz=Mz)
S1.compute_stringer_actions()
S1.compute_panel_fluxes();
```

**Axial Loads**

```python
S1.N
```

**Panel Fluxes**
evaluated to numerical values

```python
{k:sympy.N(S1.q[k].subs(datav)) for k in S1.q }
```

## Torsional moment of Inertia

```python
S1.compute_Jt()
```

```python
sympy.N(S1.Jt.subs(datav))
```

```python

```

ext: ipynb


text:
# Chapter 4 Printing

Sympyではキレイに出力することができる！

詳しい内容については:
    http://docs.sympy.org/latest/tutorial/printing.html 参照.

## 4.1 'Pretty Printing'の設定

```python
from sympy import init_printing
```

```python
init_printing()
```

--->自分の環境で最適な出力(`LaTeX`, `MathML`など)を自動的に選ぶ.

```python
from sympy import init_session
```

```python
init_session() #前章参照
```

    IPython console for SymPy 1.0 (Python 3.6.0-64-bit) (ground types: python)

    These commands were executed:
    >>> from __future__ import division
    >>> from sympy import *
    >>> x, y, z, t = symbols('x y z t')
    >>> k, m, n = symbols('k m n', integer=True)
    >>> f, g, h = symbols('f g h', cls=Function)
    >>> init_printing()

    Documentation can be found at http://docs.sympy.org/1.0/

コンソールで実行すると`LaTeX`等で出力され、`jupyter notebook`では

```python
from sympy import *
x, y, z = symbols('x y z')
init_printing()
```

```python
Integral(sqrt(1/x), x)
```

とキレイに出力される！

latexで出力させたいときは

```python
print(latex(Integral(sqrt(1/x), x)))
```

    \int \sqrt{\frac{1}{x}}\, dx

とすればよい.(ちゃんと被積分関数とdxの間に\,が入ってるのがすごい.)

次は[Chapter5 Simplifications](https://hiroyuki827.github.io/SymPy_tutorial/Chapter5_Simplifications.html)へ!