Notebook

Generating C code for the right-hand sides of Maxwell's equations, in curvilinear coordinates, using a reference metric formalism¶

Author: Ian Ruchlin¶

Formatting improvements courtesy Brandon Clark¶

This notebook outlines the generation of C code for Maxwell's equations in curvilinear coordinates, using reference metric formalism. It explores two systems: the first with fixed constraint-violating numerical error and the second employing an auxiliary variable to enhance stability by allowing constraint violations to propagate away.¶

The following formulations of Maxwell's equations, called System I and System II, are described in Illustrating Stability Properties of Numerical Relativity in Electrodynamics by Knapp et al.¶

Notebook Status: In progress

Validation Notes: This module has not yet undergone validation testing. Do *not* use it until after appropriate validation testing has been performed.

Introduction:¶

Maxwell's equations are subject to the Gauss' law constraint

$\mathcal{C} \equiv \hat{D}_{i} E^{i} - 4 \pi \rho = 0 \; ,$ where

$E^{i}$ is the electric vector field,

$\hat{D}_{i}$ is the covariant derivative associated with the reference metric

$\hat{\gamma}_{i j}$ (which is taken to represent flat space), and

$\rho$ is the electric charge density. We use

$\mathcal{C}$ as a measure of numerical error. Maxwell's equations are also required to satisfy

$\hat{D}_{i} B^{i} = 0$ , where

$B^{i}$ is the magnetic vector field. The magnetic constraint implies that the magnetic field can be expressed as

$B_{i} = \epsilon_{i j k} \hat{D}^{j} A^{k} \; ,$ where

$\epsilon_{i j k}$ is the totally antisymmetric Levi-Civita tensor and

$A^{i}$ is the vector potential field. Together with the scalar potential

$\psi$ , the electric field can be expressed in terms of the potential fields as

$E_{i} = -\hat{D}_{i} \psi - \partial_{t} A_{i} \; .$ For now, we work in vacuum, where the electric charge density and the electric current density vector both vanish (

$\rho = 0$ and

$j_{i} = 0$ ).

In addition to the Gauss constraints, the electric and magnetic fields obey two independent electromagnetic invariants

$\begin{align} \mathcal{P} &\equiv B_{i} B^{i} - E_{i} E^{i} \; , \\ \mathcal{Q} &\equiv E_{i} B^{i} \; . \end{align}$ In vacuum, these satisfy

$\mathcal{P} = \mathcal{Q} = 0$ .

Table of Contents:¶

$\label{toc}$

This notebook is organized as follows

Step 1: System I
Step 2: System II
Step 3: Output this notebook to $\LaTeX$ -formatted PDF file

Step 1: System I [Back to top]¶

$\label{sys1}$

In terms of the above definitions, the evolution Maxwell's equations take the form

$\begin{align} \partial_{t} A_{i} &= -E_{i} - \hat{D}_{i} \psi \; , \\ \partial_{t} E_{i} &= -\hat{D}_{j} \hat{D}^{j} A_{i} + \hat{D}_{i} \hat{D}_{j} A^{j}\; , \\ \partial_{t} \psi &= -\hat{D}_{i} A^{i} \; . \end{align}$ Note that this coupled system contains mixed second derivatives in the second term on the right hand side of the

$E^{i}$ evolution equation. We will revisit this fact when building System II.

It can be shown that the Gauss constraint satisfies the evolution equation

$\partial_{t} \mathcal{C} = 0 \; .$ This implies that any constraint violating numerical error remains fixed in place during the evolution. This becomes problematic when the violations grow large and spoil the physics of the simulation.

In [1]:

import NRPy_param_funcs as par   # NRPy+: parameter interface
import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import grid as gri               # NRPy+: Functions having to do with numerical grids
import finite_difference as fin  # NRPy+: Finite difference C code generation module
import reference_metric as rfm   # NRPy+: Reference metric support
from outputC import lhrh         # NRPy+: Core C code output module

par.set_parval_from_str("reference_metric::CoordSystem", "Spherical")
par.set_parval_from_str("grid::DIM", 3)

rfm.reference_metric()

# The name of this module ("maxwell") is given by __name__:
thismodule = __name__

# Step 0: Read the spatial dimension parameter as DIM.
DIM = par.parval_from_str("grid::DIM")

# Step 1: Set the finite differencing order to 4.
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", 4)

# Step 2: Register gridfunctions that are needed as input.
psi = gri.register_gridfunctions("EVOL", ["psi"])

# Step 3a: Declare the rank-1 indexed expressions E_{i}, A_{i},
#          and \partial_{i} \psi. Derivative variables like these
#          must have an underscore in them, so the finite
#          difference module can parse the variable name properly.
ED = ixp.register_gridfunctions_for_single_rank1("EVOL", "ED")
AD = ixp.register_gridfunctions_for_single_rank1("EVOL", "AD")
psi_dD = ixp.declarerank1("psi_dD")

# Step 3b: Declare the rank-2 indexed expression \partial_{j} A_{i},
#          which is not symmetric in its indices.
#          Derivative variables like these must have an underscore
#          in them, so the finite difference module can parse the
#          variable name properly.
AD_dD = ixp.declarerank2("AD_dD", "nosym")

# Step 3c: Declare the rank-3 indexed expression \partial_{jk} A_{i},
#          which is symmetric in the two {jk} indices.
AD_dDD = ixp.declarerank3("AD_dDD", "sym12")

# Step 4: Calculate first and second covariant derivatives, and the
#         necessary contractions.
# First covariant derivative
# D_{j} A_{i} = A_{i,j} - \Gamma^{k}_{ij} A_{k}
AD_dHatD = ixp.zerorank2()
for i in range(DIM):
    for j in range(DIM):
        AD_dHatD[i][j] = AD_dD[i][j]
        for k in range(DIM):
            AD_dHatD[i][j] -= rfm.GammahatUDD[k][i][j] * AD[k]

# Second covariant derivative
# D_{k} D_{j} A_{i} = \partial_{k} D_{j} A_{i} - \Gamma^{l}_{jk} D_{l} A_{i}
#                    - \Gamma^{l}_{ik} D_{j} A_{l}
#                   =  A_{i,jk}
#                    - \Gamma^{l}_{ij,k} A_{l}
#                    - \Gamma^{l}_{ij}   A_{l,k}
#                    - \Gamma^{l}_{jk}   A_{i;\hat{l}}
#                    - \Gamma^{l}_{ik}   A_{l;\hat{j}}
AD_dHatDD = ixp.zerorank3()
for i in range(DIM):
    for j in range(DIM):
        for k in range(DIM):
            AD_dHatDD[i][j][k] = AD_dDD[i][j][k]
            for l in range(DIM):
                AD_dHatDD[i][j][k] += - rfm.GammahatUDDdD[l][i][j][k]  * AD[l]  \
                                      - rfm.GammahatUDD[l][i][j] *    AD_dD[l][k] \
                                      - rfm.GammahatUDD[l][j][k] * AD_dHatD[i][l] \
                                      - rfm.GammahatUDD[l][i][k] * AD_dHatD[l][j]

# Covariant divergence
# D_{i} A^{i} = ghat^{ij} D_{j} A_{i}
DivA = 0
# Gradient of covariant divergence
# DivA_dD_{i} = ghat^{jk} A_{k;\hat{j}\hat{i}}
DivA_dD = ixp.zerorank1()
# Covariant Laplacian
# LapAD_{i} = ghat^{jk} A_{i;\hat{j}\hat{k}}
LapAD = ixp.zerorank1()
for i in range(DIM):
    for j in range(DIM):
        DivA += rfm.ghatUU[i][j] * AD_dHatD[i][j]
        for k in range(DIM):
            DivA_dD[i] += rfm.ghatUU[j][k] * AD_dHatDD[k][j][i]
            LapAD[i]   += rfm.ghatUU[j][k] * AD_dHatDD[i][j][k]

# Step 5: Define right-hand sides for the evolution.
AD_rhs = ixp.zerorank1()
ED_rhs = ixp.zerorank1()
for i in range(DIM):
    AD_rhs[i] = -ED[i] - psi_dD[i]
    ED_rhs[i] = -LapAD[i] + DivA_dD[i]
psi_rhs = -DivA

# Step 6: Generate C code for System I Maxwell's evolution equations,
#         print output to the screen (standard out, or stdout).
lhrh_list = []
for i in range(DIM):
    lhrh_list.append(lhrh(lhs=gri.gfaccess("rhs_gfs", "AD" + str(i)), rhs=AD_rhs[i]))
    lhrh_list.append(lhrh(lhs=gri.gfaccess("rhs_gfs", "ED" + str(i)), rhs=ED_rhs[i]))
lhrh_list.append(lhrh(lhs=gri.gfaccess("rhs_gfs", "psi"), rhs=psi_rhs))

fin.FD_outputC("stdout", lhrh_list)

