Notebook

Using Indexed Expressions Efficiently: Examples and Exercises¶

Author: Patrick Nelson¶

Formatting improvements courtesy Brandon Clark¶

NRPy+ Source Code for this module:¶

TODO?

Introduction:¶

NRPy+ provides several methods to declare variables and lists to represent the tensors, vectors, and scalars of general relativity. Each function has its use case, but when starting out, it is not always trivial to determine which is appropriate at any given time. This is further complicated by the fact that these use cases can change depending on whether or not the notebook at hand is going to output C code or not. So, this module will provide some pedagogy to help new users with these functions.

We will first explore the case in which one does not intend to immediately output C code. This is done in tutorials that solely generate symbolic expressions and their corresponding modules. For example, consider the tutorial Tutorial-GRHD_Equations-Cartesian.ipynb and the corresponding module. The tutorial notebook gives in-depth, $\latex$ documentation about the GRHD expressions interspersed with the python code that generates the corresponding sympy expressions. The module provides the same python code, but in a format that makes it easy to import it into other Jupyter notebooks and modules. So, the tutorial also performs a self-validation check to guarantee that the notebook and module are generating identical sympy expressions.

Then, we will explore the case in which one is outputting C code. In this case, one will import the python modules described above to generate any symbolic expressions needed. These expressions can then be passed to FD_outputC(). But this function requires that every symbol in the expressions passed to it must be either a gridfunction or C parameter to help make sure that the resulting C code compiles correctly.

This tutorial assumes that the reader has previously looked over the previous module about indexed expressions in NRPy+, available here. If you have not read this, please do so before continuing to ensure that you are familiar with the basic syntax for the functions provided by indexedexp.py. This module is focused more on how and when to use each of those functions.

Table of Contents¶

$$\label{toc}$$

This notebook is organized as follows

Step 1: Zero-rank, Declare-rank, Symbols
1. Step 1.a Exercises
Step 2: Gridfunctions and C parameters
1. Step 2.a Exercises
Step 3: Applications
Step 4 Solutions for Exercises in Steps 2 & 3

Step 1: Zero-rank, Declare-rank, Symbols [Back to top]¶

$$\label{no_c}$$

If we do not intend for a notebook to output C code, we only need to consider two families. In the first, we declare something symbolically; in the second, we set it to zero. For example, consider a simple index-lowering operation of a three-vector $A^i$ using the three-metric $\gamma_{ij}$: $$ A_i = \gamma_{ij} A^j. $$ We've written this in a way that suggests we will define new quantities, $A_i$, in terms of a known quantities, $A^i$ and $\gamma_{ij}$. The known quantities must be declared; depending on the rank of the quantity in question, this will use one of the following functions:

sympy.symbols()
indexedexp.declarerank1()
indexedexp.declarerank2()
indexedexp.declarerank3()
indexedexp.declarerank4()

The parameters that these take are detailed in a previous tutorial. The new quantities belong to the other family. These must be zeroed before we can add things to them; this is analogous to initializing a variable in C. This is done with one of the following functions, depending on the rank of the quantities:

sympy.sympify(0)
indexedexp.zerorank1()
indexedexp.zerorank2()
indexedexp.zerorank3()
indexedexp.zerorank4()

These all generate either zero or lists of zero that are compatible with sympify.

Let us again consider the example $$ A_i = \gamma_{ij} A^j. $$

In [3]:

import indexedexp as ixp         # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support
import sympy as sp

# First we'll declare the knowns:
AU = ixp.declarerank1("AU",DIM=3)
gammaDD = ixp.declarerank2("gammaDD","sym01",DIM=3) # The metric tensor is symmetric, so we set it as such.

# Now, let's initialize the new quantity:
AD = ixp.zerorank1()

# Finally, we can loop over all the indices in the implied sum to express the new
# quantity in terms of the knowns:
for i in range(3):
    for j in range(3):
        AD[i] += gammaDD[i][j] * AU[j]

Step 1.a: Exercises [Back to top]¶

$$\label{no_c_ex}$$

Code the following expressions, paying particular attention to the initial declarations and zeroings of tensors and scalars.

$$ a = \gamma_{ij} A^i A^j $$

In [ ]:

# What are we trying to compute? What known quantities is that in terms of?
# Declare the knowns:

# And initialize the new quantity you want to compute to zero:

# Then loop through the indices to build the quantity:

$$ T^\xi_{\ \ \mu\nu} = k L^\xi M_\mu N_\nu $$

In [ ]:

# Remember that Greek letters denote four-dimensional spacetime quantities!

Step 2: Gridfunctions and C parameters [Back to top]¶

$$\label{gf_c}$$

Once we are ready to write a start-to-finish module, in which we generate, compile, and then run C code, we will also need to start registering gridfunctions and C parameters. Any expression passed to FD_outputC() must be entirely in terms of gridfunctions and C parameters; this requirement was put in place to help make sure that the automatically generated files cover everything that they need to do in order to minimize the amount of hand-coding in C that will need to be done. So, we will need to use the following functions:

NRPy_param_funcs.Cparameters()
grid.register_gridfunctions()
indexedexp.register_gridfunctions_for_single_rank1()
indexedexp.register_gridfunctions_for_single_rank2()

The zero-rank and declare-rank functions still have uses here, though. We will use declarerank whenever we want to take a finite-difference derivative; by appropriately naming it, NRPy+ will automatically generate code to differentiate a specified gridfunction as detailed in its tutorial. The zerorank functions find a use for intermediate expressions and for expressions whose variable name does not match the gridfunction's name (e.g. the right-hand side of an evolution equation. We will demonstrate these points with a mock PDE. Suppose that there is some vector $J_i$ that evolves according to the following equation: $$ \partial_t J_i = \partial_j J^j K_i - G \Lambda_{i}, $$

where $G$ is Newton's gravitational constant and $\Lambda_{i} = J_i K_j K^j$. Suppose that, in this simulation, the quantities we want to evolve (and thus store) are $J_i$ and $K^i$

(N.B. This equation is not intended to correspond to any particular physics, merely to demonstrate many different use cases of NRPy+'s indexed expressions utilities.)

In [6]:

import NRPy_param_funcs as par   # NRPy+: Parameter interface
import grid as gri               # NRPy+: Functions having to do with numerical grids

Step 2.a: Exercises [Back to top]¶

$$\label{gf_c_ex}$$

In [ ]:

Step 3: Applications [Back to top]¶

$$\label{applications}$$

Step 4: Solutions for Exercises in Steps 2 & 3 [Back to top]¶

$$\label{key}$$$$ a = \gamma_{ij} A^i A^j $$

In [4]:

gammaDD = ixp.declarerank2("gammaDD","sym01",DIM=3)
AU = ixp.declarerank1("AU",DIM=3)

a = sp.sympify(0)

for i in range(3):
    for j in range(3):
        a += gammaDD[i][j] * AU[i] * AU[j]

$$ T^\xi_{\ \ \mu\nu} = k L^\xi M_\mu N_\nu $$

In [ ]:

k = sp.symbols("k",real=True)
LU = ixp.declarerank1("LU",DIM=4)
MD = ixp.declarerank1("MD",DIM=4)
ND = ixp.declarerank1("ND",DIM=4)

TUDD = ixp.zerorank3()

for xi in range(4):
    for mu in range(4):
        for nu in range(4):
            TUDD[xi][mu][nu] += k * LU[xi] * MD[mu] * ND[nu]