{
   /*
    * NRPy+ Finite Difference Code Generation, Step 1 of 2: Read from main memory and compute finite difference stencils:
    */
   /*
    *  Original SymPy expressions:
    *  "[const double AD_dD00 = invdx0*(-2*AD0_i0m1_i1_i2/3 + AD0_i0m2_i1_i2/12 + 2*AD0_i0p1_i1_i2/3 - AD0_i0p2_i1_i2/12),
    *    const double AD_dD01 = invdx1*(-2*AD0_i0_i1m1_i2/3 + AD0_i0_i1m2_i2/12 + 2*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),
    *    const double AD_dD02 = invdx2*(-2*AD0_i0_i1_i2m1/3 + AD0_i0_i1_i2m2/12 + 2*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),
    *    const double AD_dD10 = invdx0*(-2*AD1_i0m1_i1_i2/3 + AD1_i0m2_i1_i2/12 + 2*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),
    *    const double AD_dD11 = invdx1*(-2*AD1_i0_i1m1_i2/3 + AD1_i0_i1m2_i2/12 + 2*AD1_i0_i1p1_i2/3 - AD1_i0_i1p2_i2/12),
    *    const double AD_dD12 = invdx2*(-2*AD1_i0_i1_i2m1/3 + AD1_i0_i1_i2m2/12 + 2*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),
    *    const double AD_dD20 = invdx0*(-2*AD2_i0m1_i1_i2/3 + AD2_i0m2_i1_i2/12 + 2*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),
    *    const double AD_dD21 = invdx1*(-2*AD2_i0_i1m1_i2/3 + AD2_i0_i1m2_i2/12 + 2*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),
    *    const double AD_dD22 = invdx2*(-2*AD2_i0_i1_i2m1/3 + AD2_i0_i1_i2m2/12 + 2*AD2_i0_i1_i2p1/3 - AD2_i0_i1_i2p2/12),
    *    const double AD_dDD001 = invdx0*invdx1*(4*AD0_i0m1_i1m1_i2/9 - AD0_i0m1_i1m2_i2/18 - 4*AD0_i0m1_i1p1_i2/9 + AD0_i0m1_i1p2_i2/18 - AD0_i0m2_i1m1_i2/18 + AD0_i0m2_i1m2_i2/144 + AD0_i0m2_i1p1_i2/18 - AD0_i0m2_i1p2_i2/144 - 4*AD0_i0p1_i1m1_i2/9 + AD0_i0p1_i1m2_i2/18 + 4*AD0_i0p1_i1p1_i2/9 - AD0_i0p1_i1p2_i2/18 + AD0_i0p2_i1m1_i2/18 - AD0_i0p2_i1m2_i2/144 - AD0_i0p2_i1p1_i2/18 + AD0_i0p2_i1p2_i2/144),
    *    const double AD_dDD002 = invdx0*invdx2*(4*AD0_i0m1_i1_i2m1/9 - AD0_i0m1_i1_i2m2/18 - 4*AD0_i0m1_i1_i2p1/9 + AD0_i0m1_i1_i2p2/18 - AD0_i0m2_i1_i2m1/18 + AD0_i0m2_i1_i2m2/144 + AD0_i0m2_i1_i2p1/18 - AD0_i0m2_i1_i2p2/144 - 4*AD0_i0p1_i1_i2m1/9 + AD0_i0p1_i1_i2m2/18 + 4*AD0_i0p1_i1_i2p1/9 - AD0_i0p1_i1_i2p2/18 + AD0_i0p2_i1_i2m1/18 - AD0_i0p2_i1_i2m2/144 - AD0_i0p2_i1_i2p1/18 + AD0_i0p2_i1_i2p2/144),
    *    const double AD_dDD011 = invdx1**2*(-5*AD0/2 + 4*AD0_i0_i1m1_i2/3 - AD0_i0_i1m2_i2/12 + 4*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),
    *    const double AD_dDD022 = invdx2**2*(-5*AD0/2 + 4*AD0_i0_i1_i2m1/3 - AD0_i0_i1_i2m2/12 + 4*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),
    *    const double AD_dDD100 = invdx0**2*(-5*AD1/2 + 4*AD1_i0m1_i1_i2/3 - AD1_i0m2_i1_i2/12 + 4*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),
    *    const double AD_dDD101 = invdx0*invdx1*(4*AD1_i0m1_i1m1_i2/9 - AD1_i0m1_i1m2_i2/18 - 4*AD1_i0m1_i1p1_i2/9 + AD1_i0m1_i1p2_i2/18 - AD1_i0m2_i1m1_i2/18 + AD1_i0m2_i1m2_i2/144 + AD1_i0m2_i1p1_i2/18 - AD1_i0m2_i1p2_i2/144 - 4*AD1_i0p1_i1m1_i2/9 + AD1_i0p1_i1m2_i2/18 + 4*AD1_i0p1_i1p1_i2/9 - AD1_i0p1_i1p2_i2/18 + AD1_i0p2_i1m1_i2/18 - AD1_i0p2_i1m2_i2/144 - AD1_i0p2_i1p1_i2/18 + AD1_i0p2_i1p2_i2/144),
    *    const double AD_dDD112 = invdx1*invdx2*(4*AD1_i0_i1m1_i2m1/9 - AD1_i0_i1m1_i2m2/18 - 4*AD1_i0_i1m1_i2p1/9 + AD1_i0_i1m1_i2p2/18 - AD1_i0_i1m2_i2m1/18 + AD1_i0_i1m2_i2m2/144 + AD1_i0_i1m2_i2p1/18 - AD1_i0_i1m2_i2p2/144 - 4*AD1_i0_i1p1_i2m1/9 + AD1_i0_i1p1_i2m2/18 + 4*AD1_i0_i1p1_i2p1/9 - AD1_i0_i1p1_i2p2/18 + AD1_i0_i1p2_i2m1/18 - AD1_i0_i1p2_i2m2/144 - AD1_i0_i1p2_i2p1/18 + AD1_i0_i1p2_i2p2/144),
    *    const double AD_dDD122 = invdx2**2*(-5*AD1/2 + 4*AD1_i0_i1_i2m1/3 - AD1_i0_i1_i2m2/12 + 4*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),
    *    const double AD_dDD200 = invdx0**2*(-5*AD2/2 + 4*AD2_i0m1_i1_i2/3 - AD2_i0m2_i1_i2/12 + 4*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),
    *    const double AD_dDD202 = invdx0*invdx2*(4*AD2_i0m1_i1_i2m1/9 - AD2_i0m1_i1_i2m2/18 - 4*AD2_i0m1_i1_i2p1/9 + AD2_i0m1_i1_i2p2/18 - AD2_i0m2_i1_i2m1/18 + AD2_i0m2_i1_i2m2/144 + AD2_i0m2_i1_i2p1/18 - AD2_i0m2_i1_i2p2/144 - 4*AD2_i0p1_i1_i2m1/9 + AD2_i0p1_i1_i2m2/18 + 4*AD2_i0p1_i1_i2p1/9 - AD2_i0p1_i1_i2p2/18 + AD2_i0p2_i1_i2m1/18 - AD2_i0p2_i1_i2m2/144 - AD2_i0p2_i1_i2p1/18 + AD2_i0p2_i1_i2p2/144),
    *    const double AD_dDD211 = invdx1**2*(-5*AD2/2 + 4*AD2_i0_i1m1_i2/3 - AD2_i0_i1m2_i2/12 + 4*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),
    *    const double AD_dDD212 = invdx1*invdx2*(4*AD2_i0_i1m1_i2m1/9 - AD2_i0_i1m1_i2m2/18 - 4*AD2_i0_i1m1_i2p1/9 + AD2_i0_i1m1_i2p2/18 - AD2_i0_i1m2_i2m1/18 + AD2_i0_i1m2_i2m2/144 + AD2_i0_i1m2_i2p1/18 - AD2_i0_i1m2_i2p2/144 - 4*AD2_i0_i1p1_i2m1/9 + AD2_i0_i1p1_i2m2/18 + 4*AD2_i0_i1p1_i2p1/9 - AD2_i0_i1p1_i2p2/18 + AD2_i0_i1p2_i2m1/18 - AD2_i0_i1p2_i2m2/144 - AD2_i0_i1p2_i2p1/18 + AD2_i0_i1p2_i2p2/144),
    *    const double psi_dD0 = invdx0*(-2*psi_i0m1_i1_i2/3 + psi_i0m2_i1_i2/12 + 2*psi_i0p1_i1_i2/3 - psi_i0p2_i1_i2/12),
    *    const double psi_dD1 = invdx1*(-2*psi_i0_i1m1_i2/3 + psi_i0_i1m2_i2/12 + 2*psi_i0_i1p1_i2/3 - psi_i0_i1p2_i2/12),
    *    const double psi_dD2 = invdx2*(-2*psi_i0_i1_i2m1/3 + psi_i0_i1_i2m2/12 + 2*psi_i0_i1_i2p1/3 - psi_i0_i1_i2p2/12)]"
    */
   const double psi_i0_i1_i2m2 = in_gfs[IDX4(PSIGF, i0,i1,i2-2)];
   const double psi_i0_i1_i2m1 = in_gfs[IDX4(PSIGF, i0,i1,i2-1)];
   const double psi_i0_i1m2_i2 = in_gfs[IDX4(PSIGF, i0,i1-2,i2)];
   const double psi_i0_i1m1_i2 = in_gfs[IDX4(PSIGF, i0,i1-1,i2)];
   const double psi_i0m2_i1_i2 = in_gfs[IDX4(PSIGF, i0-2,i1,i2)];
   const double psi_i0m1_i1_i2 = in_gfs[IDX4(PSIGF, i0-1,i1,i2)];
   const double psi_i0p1_i1_i2 = in_gfs[IDX4(PSIGF, i0+1,i1,i2)];
   const double psi_i0p2_i1_i2 = in_gfs[IDX4(PSIGF, i0+2,i1,i2)];
   const double psi_i0_i1p1_i2 = in_gfs[IDX4(PSIGF, i0,i1+1,i2)];
   const double psi_i0_i1p2_i2 = in_gfs[IDX4(PSIGF, i0,i1+2,i2)];
   const double psi_i0_i1_i2p1 = in_gfs[IDX4(PSIGF, i0,i1,i2+1)];
   const double psi_i0_i1_i2p2 = in_gfs[IDX4(PSIGF, i0,i1,i2+2)];
   const double ED0 = in_gfs[IDX4(ED0GF, i0,i1,i2)];
   const double ED1 = in_gfs[IDX4(ED1GF, i0,i1,i2)];
   const double ED2 = in_gfs[IDX4(ED2GF, i0,i1,i2)];
   const double AD0_i0m2_i1_i2m2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2-2)];
   const double AD0_i0m1_i1_i2m2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2-2)];
   const double AD0_i0_i1_i2m2 = in_gfs[IDX4(AD0GF, i0,i1,i2-2)];
   const double AD0_i0p1_i1_i2m2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2-2)];
   const double AD0_i0p2_i1_i2m2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2-2)];
   const double AD0_i0m2_i1_i2m1 = in_gfs[IDX4(AD0GF, i0-2,i1,i2-1)];
   const double AD0_i0m1_i1_i2m1 = in_gfs[IDX4(AD0GF, i0-1,i1,i2-1)];
   const double AD0_i0_i1_i2m1 = in_gfs[IDX4(AD0GF, i0,i1,i2-1)];
   const double AD0_i0p1_i1_i2m1 = in_gfs[IDX4(AD0GF, i0+1,i1,i2-1)];
   const double AD0_i0p2_i1_i2m1 = in_gfs[IDX4(AD0GF, i0+2,i1,i2-1)];
   const double AD0_i0m2_i1m2_i2 = in_gfs[IDX4(AD0GF, i0-2,i1-2,i2)];
   const double AD0_i0m1_i1m2_i2 = in_gfs[IDX4(AD0GF, i0-1,i1-2,i2)];
   const double AD0_i0_i1m2_i2 = in_gfs[IDX4(AD0GF, i0,i1-2,i2)];
   const double AD0_i0p1_i1m2_i2 = in_gfs[IDX4(AD0GF, i0+1,i1-2,i2)];
   const double AD0_i0p2_i1m2_i2 = in_gfs[IDX4(AD0GF, i0+2,i1-2,i2)];
   const double AD0_i0m2_i1m1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1-1,i2)];
   const double AD0_i0m1_i1m1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1-1,i2)];
   const double AD0_i0_i1m1_i2 = in_gfs[IDX4(AD0GF, i0,i1-1,i2)];
   const double AD0_i0p1_i1m1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1-1,i2)];
   const double AD0_i0p2_i1m1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1-1,i2)];
   const double AD0_i0m2_i1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2)];
   const double AD0_i0m1_i1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2)];
   const double AD0 = in_gfs[IDX4(AD0GF, i0,i1,i2)];
   const double AD0_i0p1_i1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2)];
   const double AD0_i0p2_i1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2)];
   const double AD0_i0m2_i1p1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1+1,i2)];
   const double AD0_i0m1_i1p1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1+1,i2)];
   const double AD0_i0_i1p1_i2 = in_gfs[IDX4(AD0GF, i0,i1+1,i2)];
   const double AD0_i0p1_i1p1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1+1,i2)];
   const double AD0_i0p2_i1p1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1+1,i2)];
   const double AD0_i0m2_i1p2_i2 = in_gfs[IDX4(AD0GF, i0-2,i1+2,i2)];
   const double AD0_i0m1_i1p2_i2 = in_gfs[IDX4(AD0GF, i0-1,i1+2,i2)];
   const double AD0_i0_i1p2_i2 = in_gfs[IDX4(AD0GF, i0,i1+2,i2)];
   const double AD0_i0p1_i1p2_i2 = in_gfs[IDX4(AD0GF, i0+1,i1+2,i2)];
   const double AD0_i0p2_i1p2_i2 = in_gfs[IDX4(AD0GF, i0+2,i1+2,i2)];
   const double AD0_i0m2_i1_i2p1 = in_gfs[IDX4(AD0GF, i0-2,i1,i2+1)];
   const double AD0_i0m1_i1_i2p1 = in_gfs[IDX4(AD0GF, i0-1,i1,i2+1)];
   const double AD0_i0_i1_i2p1 = in_gfs[IDX4(AD0GF, i0,i1,i2+1)];
   const double AD0_i0p1_i1_i2p1 = in_gfs[IDX4(AD0GF, i0+1,i1,i2+1)];
   const double AD0_i0p2_i1_i2p1 = in_gfs[IDX4(AD0GF, i0+2,i1,i2+1)];
   const double AD0_i0m2_i1_i2p2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2+2)];
   const double AD0_i0m1_i1_i2p2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2+2)];
   const double AD0_i0_i1_i2p2 = in_gfs[IDX4(AD0GF, i0,i1,i2+2)];
   const double AD0_i0p1_i1_i2p2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2+2)];
   const double AD0_i0p2_i1_i2p2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2+2)];
   const double AD1_i0_i1m2_i2m2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2-2)];
   const double AD1_i0_i1m1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2-2)];
   const double AD1_i0_i1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1,i2-2)];
   const double AD1_i0_i1p1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2-2)];
   const double AD1_i0_i1p2_i2m2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2-2)];
   const double AD1_i0_i1m2_i2m1 = in_gfs[IDX4(AD1GF, i0,i1-2,i2-1)];
   const double AD1_i0_i1m1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1-1,i2-1)];
   const double AD1_i0_i1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1,i2-1)];
   const double AD1_i0_i1p1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1+1,i2-1)];
   const double AD1_i0_i1p2_i2m1 = in_gfs[IDX4(AD1GF, i0,i1+2,i2-1)];
   const double AD1_i0m2_i1m2_i2 = in_gfs[IDX4(AD1GF, i0-2,i1-2,i2)];
   const double AD1_i0m1_i1m2_i2 = in_gfs[IDX4(AD1GF, i0-1,i1-2,i2)];
   const double AD1_i0_i1m2_i2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2)];
   const double AD1_i0p1_i1m2_i2 = in_gfs[IDX4(AD1GF, i0+1,i1-2,i2)];
   const double AD1_i0p2_i1m2_i2 = in_gfs[IDX4(AD1GF, i0+2,i1-2,i2)];
   const double AD1_i0m2_i1m1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1-1,i2)];
   const double AD1_i0m1_i1m1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1-1,i2)];
   const double AD1_i0_i1m1_i2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2)];
   const double AD1_i0p1_i1m1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1-1,i2)];
   const double AD1_i0p2_i1m1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1-1,i2)];
   const double AD1_i0m2_i1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1,i2)];
   const double AD1_i0m1_i1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1,i2)];
   const double AD1 = in_gfs[IDX4(AD1GF, i0,i1,i2)];
   const double AD1_i0p1_i1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1,i2)];
   const double AD1_i0p2_i1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1,i2)];
   const double AD1_i0m2_i1p1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1+1,i2)];
   const double AD1_i0m1_i1p1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1+1,i2)];
   const double AD1_i0_i1p1_i2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2)];
   const double AD1_i0p1_i1p1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1+1,i2)];
   const double AD1_i0p2_i1p1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1+1,i2)];
   const double AD1_i0m2_i1p2_i2 = in_gfs[IDX4(AD1GF, i0-2,i1+2,i2)];
   const double AD1_i0m1_i1p2_i2 = in_gfs[IDX4(AD1GF, i0-1,i1+2,i2)];
   const double AD1_i0_i1p2_i2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2)];
   const double AD1_i0p1_i1p2_i2 = in_gfs[IDX4(AD1GF, i0+1,i1+2,i2)];
   const double AD1_i0p2_i1p2_i2 = in_gfs[IDX4(AD1GF, i0+2,i1+2,i2)];
   const double AD1_i0_i1m2_i2p1 = in_gfs[IDX4(AD1GF, i0,i1-2,i2+1)];
   const double AD1_i0_i1m1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1-1,i2+1)];
   const double AD1_i0_i1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1,i2+1)];
   const double AD1_i0_i1p1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1+1,i2+1)];
   const double AD1_i0_i1p2_i2p1 = in_gfs[IDX4(AD1GF, i0,i1+2,i2+1)];
   const double AD1_i0_i1m2_i2p2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2+2)];
   const double AD1_i0_i1m1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2+2)];
   const double AD1_i0_i1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1,i2+2)];
   const double AD1_i0_i1p1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2+2)];
   const double AD1_i0_i1p2_i2p2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2+2)];
   const double AD2_i0_i1m2_i2m2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2-2)];
   const double AD2_i0_i1m1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2-2)];
   const double AD2_i0m2_i1_i2m2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2-2)];
   const double AD2_i0m1_i1_i2m2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2-2)];
   const double AD2_i0_i1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1,i2-2)];
   const double AD2_i0p1_i1_i2m2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2-2)];
   const double AD2_i0p2_i1_i2m2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2-2)];
   const double AD2_i0_i1p1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2-2)];
   const double AD2_i0_i1p2_i2m2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2-2)];
   const double AD2_i0_i1m2_i2m1 = in_gfs[IDX4(AD2GF, i0,i1-2,i2-1)];
   const double AD2_i0_i1m1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1-1,i2-1)];
   const double AD2_i0m2_i1_i2m1 = in_gfs[IDX4(AD2GF, i0-2,i1,i2-1)];
   const double AD2_i0m1_i1_i2m1 = in_gfs[IDX4(AD2GF, i0-1,i1,i2-1)];
   const double AD2_i0_i1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1,i2-1)];
   const double AD2_i0p1_i1_i2m1 = in_gfs[IDX4(AD2GF, i0+1,i1,i2-1)];
   const double AD2_i0p2_i1_i2m1 = in_gfs[IDX4(AD2GF, i0+2,i1,i2-1)];
   const double AD2_i0_i1p1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1+1,i2-1)];
   const double AD2_i0_i1p2_i2m1 = in_gfs[IDX4(AD2GF, i0,i1+2,i2-1)];
   const double AD2_i0_i1m2_i2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2)];
   const double AD2_i0_i1m1_i2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2)];
   const double AD2_i0m2_i1_i2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2)];
   const double AD2_i0m1_i1_i2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2)];
   const double AD2 = in_gfs[IDX4(AD2GF, i0,i1,i2)];
   const double AD2_i0p1_i1_i2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2)];
   const double AD2_i0p2_i1_i2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2)];
   const double AD2_i0_i1p1_i2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2)];
   const double AD2_i0_i1p2_i2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2)];
   const double AD2_i0_i1m2_i2p1 = in_gfs[IDX4(AD2GF, i0,i1-2,i2+1)];
   const double AD2_i0_i1m1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1-1,i2+1)];
   const double AD2_i0m2_i1_i2p1 = in_gfs[IDX4(AD2GF, i0-2,i1,i2+1)];
   const double AD2_i0m1_i1_i2p1 = in_gfs[IDX4(AD2GF, i0-1,i1,i2+1)];
   const double AD2_i0_i1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1,i2+1)];
   const double AD2_i0p1_i1_i2p1 = in_gfs[IDX4(AD2GF, i0+1,i1,i2+1)];
   const double AD2_i0p2_i1_i2p1 = in_gfs[IDX4(AD2GF, i0+2,i1,i2+1)];
   const double AD2_i0_i1p1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1+1,i2+1)];
   const double AD2_i0_i1p2_i2p1 = in_gfs[IDX4(AD2GF, i0,i1+2,i2+1)];
   const double AD2_i0_i1m2_i2p2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2+2)];
   const double AD2_i0_i1m1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2+2)];
   const double AD2_i0m2_i1_i2p2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2+2)];
   const double AD2_i0m1_i1_i2p2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2+2)];
   const double AD2_i0_i1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1,i2+2)];
   const double AD2_i0p1_i1_i2p2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2+2)];
   const double AD2_i0p2_i1_i2p2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2+2)];
   const double AD2_i0_i1p1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2+2)];
   const double AD2_i0_i1p2_i2p2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2+2)];
   const double FDPart1_Rational_2_3 = 2.0/3.0;
   const double FDPart1_Rational_1_12 = 1.0/12.0;
   const double FDPart1_Rational_4_9 = 4.0/9.0;
   const double FDPart1_Rational_1_18 = 1.0/18.0;
   const double FDPart1_Rational_1_144 = 1.0/144.0;
   const double FDPart1_Rational_5_2 = 5.0/2.0;
   const double FDPart1_Rational_4_3 = 4.0/3.0;
   const double FDPart1_1 = -AD0_i0_i1_i2p2;
   const double FDPart1_9 = -AD0*FDPart1_Rational_5_2;
   const double FDPart1_12 = -AD1*FDPart1_Rational_5_2;
   const double FDPart1_14 = -AD2*FDPart1_Rational_5_2;
   const double AD_dD00 = invdx0*(FDPart1_Rational_1_12*(AD0_i0m2_i1_i2 - AD0_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD0_i0m1_i1_i2 + AD0_i0p1_i1_i2));
   const double AD_dD01 = invdx1*(FDPart1_Rational_1_12*(AD0_i0_i1m2_i2 - AD0_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));
   const double AD_dD02 = invdx2*(FDPart1_Rational_1_12*(AD0_i0_i1_i2m2 + FDPart1_1) + FDPart1_Rational_2_3*(-AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));
   const double AD_dD10 = invdx0*(FDPart1_Rational_1_12*(AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));
   const double AD_dD11 = invdx1*(FDPart1_Rational_1_12*(AD1_i0_i1m2_i2 - AD1_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD1_i0_i1m1_i2 + AD1_i0_i1p1_i2));
   const double AD_dD12 = invdx2*(FDPart1_Rational_1_12*(AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));
   const double AD_dD20 = invdx0*(FDPart1_Rational_1_12*(AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));
   const double AD_dD21 = invdx1*(FDPart1_Rational_1_12*(AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));
   const double AD_dD22 = invdx2*(FDPart1_Rational_1_12*(AD2_i0_i1_i2m2 - AD2_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD2_i0_i1_i2m1 + AD2_i0_i1_i2p1));
   const double AD_dDD001 = invdx0*invdx1*(FDPart1_Rational_1_144*(AD0_i0m2_i1m2_i2 - AD0_i0m2_i1p2_i2 - AD0_i0p2_i1m2_i2 + AD0_i0p2_i1p2_i2) + FDPart1_Rational_1_18*(-AD0_i0m1_i1m2_i2 + AD0_i0m1_i1p2_i2 - AD0_i0m2_i1m1_i2 + AD0_i0m2_i1p1_i2 + AD0_i0p1_i1m2_i2 - AD0_i0p1_i1p2_i2 + AD0_i0p2_i1m1_i2 - AD0_i0p2_i1p1_i2) + FDPart1_Rational_4_9*(AD0_i0m1_i1m1_i2 - AD0_i0m1_i1p1_i2 - AD0_i0p1_i1m1_i2 + AD0_i0p1_i1p1_i2));
   const double AD_dDD002 = invdx0*invdx2*(FDPart1_Rational_1_144*(AD0_i0m2_i1_i2m2 - AD0_i0m2_i1_i2p2 - AD0_i0p2_i1_i2m2 + AD0_i0p2_i1_i2p2) + FDPart1_Rational_1_18*(-AD0_i0m1_i1_i2m2 + AD0_i0m1_i1_i2p2 - AD0_i0m2_i1_i2m1 + AD0_i0m2_i1_i2p1 + AD0_i0p1_i1_i2m2 - AD0_i0p1_i1_i2p2 + AD0_i0p2_i1_i2m1 - AD0_i0p2_i1_i2p1) + FDPart1_Rational_4_9*(AD0_i0m1_i1_i2m1 - AD0_i0m1_i1_i2p1 - AD0_i0p1_i1_i2m1 + AD0_i0p1_i1_i2p1));
   const double AD_dDD011 = ((invdx1)*(invdx1))*(FDPart1_9 + FDPart1_Rational_1_12*(-AD0_i0_i1m2_i2 - AD0_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));
   const double AD_dDD022 = ((invdx2)*(invdx2))*(FDPart1_9 + FDPart1_Rational_1_12*(-AD0_i0_i1_i2m2 + FDPart1_1) + FDPart1_Rational_4_3*(AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));
   const double AD_dDD100 = ((invdx0)*(invdx0))*(FDPart1_12 + FDPart1_Rational_1_12*(-AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));
   const double AD_dDD101 = invdx0*invdx1*(FDPart1_Rational_1_144*(AD1_i0m2_i1m2_i2 - AD1_i0m2_i1p2_i2 - AD1_i0p2_i1m2_i2 + AD1_i0p2_i1p2_i2) + FDPart1_Rational_1_18*(-AD1_i0m1_i1m2_i2 + AD1_i0m1_i1p2_i2 - AD1_i0m2_i1m1_i2 + AD1_i0m2_i1p1_i2 + AD1_i0p1_i1m2_i2 - AD1_i0p1_i1p2_i2 + AD1_i0p2_i1m1_i2 - AD1_i0p2_i1p1_i2) + FDPart1_Rational_4_9*(AD1_i0m1_i1m1_i2 - AD1_i0m1_i1p1_i2 - AD1_i0p1_i1m1_i2 + AD1_i0p1_i1p1_i2));
   const double AD_dDD112 = invdx1*invdx2*(FDPart1_Rational_1_144*(AD1_i0_i1m2_i2m2 - AD1_i0_i1m2_i2p2 - AD1_i0_i1p2_i2m2 + AD1_i0_i1p2_i2p2) + FDPart1_Rational_1_18*(-AD1_i0_i1m1_i2m2 + AD1_i0_i1m1_i2p2 - AD1_i0_i1m2_i2m1 + AD1_i0_i1m2_i2p1 + AD1_i0_i1p1_i2m2 - AD1_i0_i1p1_i2p2 + AD1_i0_i1p2_i2m1 - AD1_i0_i1p2_i2p1) + FDPart1_Rational_4_9*(AD1_i0_i1m1_i2m1 - AD1_i0_i1m1_i2p1 - AD1_i0_i1p1_i2m1 + AD1_i0_i1p1_i2p1));
   const double AD_dDD122 = ((invdx2)*(invdx2))*(FDPart1_12 + FDPart1_Rational_1_12*(-AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));
   const double AD_dDD200 = ((invdx0)*(invdx0))*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));
   const double AD_dDD202 = invdx0*invdx2*(FDPart1_Rational_1_144*(AD2_i0m2_i1_i2m2 - AD2_i0m2_i1_i2p2 - AD2_i0p2_i1_i2m2 + AD2_i0p2_i1_i2p2) + FDPart1_Rational_1_18*(-AD2_i0m1_i1_i2m2 + AD2_i0m1_i1_i2p2 - AD2_i0m2_i1_i2m1 + AD2_i0m2_i1_i2p1 + AD2_i0p1_i1_i2m2 - AD2_i0p1_i1_i2p2 + AD2_i0p2_i1_i2m1 - AD2_i0p2_i1_i2p1) + FDPart1_Rational_4_9*(AD2_i0m1_i1_i2m1 - AD2_i0m1_i1_i2p1 - AD2_i0p1_i1_i2m1 + AD2_i0p1_i1_i2p1));
   const double AD_dDD211 = ((invdx1)*(invdx1))*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));
   const double AD_dDD212 = invdx1*invdx2*(FDPart1_Rational_1_144*(AD2_i0_i1m2_i2m2 - AD2_i0_i1m2_i2p2 - AD2_i0_i1p2_i2m2 + AD2_i0_i1p2_i2p2) + FDPart1_Rational_1_18*(-AD2_i0_i1m1_i2m2 + AD2_i0_i1m1_i2p2 - AD2_i0_i1m2_i2m1 + AD2_i0_i1m2_i2p1 + AD2_i0_i1p1_i2m2 - AD2_i0_i1p1_i2p2 + AD2_i0_i1p2_i2m1 - AD2_i0_i1p2_i2p1) + FDPart1_Rational_4_9*(AD2_i0_i1m1_i2m1 - AD2_i0_i1m1_i2p1 - AD2_i0_i1p1_i2m1 + AD2_i0_i1p1_i2p1));
   const double psi_dD0 = invdx0*(FDPart1_Rational_1_12*(psi_i0m2_i1_i2 - psi_i0p2_i1_i2) + FDPart1_Rational_2_3*(-psi_i0m1_i1_i2 + psi_i0p1_i1_i2));
   const double psi_dD1 = invdx1*(FDPart1_Rational_1_12*(psi_i0_i1m2_i2 - psi_i0_i1p2_i2) + FDPart1_Rational_2_3*(-psi_i0_i1m1_i2 + psi_i0_i1p1_i2));
   const double psi_dD2 = invdx2*(FDPart1_Rational_1_12*(psi_i0_i1_i2m2 - psi_i0_i1_i2p2) + FDPart1_Rational_2_3*(-psi_i0_i1_i2m1 + psi_i0_i1_i2p1));
   /*
    * NRPy+ Finite Difference Code Generation, Step 2 of 2: Evaluate SymPy expressions and write to main memory:
    */
   /*
    *  Original SymPy expressions:
    *  "[rhs_gfs[IDX4(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0,
    *    rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = (AD0 + AD_dD00*xx0 + AD_dDD101 - 2*(AD0*xx0 + AD_dD11)/xx0)/xx0**2 - (AD_dD00*xx0 - AD_dD11/xx0 + AD_dDD011 - (AD0*xx0 + AD_dD11)/xx0)/xx0**2 + (AD0*sin(xx1)**2 + AD_dD00*xx0*sin(xx1)**2 + AD_dD10*sin(2*xx1)/2 + AD_dDD202 - 2*(AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/xx0)/(xx0**2*sin(xx1)**2) - (AD_dD00*xx0*sin(xx1)**2 - AD_dD22/xx0 + AD_dDD022 + (-AD1/xx0 + AD_dD01)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/xx0)/(xx0**2*sin(xx1)**2),
    *    rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1,
    *    rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1/xx0**2 + AD_dD10/xx0 + AD_dDD001 - AD_dDD100 - (-AD1/xx0 + AD_dD01)/xx0 - (-AD_dD22*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD122 + xx0*(-AD1/xx0 + AD_dD10)*sin(xx1)**2 + (AD0*xx0 + AD_dD11)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)*sin(2*xx1)/(2*sin(xx1)**2))/(xx0**2*sin(xx1)**2) + (2*AD0*xx0*sin(xx1)*cos(xx1) + AD1*cos(2*xx1) + AD_dD01*xx0*sin(xx1)**2 + AD_dD11*sin(2*xx1)/2 + AD_dDD212 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)*sin(2*xx1)/sin(xx1)**2)/(xx0**2*sin(xx1)**2),
    *    rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2,
    *    rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2/xx0**2 + AD_dD20/xx0 + AD_dDD002 - AD_dDD200 - (-AD2/xx0 + AD_dD02)/xx0 + (AD_dD02*xx0 + AD_dDD112 - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD12)*sin(2*xx1)/(2*sin(xx1)**2) - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/(2*sin(xx1)**2))/xx0**2 - (AD2*(-cos(2*xx1)/sin(xx1)**2 + sin(2*xx1)*cos(xx1)/sin(xx1)**3) - AD_dD21*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD211 + xx0*(-AD2/xx0 + AD_dD20) - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/(2*sin(xx1)**2))/xx0**2,
    *    rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -AD_dD00 - (AD0*xx0 + AD_dD11)/xx0**2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/(xx0**2*sin(xx1)**2)]"
    */
   const double FDPart3_0 = (1.0/((xx0)*(xx0)));
   const double FDPart3_1 = AD_dD00*xx0;
   const double FDPart3_2 = (1.0/(xx0));
   const double FDPart3_3 = AD0*xx0;
   const double FDPart3_4 = AD_dD11 + FDPart3_3;
   const double FDPart3_6 = sin(xx1);
   const double FDPart3_7 = ((FDPart3_6)*(FDPart3_6));
   const double FDPart3_9 = sin(2*xx1);
   const double FDPart3_10 = (1.0/2.0)*FDPart3_9;
   const double FDPart3_12 = AD1*FDPart3_10 + AD_dD22 + FDPart3_3*FDPart3_7;
   const double FDPart3_14 = (1.0/(FDPart3_7));
   const double FDPart3_15 = FDPart3_0*FDPart3_14;
   const double FDPart3_16 = -AD1*FDPart3_2;
   const double FDPart3_18 = cos(2*xx1);
   const double FDPart3_20 = cos(xx1);
   const double FDPart3_22 = FDPart3_10*FDPart3_14;
   const double FDPart3_23 = -AD2*FDPart3_2;
   const double FDPart3_24 = -AD2*FDPart3_22;
   const double FDPart3_25 = -FDPart3_22*(AD_dD21 + FDPart3_24);
   rhs_gfs[IDX4(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0;
   rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = FDPart3_0*(AD0 + AD_dDD101 + FDPart3_1 - 2*FDPart3_2*FDPart3_4) - FDPart3_0*(-AD_dD11*FDPart3_2 + AD_dDD011 + FDPart3_1 - FDPart3_2*FDPart3_4) + FDPart3_15*(AD0*FDPart3_7 + AD_dD10*FDPart3_10 + AD_dDD202 + FDPart3_1*FDPart3_7 - 2*FDPart3_12*FDPart3_2) - FDPart3_15*(-AD_dD22*FDPart3_2 + AD_dDD022 + FDPart3_1*FDPart3_7 + FDPart3_10*(AD_dD01 + FDPart3_16) - FDPart3_12*FDPart3_2);
   rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1;
   rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1*FDPart3_0 + AD_dD10*FDPart3_2 + AD_dDD001 - AD_dDD100 - FDPart3_15*(-AD_dD22*FDPart3_22 + AD_dDD122 - FDPart3_10*FDPart3_12*FDPart3_14 + FDPart3_10*FDPart3_4 + FDPart3_7*xx0*(AD_dD10 + FDPart3_16)) + FDPart3_15*(AD1*FDPart3_18 + AD_dD01*FDPart3_7*xx0 + AD_dD11*FDPart3_10 + AD_dDD212 - FDPart3_12*FDPart3_14*FDPart3_9 + 2*FDPart3_20*FDPart3_3*FDPart3_6) - FDPart3_2*(AD_dD01 + FDPart3_16);
   rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2;
   rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2*FDPart3_0 + AD_dD20*FDPart3_2 + AD_dDD002 - AD_dDD200 + FDPart3_0*(AD_dD02*xx0 + AD_dDD112 - FDPart3_22*(AD_dD12 + FDPart3_24) + FDPart3_25) - FDPart3_0*(AD2*(-FDPart3_14*FDPart3_18 + FDPart3_20*FDPart3_9/((FDPart3_6)*(FDPart3_6)*(FDPart3_6))) - AD_dD21*FDPart3_22 + AD_dDD211 + FDPart3_25 + xx0*(AD_dD20 + FDPart3_23)) - FDPart3_2*(AD_dD02 + FDPart3_23);
   rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -AD_dD00 - FDPart3_0*FDPart3_4 - FDPart3_12*FDPart3_15;
}

Step 2: System II [Back to top]¶

$\label{sys2}$

Define the auxiliary variable

$\Gamma \equiv \hat{D}_{i} A^{i} \; .$ Substituting this into Maxwell's equations yields the system

$\begin{align} \partial_{t} A_{i} &= -E_{i} - \hat{D}_{i} \psi \; , \\ \partial_{t} E_{i} &= -\hat{D}_{j} \hat{D}^{j} A_{i} + \hat{D}_{i} \Gamma \; , \\ \partial_{t} \psi &= -\Gamma \; , \\ \partial_{t} \Gamma &= -\hat{D}_{i} \hat{D}^{i} \psi \; . \end{align}$

It can be shown that the Gauss constraint now satisfies the wave equation

$\partial_{t}^{2} \mathcal{C} = \hat{D}_{i} \hat{D}^{i} \mathcal{C} \; .$ Thus, any constraint violation introduced by numerical error propagates away at the speed of light. This property increases the stability of the simulation, compared to System I above. A similar trick is used in the BSSN formulation of Einstein's equations.

In [2]:

# We inherit here all of the definitions from System I, above

# Step 7a: Register the scalar auxiliary variable \Gamma
Gamma = gri.register_gridfunctions("EVOL", ["Gamma"])

# Step 7b: Declare the ordinary gradient \partial_{i} \Gamma
Gamma_dD = ixp.declarerank1("Gamma_dD")

# Step 8a: Construct the second covariant derivative of the scalar \psi
# \psi_{;\hat{i}\hat{j}} = \psi_{,i;\hat{j}}
#                        = \psi_{,ij} - \Gamma^{k}_{ij} \psi_{,k}
psi_dDD = ixp.declarerank2("psi_dDD", "sym01")
psi_dHatDD = ixp.zerorank2()
for i in range(DIM):
    for j in range(DIM):
        psi_dHatDD[i][j] = psi_dDD[i][j]
        for k in range(DIM):
            psi_dHatDD[i][j] += - rfm.GammahatUDD[k][i][j] * psi_dD[k]

# Step 8b: Construct the covariant Laplacian of \psi
# Lappsi = ghat^{ij} D_{j} D_{i} \psi
Lappsi = 0
for i in range(DIM):
    for j in range(DIM):
        Lappsi += rfm.ghatUU[i][j] * psi_dHatDD[i][j]

# Step 9: Define right-hand sides for the evolution.
AD_rhs = ixp.zerorank1()
ED_rhs = ixp.zerorank1()
for i in range(DIM):
    AD_rhs[i] = -ED[i] - psi_dD[i]
    ED_rhs[i] = -LapAD[i] + Gamma_dD[i]
psi_rhs = -Gamma
Gamma_rhs = -Lappsi

# Step 10: Generate C code for System II Maxwell's evolution equations,
#          print output to the screen (standard out, or stdout).
lhrh_list = []
for i in range(DIM):
    lhrh_list.append(lhrh(lhs=gri.gfaccess("rhs_gfs", "AD" + str(i)), rhs=AD_rhs[i]))
    lhrh_list.append(lhrh(lhs=gri.gfaccess("rhs_gfs", "ED" + str(i)), rhs=ED_rhs[i]))
lhrh_list.append(lhrh(lhs=gri.gfaccess("rhs_gfs", "psi"), rhs=psi_rhs))
lhrh_list.append(lhrh(lhs=gri.gfaccess("rhs_gfs", "Gamma"), rhs=Gamma_rhs))

fin.FD_outputC("stdout", lhrh_list)

{
   /*
    * NRPy+ Finite Difference Code Generation, Step 1 of 2: Read from main memory and compute finite difference stencils:
    */
   /*
    *  Original SymPy expressions:
    *  "[const double AD_dD00 = invdx0*(-2*AD0_i0m1_i1_i2/3 + AD0_i0m2_i1_i2/12 + 2*AD0_i0p1_i1_i2/3 - AD0_i0p2_i1_i2/12),
    *    const double AD_dD01 = invdx1*(-2*AD0_i0_i1m1_i2/3 + AD0_i0_i1m2_i2/12 + 2*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),
    *    const double AD_dD02 = invdx2*(-2*AD0_i0_i1_i2m1/3 + AD0_i0_i1_i2m2/12 + 2*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),
    *    const double AD_dD10 = invdx0*(-2*AD1_i0m1_i1_i2/3 + AD1_i0m2_i1_i2/12 + 2*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),
    *    const double AD_dD11 = invdx1*(-2*AD1_i0_i1m1_i2/3 + AD1_i0_i1m2_i2/12 + 2*AD1_i0_i1p1_i2/3 - AD1_i0_i1p2_i2/12),
    *    const double AD_dD12 = invdx2*(-2*AD1_i0_i1_i2m1/3 + AD1_i0_i1_i2m2/12 + 2*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),
    *    const double AD_dD20 = invdx0*(-2*AD2_i0m1_i1_i2/3 + AD2_i0m2_i1_i2/12 + 2*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),
    *    const double AD_dD21 = invdx1*(-2*AD2_i0_i1m1_i2/3 + AD2_i0_i1m2_i2/12 + 2*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),
    *    const double AD_dD22 = invdx2*(-2*AD2_i0_i1_i2m1/3 + AD2_i0_i1_i2m2/12 + 2*AD2_i0_i1_i2p1/3 - AD2_i0_i1_i2p2/12),
    *    const double AD_dDD000 = invdx0**2*(-5*AD0/2 + 4*AD0_i0m1_i1_i2/3 - AD0_i0m2_i1_i2/12 + 4*AD0_i0p1_i1_i2/3 - AD0_i0p2_i1_i2/12),
    *    const double AD_dDD011 = invdx1**2*(-5*AD0/2 + 4*AD0_i0_i1m1_i2/3 - AD0_i0_i1m2_i2/12 + 4*AD0_i0_i1p1_i2/3 - AD0_i0_i1p2_i2/12),
    *    const double AD_dDD022 = invdx2**2*(-5*AD0/2 + 4*AD0_i0_i1_i2m1/3 - AD0_i0_i1_i2m2/12 + 4*AD0_i0_i1_i2p1/3 - AD0_i0_i1_i2p2/12),
    *    const double AD_dDD100 = invdx0**2*(-5*AD1/2 + 4*AD1_i0m1_i1_i2/3 - AD1_i0m2_i1_i2/12 + 4*AD1_i0p1_i1_i2/3 - AD1_i0p2_i1_i2/12),
    *    const double AD_dDD111 = invdx1**2*(-5*AD1/2 + 4*AD1_i0_i1m1_i2/3 - AD1_i0_i1m2_i2/12 + 4*AD1_i0_i1p1_i2/3 - AD1_i0_i1p2_i2/12),
    *    const double AD_dDD122 = invdx2**2*(-5*AD1/2 + 4*AD1_i0_i1_i2m1/3 - AD1_i0_i1_i2m2/12 + 4*AD1_i0_i1_i2p1/3 - AD1_i0_i1_i2p2/12),
    *    const double AD_dDD200 = invdx0**2*(-5*AD2/2 + 4*AD2_i0m1_i1_i2/3 - AD2_i0m2_i1_i2/12 + 4*AD2_i0p1_i1_i2/3 - AD2_i0p2_i1_i2/12),
    *    const double AD_dDD211 = invdx1**2*(-5*AD2/2 + 4*AD2_i0_i1m1_i2/3 - AD2_i0_i1m2_i2/12 + 4*AD2_i0_i1p1_i2/3 - AD2_i0_i1p2_i2/12),
    *    const double AD_dDD222 = invdx2**2*(-5*AD2/2 + 4*AD2_i0_i1_i2m1/3 - AD2_i0_i1_i2m2/12 + 4*AD2_i0_i1_i2p1/3 - AD2_i0_i1_i2p2/12),
    *    const double Gamma_dD0 = invdx0*(-2*Gamma_i0m1_i1_i2/3 + Gamma_i0m2_i1_i2/12 + 2*Gamma_i0p1_i1_i2/3 - Gamma_i0p2_i1_i2/12),
    *    const double Gamma_dD1 = invdx1*(-2*Gamma_i0_i1m1_i2/3 + Gamma_i0_i1m2_i2/12 + 2*Gamma_i0_i1p1_i2/3 - Gamma_i0_i1p2_i2/12),
    *    const double Gamma_dD2 = invdx2*(-2*Gamma_i0_i1_i2m1/3 + Gamma_i0_i1_i2m2/12 + 2*Gamma_i0_i1_i2p1/3 - Gamma_i0_i1_i2p2/12),
    *    const double psi_dD0 = invdx0*(-2*psi_i0m1_i1_i2/3 + psi_i0m2_i1_i2/12 + 2*psi_i0p1_i1_i2/3 - psi_i0p2_i1_i2/12),
    *    const double psi_dD1 = invdx1*(-2*psi_i0_i1m1_i2/3 + psi_i0_i1m2_i2/12 + 2*psi_i0_i1p1_i2/3 - psi_i0_i1p2_i2/12),
    *    const double psi_dD2 = invdx2*(-2*psi_i0_i1_i2m1/3 + psi_i0_i1_i2m2/12 + 2*psi_i0_i1_i2p1/3 - psi_i0_i1_i2p2/12),
    *    const double psi_dDD00 = invdx0**2*(-5*psi/2 + 4*psi_i0m1_i1_i2/3 - psi_i0m2_i1_i2/12 + 4*psi_i0p1_i1_i2/3 - psi_i0p2_i1_i2/12),
    *    const double psi_dDD11 = invdx1**2*(-5*psi/2 + 4*psi_i0_i1m1_i2/3 - psi_i0_i1m2_i2/12 + 4*psi_i0_i1p1_i2/3 - psi_i0_i1p2_i2/12),
    *    const double psi_dDD22 = invdx2**2*(-5*psi/2 + 4*psi_i0_i1_i2m1/3 - psi_i0_i1_i2m2/12 + 4*psi_i0_i1_i2p1/3 - psi_i0_i1_i2p2/12)]"
    */
   const double psi_i0_i1_i2m2 = in_gfs[IDX4(PSIGF, i0,i1,i2-2)];
   const double psi_i0_i1_i2m1 = in_gfs[IDX4(PSIGF, i0,i1,i2-1)];
   const double psi_i0_i1m2_i2 = in_gfs[IDX4(PSIGF, i0,i1-2,i2)];
   const double psi_i0_i1m1_i2 = in_gfs[IDX4(PSIGF, i0,i1-1,i2)];
   const double psi_i0m2_i1_i2 = in_gfs[IDX4(PSIGF, i0-2,i1,i2)];
   const double psi_i0m1_i1_i2 = in_gfs[IDX4(PSIGF, i0-1,i1,i2)];
   const double psi = in_gfs[IDX4(PSIGF, i0,i1,i2)];
   const double psi_i0p1_i1_i2 = in_gfs[IDX4(PSIGF, i0+1,i1,i2)];
   const double psi_i0p2_i1_i2 = in_gfs[IDX4(PSIGF, i0+2,i1,i2)];
   const double psi_i0_i1p1_i2 = in_gfs[IDX4(PSIGF, i0,i1+1,i2)];
   const double psi_i0_i1p2_i2 = in_gfs[IDX4(PSIGF, i0,i1+2,i2)];
   const double psi_i0_i1_i2p1 = in_gfs[IDX4(PSIGF, i0,i1,i2+1)];
   const double psi_i0_i1_i2p2 = in_gfs[IDX4(PSIGF, i0,i1,i2+2)];
   const double ED0 = in_gfs[IDX4(ED0GF, i0,i1,i2)];
   const double ED1 = in_gfs[IDX4(ED1GF, i0,i1,i2)];
   const double ED2 = in_gfs[IDX4(ED2GF, i0,i1,i2)];
   const double AD0_i0_i1_i2m2 = in_gfs[IDX4(AD0GF, i0,i1,i2-2)];
   const double AD0_i0_i1_i2m1 = in_gfs[IDX4(AD0GF, i0,i1,i2-1)];
   const double AD0_i0_i1m2_i2 = in_gfs[IDX4(AD0GF, i0,i1-2,i2)];
   const double AD0_i0_i1m1_i2 = in_gfs[IDX4(AD0GF, i0,i1-1,i2)];
   const double AD0_i0m2_i1_i2 = in_gfs[IDX4(AD0GF, i0-2,i1,i2)];
   const double AD0_i0m1_i1_i2 = in_gfs[IDX4(AD0GF, i0-1,i1,i2)];
   const double AD0 = in_gfs[IDX4(AD0GF, i0,i1,i2)];
   const double AD0_i0p1_i1_i2 = in_gfs[IDX4(AD0GF, i0+1,i1,i2)];
   const double AD0_i0p2_i1_i2 = in_gfs[IDX4(AD0GF, i0+2,i1,i2)];
   const double AD0_i0_i1p1_i2 = in_gfs[IDX4(AD0GF, i0,i1+1,i2)];
   const double AD0_i0_i1p2_i2 = in_gfs[IDX4(AD0GF, i0,i1+2,i2)];
   const double AD0_i0_i1_i2p1 = in_gfs[IDX4(AD0GF, i0,i1,i2+1)];
   const double AD0_i0_i1_i2p2 = in_gfs[IDX4(AD0GF, i0,i1,i2+2)];
   const double AD1_i0_i1_i2m2 = in_gfs[IDX4(AD1GF, i0,i1,i2-2)];
   const double AD1_i0_i1_i2m1 = in_gfs[IDX4(AD1GF, i0,i1,i2-1)];
   const double AD1_i0_i1m2_i2 = in_gfs[IDX4(AD1GF, i0,i1-2,i2)];
   const double AD1_i0_i1m1_i2 = in_gfs[IDX4(AD1GF, i0,i1-1,i2)];
   const double AD1_i0m2_i1_i2 = in_gfs[IDX4(AD1GF, i0-2,i1,i2)];
   const double AD1_i0m1_i1_i2 = in_gfs[IDX4(AD1GF, i0-1,i1,i2)];
   const double AD1 = in_gfs[IDX4(AD1GF, i0,i1,i2)];
   const double AD1_i0p1_i1_i2 = in_gfs[IDX4(AD1GF, i0+1,i1,i2)];
   const double AD1_i0p2_i1_i2 = in_gfs[IDX4(AD1GF, i0+2,i1,i2)];
   const double AD1_i0_i1p1_i2 = in_gfs[IDX4(AD1GF, i0,i1+1,i2)];
   const double AD1_i0_i1p2_i2 = in_gfs[IDX4(AD1GF, i0,i1+2,i2)];
   const double AD1_i0_i1_i2p1 = in_gfs[IDX4(AD1GF, i0,i1,i2+1)];
   const double AD1_i0_i1_i2p2 = in_gfs[IDX4(AD1GF, i0,i1,i2+2)];
   const double AD2_i0_i1_i2m2 = in_gfs[IDX4(AD2GF, i0,i1,i2-2)];
   const double AD2_i0_i1_i2m1 = in_gfs[IDX4(AD2GF, i0,i1,i2-1)];
   const double AD2_i0_i1m2_i2 = in_gfs[IDX4(AD2GF, i0,i1-2,i2)];
   const double AD2_i0_i1m1_i2 = in_gfs[IDX4(AD2GF, i0,i1-1,i2)];
   const double AD2_i0m2_i1_i2 = in_gfs[IDX4(AD2GF, i0-2,i1,i2)];
   const double AD2_i0m1_i1_i2 = in_gfs[IDX4(AD2GF, i0-1,i1,i2)];
   const double AD2 = in_gfs[IDX4(AD2GF, i0,i1,i2)];
   const double AD2_i0p1_i1_i2 = in_gfs[IDX4(AD2GF, i0+1,i1,i2)];
   const double AD2_i0p2_i1_i2 = in_gfs[IDX4(AD2GF, i0+2,i1,i2)];
   const double AD2_i0_i1p1_i2 = in_gfs[IDX4(AD2GF, i0,i1+1,i2)];
   const double AD2_i0_i1p2_i2 = in_gfs[IDX4(AD2GF, i0,i1+2,i2)];
   const double AD2_i0_i1_i2p1 = in_gfs[IDX4(AD2GF, i0,i1,i2+1)];
   const double AD2_i0_i1_i2p2 = in_gfs[IDX4(AD2GF, i0,i1,i2+2)];
   const double Gamma_i0_i1_i2m2 = in_gfs[IDX4(GAMMAGF, i0,i1,i2-2)];
   const double Gamma_i0_i1_i2m1 = in_gfs[IDX4(GAMMAGF, i0,i1,i2-1)];
   const double Gamma_i0_i1m2_i2 = in_gfs[IDX4(GAMMAGF, i0,i1-2,i2)];
   const double Gamma_i0_i1m1_i2 = in_gfs[IDX4(GAMMAGF, i0,i1-1,i2)];
   const double Gamma_i0m2_i1_i2 = in_gfs[IDX4(GAMMAGF, i0-2,i1,i2)];
   const double Gamma_i0m1_i1_i2 = in_gfs[IDX4(GAMMAGF, i0-1,i1,i2)];
   const double Gamma = in_gfs[IDX4(GAMMAGF, i0,i1,i2)];
   const double Gamma_i0p1_i1_i2 = in_gfs[IDX4(GAMMAGF, i0+1,i1,i2)];
   const double Gamma_i0p2_i1_i2 = in_gfs[IDX4(GAMMAGF, i0+2,i1,i2)];
   const double Gamma_i0_i1p1_i2 = in_gfs[IDX4(GAMMAGF, i0,i1+1,i2)];
   const double Gamma_i0_i1p2_i2 = in_gfs[IDX4(GAMMAGF, i0,i1+2,i2)];
   const double Gamma_i0_i1_i2p1 = in_gfs[IDX4(GAMMAGF, i0,i1,i2+1)];
   const double Gamma_i0_i1_i2p2 = in_gfs[IDX4(GAMMAGF, i0,i1,i2+2)];
   const double FDPart1_Rational_2_3 = 2.0/3.0;
   const double FDPart1_Rational_1_12 = 1.0/12.0;
   const double FDPart1_Rational_5_2 = 5.0/2.0;
   const double FDPart1_Rational_4_3 = 4.0/3.0;
   const double FDPart1_1 = -AD0_i0_i1p2_i2;
   const double FDPart1_9 = ((invdx0)*(invdx0));
   const double FDPart1_10 = -AD0*FDPart1_Rational_5_2;
   const double FDPart1_11 = ((invdx1)*(invdx1));
   const double FDPart1_12 = ((invdx2)*(invdx2));
   const double FDPart1_13 = -AD1*FDPart1_Rational_5_2;
   const double FDPart1_14 = -AD2*FDPart1_Rational_5_2;
   const double FDPart1_18 = -FDPart1_Rational_5_2*psi;
   const double AD_dD00 = invdx0*(FDPart1_Rational_1_12*(AD0_i0m2_i1_i2 - AD0_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD0_i0m1_i1_i2 + AD0_i0p1_i1_i2));
   const double AD_dD01 = invdx1*(FDPart1_Rational_1_12*(AD0_i0_i1m2_i2 + FDPart1_1) + FDPart1_Rational_2_3*(-AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));
   const double AD_dD02 = invdx2*(FDPart1_Rational_1_12*(AD0_i0_i1_i2m2 - AD0_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));
   const double AD_dD10 = invdx0*(FDPart1_Rational_1_12*(AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));
   const double AD_dD11 = invdx1*(FDPart1_Rational_1_12*(AD1_i0_i1m2_i2 - AD1_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD1_i0_i1m1_i2 + AD1_i0_i1p1_i2));
   const double AD_dD12 = invdx2*(FDPart1_Rational_1_12*(AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));
   const double AD_dD20 = invdx0*(FDPart1_Rational_1_12*(AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_2_3*(-AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));
   const double AD_dD21 = invdx1*(FDPart1_Rational_1_12*(AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_2_3*(-AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));
   const double AD_dD22 = invdx2*(FDPart1_Rational_1_12*(AD2_i0_i1_i2m2 - AD2_i0_i1_i2p2) + FDPart1_Rational_2_3*(-AD2_i0_i1_i2m1 + AD2_i0_i1_i2p1));
   const double AD_dDD000 = FDPart1_9*(FDPart1_10 + FDPart1_Rational_1_12*(-AD0_i0m2_i1_i2 - AD0_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD0_i0m1_i1_i2 + AD0_i0p1_i1_i2));
   const double AD_dDD011 = FDPart1_11*(FDPart1_10 + FDPart1_Rational_1_12*(-AD0_i0_i1m2_i2 + FDPart1_1) + FDPart1_Rational_4_3*(AD0_i0_i1m1_i2 + AD0_i0_i1p1_i2));
   const double AD_dDD022 = FDPart1_12*(FDPart1_10 + FDPart1_Rational_1_12*(-AD0_i0_i1_i2m2 - AD0_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD0_i0_i1_i2m1 + AD0_i0_i1_i2p1));
   const double AD_dDD100 = FDPart1_9*(FDPart1_13 + FDPart1_Rational_1_12*(-AD1_i0m2_i1_i2 - AD1_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD1_i0m1_i1_i2 + AD1_i0p1_i1_i2));
   const double AD_dDD111 = FDPart1_11*(FDPart1_13 + FDPart1_Rational_1_12*(-AD1_i0_i1m2_i2 - AD1_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD1_i0_i1m1_i2 + AD1_i0_i1p1_i2));
   const double AD_dDD122 = FDPart1_12*(FDPart1_13 + FDPart1_Rational_1_12*(-AD1_i0_i1_i2m2 - AD1_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD1_i0_i1_i2m1 + AD1_i0_i1_i2p1));
   const double AD_dDD200 = FDPart1_9*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0m2_i1_i2 - AD2_i0p2_i1_i2) + FDPart1_Rational_4_3*(AD2_i0m1_i1_i2 + AD2_i0p1_i1_i2));
   const double AD_dDD211 = FDPart1_11*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0_i1m2_i2 - AD2_i0_i1p2_i2) + FDPart1_Rational_4_3*(AD2_i0_i1m1_i2 + AD2_i0_i1p1_i2));
   const double AD_dDD222 = FDPart1_12*(FDPart1_14 + FDPart1_Rational_1_12*(-AD2_i0_i1_i2m2 - AD2_i0_i1_i2p2) + FDPart1_Rational_4_3*(AD2_i0_i1_i2m1 + AD2_i0_i1_i2p1));
   const double Gamma_dD0 = invdx0*(FDPart1_Rational_1_12*(Gamma_i0m2_i1_i2 - Gamma_i0p2_i1_i2) + FDPart1_Rational_2_3*(-Gamma_i0m1_i1_i2 + Gamma_i0p1_i1_i2));
   const double Gamma_dD1 = invdx1*(FDPart1_Rational_1_12*(Gamma_i0_i1m2_i2 - Gamma_i0_i1p2_i2) + FDPart1_Rational_2_3*(-Gamma_i0_i1m1_i2 + Gamma_i0_i1p1_i2));
   const double Gamma_dD2 = invdx2*(FDPart1_Rational_1_12*(Gamma_i0_i1_i2m2 - Gamma_i0_i1_i2p2) + FDPart1_Rational_2_3*(-Gamma_i0_i1_i2m1 + Gamma_i0_i1_i2p1));
   const double psi_dD0 = invdx0*(FDPart1_Rational_1_12*(psi_i0m2_i1_i2 - psi_i0p2_i1_i2) + FDPart1_Rational_2_3*(-psi_i0m1_i1_i2 + psi_i0p1_i1_i2));
   const double psi_dD1 = invdx1*(FDPart1_Rational_1_12*(psi_i0_i1m2_i2 - psi_i0_i1p2_i2) + FDPart1_Rational_2_3*(-psi_i0_i1m1_i2 + psi_i0_i1p1_i2));
   const double psi_dD2 = invdx2*(FDPart1_Rational_1_12*(psi_i0_i1_i2m2 - psi_i0_i1_i2p2) + FDPart1_Rational_2_3*(-psi_i0_i1_i2m1 + psi_i0_i1_i2p1));
   const double psi_dDD00 = FDPart1_9*(FDPart1_18 + FDPart1_Rational_1_12*(-psi_i0m2_i1_i2 - psi_i0p2_i1_i2) + FDPart1_Rational_4_3*(psi_i0m1_i1_i2 + psi_i0p1_i1_i2));
   const double psi_dDD11 = FDPart1_11*(FDPart1_18 + FDPart1_Rational_1_12*(-psi_i0_i1m2_i2 - psi_i0_i1p2_i2) + FDPart1_Rational_4_3*(psi_i0_i1m1_i2 + psi_i0_i1p1_i2));
   const double psi_dDD22 = FDPart1_12*(FDPart1_18 + FDPart1_Rational_1_12*(-psi_i0_i1_i2m2 - psi_i0_i1_i2p2) + FDPart1_Rational_4_3*(psi_i0_i1_i2m1 + psi_i0_i1_i2p1));
   /*
    * NRPy+ Finite Difference Code Generation, Step 2 of 2: Evaluate SymPy expressions and write to main memory:
    */
   /*
    *  Original SymPy expressions:
    *  "[rhs_gfs[IDX4(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0,
    *    rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = -AD_dDD000 + Gamma_dD0 - (AD_dD00*xx0 - AD_dD11/xx0 + AD_dDD011 - (AD0*xx0 + AD_dD11)/xx0)/xx0**2 - (AD_dD00*xx0*sin(xx1)**2 - AD_dD22/xx0 + AD_dDD022 + (-AD1/xx0 + AD_dD01)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)/xx0)/(xx0**2*sin(xx1)**2),
    *    rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1,
    *    rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1/xx0**2 + AD_dD10/xx0 - AD_dDD100 + Gamma_dD1 + (-AD1/xx0 + AD_dD10)/xx0 - (AD_dD01*xx0 + AD_dDD111 + xx0*(-AD1/xx0 + AD_dD01) + xx0*(-AD1/xx0 + AD_dD10))/xx0**2 - (-AD_dD22*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD122 + xx0*(-AD1/xx0 + AD_dD10)*sin(xx1)**2 + (AD0*xx0 + AD_dD11)*sin(2*xx1)/2 - (AD0*xx0*sin(xx1)**2 + AD1*sin(2*xx1)/2 + AD_dD22)*sin(2*xx1)/(2*sin(xx1)**2))/(xx0**2*sin(xx1)**2),
    *    rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2,
    *    rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2/xx0**2 + AD_dD20/xx0 - AD_dDD200 + Gamma_dD2 + (-AD2/xx0 + AD_dD20)/xx0 - (AD2*(-cos(2*xx1)/sin(xx1)**2 + sin(2*xx1)*cos(xx1)/sin(xx1)**3) - AD_dD21*sin(2*xx1)/(2*sin(xx1)**2) + AD_dDD211 + xx0*(-AD2/xx0 + AD_dD20) - (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/(2*sin(xx1)**2))/xx0**2 - (AD_dD02*xx0*sin(xx1)**2 + AD_dD12*sin(2*xx1)/2 + AD_dDD222 + xx0*(-AD2/xx0 + AD_dD02)*sin(xx1)**2 + xx0*(-AD2/xx0 + AD_dD20)*sin(xx1)**2 + (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD12)*sin(2*xx1)/2 + (-AD2*sin(2*xx1)/(2*sin(xx1)**2) + AD_dD21)*sin(2*xx1)/2)/(xx0**2*sin(xx1)**2),
    *    rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -Gamma,
    *    rhs_gfs[IDX4(GAMMAGF, i0, i1, i2)] = -psi_dDD00 - (psi_dD0*xx0 + psi_dDD11)/xx0**2 - (psi_dD0*xx0*sin(xx1)**2 + psi_dD1*sin(2*xx1)/2 + psi_dDD22)/(xx0**2*sin(xx1)**2)]"
    */
   const double FDPart3_0 = (1.0/((xx0)*(xx0)));
   const double FDPart3_2 = (1.0/(xx0));
   const double FDPart3_4 = AD0*xx0 + AD_dD11;
   const double FDPart3_5 = sin(xx1);
   const double FDPart3_6 = ((FDPart3_5)*(FDPart3_5));
   const double FDPart3_7 = -AD1*FDPart3_2;
   const double FDPart3_10 = sin(2*xx1);
   const double FDPart3_11 = (1.0/2.0)*FDPart3_10;
   const double FDPart3_12 = AD0*FDPart3_6*xx0 + AD1*FDPart3_11 + AD_dD22;
   const double FDPart3_13 = (1.0/(FDPart3_6));
   const double FDPart3_14 = FDPart3_0*FDPart3_13;
   const double FDPart3_16 = xx0*(AD_dD10 + FDPart3_7);
   const double FDPart3_17 = FDPart3_11*FDPart3_13;
   const double FDPart3_18 = -AD2*FDPart3_2;
   const double FDPart3_20 = xx0*(AD_dD20 + FDPart3_18);
   const double FDPart3_21 = -AD2*FDPart3_17;
   const double FDPart3_22 = FDPart3_11*(AD_dD21 + FDPart3_21);
   rhs_gfs[IDX4(AD0GF, i0, i1, i2)] = -ED0 - psi_dD0;
   rhs_gfs[IDX4(ED0GF, i0, i1, i2)] = -AD_dDD000 - FDPart3_0*(AD_dD00*xx0 - AD_dD11*FDPart3_2 + AD_dDD011 - FDPart3_2*FDPart3_4) - FDPart3_14*(AD_dD00*FDPart3_6*xx0 - AD_dD22*FDPart3_2 + AD_dDD022 + FDPart3_11*(AD_dD01 + FDPart3_7) - FDPart3_12*FDPart3_2) + Gamma_dD0;
   rhs_gfs[IDX4(AD1GF, i0, i1, i2)] = -ED1 - psi_dD1;
   rhs_gfs[IDX4(ED1GF, i0, i1, i2)] = -AD1*FDPart3_0 + AD_dD10*FDPart3_2 - AD_dDD100 - FDPart3_0*(AD_dD01*xx0 + AD_dDD111 + FDPart3_16 + xx0*(AD_dD01 + FDPart3_7)) - FDPart3_14*(-AD_dD22*FDPart3_17 + AD_dDD122 + FDPart3_11*FDPart3_4 - FDPart3_12*FDPart3_17 + FDPart3_16*FDPart3_6) + FDPart3_2*(AD_dD10 + FDPart3_7) + Gamma_dD1;
   rhs_gfs[IDX4(AD2GF, i0, i1, i2)] = -ED2 - psi_dD2;
   rhs_gfs[IDX4(ED2GF, i0, i1, i2)] = -AD2*FDPart3_0 + AD_dD20*FDPart3_2 - AD_dDD200 - FDPart3_0*(AD2*(FDPart3_10*cos(xx1)/((FDPart3_5)*(FDPart3_5)*(FDPart3_5)) - FDPart3_13*cos(2*xx1)) - AD_dD21*FDPart3_17 + AD_dDD211 - FDPart3_13*FDPart3_22 + FDPart3_20) - FDPart3_14*(AD_dD02*FDPart3_6*xx0 + AD_dD12*FDPart3_11 + AD_dDD222 + FDPart3_11*(AD_dD12 + FDPart3_21) + FDPart3_20*FDPart3_6 + FDPart3_22 + FDPart3_6*xx0*(AD_dD02 + FDPart3_18)) + FDPart3_2*(AD_dD20 + FDPart3_18) + Gamma_dD2;
   rhs_gfs[IDX4(PSIGF, i0, i1, i2)] = -Gamma;
   rhs_gfs[IDX4(GAMMAGF, i0, i1, i2)] = -FDPart3_0*(psi_dD0*xx0 + psi_dDD11) - FDPart3_14*(FDPart3_11*psi_dD1 + FDPart3_6*psi_dD0*xx0 + psi_dDD22) - psi_dDD00;
}

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

In [3]:

import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface
cmd.output_Jupyter_notebook_to_LaTeXed_PDF("Tutorial-MaxwellCurvilinear")

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

Generating C code for the right-hand sides of Maxwell's equations, in *curvilinear* coordinates, using a reference metric formalism¶