#!/usr/bin/env python
# coding: utf-8

# <script async src="https://www.googletagmanager.com/gtag/js?id=UA-59152712-8"></script>
# <script>
#   window.dataLayer = window.dataLayer || [];
#   function gtag(){dataLayer.push(arguments);}
#   gtag('js', new Date());
# 
#   gtag('config', 'UA-59152712-8');
# </script>
# 
# # NRPy+ 10-Minute Overview
# 
# ## Author: Zach Etienne
# 
# ## This notebook presents NRPy+, a Python framework that converts intricate SymPy expressions into optimized C/C++ code for use in numerical relativity. It emphasizes the computation of 3-Christoffel symbols $\Gamma^{i}_{jk}$ and numerical derivative designations within expressions. 
# [comment]: <> (added this abstract in case. omit of it's unnecessary.)

# <a id='toc'></a>
# 
# # Table of Contents
# $$\label{toc}$$ 
# 
# This notebook is organized as follows
# 
# 1. [Part 1](#why): Why NRPy+?
# 1. [Part 2](#christoffel_symbols): Constructing 3-Christoffels $\Gamma^i_{jk}$ as symbolic expressions in terms of 3-metric $\gamma_{ij}$ and derivatives
# 1. [Part 3](#outputc): `outputC()` example: Output $\Gamma^i_{jk}$ expressions as optimized C code, assuming derivatives already specified
# 1. [Part 4](#fd_outputc): `FD_outputC()` example: Specify numerical derivatives within $\Gamma^i_{jk}$ expressions as finite differences
# 1. [Part 5](#what_next): What next? Navigating the NRPy+ tutorial
# 1. [Part 6](#latex_pdf_output): Output this notebook to $\LaTeX$-formatted PDF    

# <a id='why'></a>
# 
# # Part 1: Why NRPy+? \[Back to [top](#toc)\]
# $$\label{why}$$
# 
# A core problem faced by numerical relativity is not that the techniques we use are particularly complex, but that *the equations numerical relativists solve are very complex*.
# 
# [Einstein notation](https://en.wikipedia.org/wiki/Einstein_notation) certainly makes the complexity more manageable, and is the most common way to perform tensorial mathematics in numerical relativity.
# 
# NRPy+ is built upon the idea that equations written in Einstein notation *are a form of code*, and when this code is expressed in NRPy+'s native (Python) language
# 
# * rank-0 tensors (scalars) are *variables*
# * rank-1 tensors are *lists*
# * rank-2 tensors are *lists of lists*
# * ... and so forth.
# 
# Further, implied tensor summations are simply *loops*.
# 
# NRPy+ combines the above ideas with the incredibly powerful [SymPy](https://www.sympy.org/) computer algebra package (think: Mathematica for Python, but fully free & open source), with a custom-built code generation infrastructure for converting complex SymPy expressions directly into highly-optimized C/C++ code.
# 
# Importantly, NRPy+ builds on the idea of *learning by example*. The core NRPy+ repository contains more than 100 pedagogical and well-formatted Jupyter notebook tutorials, covering topics of core relevance to the field of numerical relativity. About 25 of these tutorials generate *complete C codes* capable of e.g., evolving the scalar wave equation, Maxwell's equations, and Einstein's equations of general relativity -- all in a variety of coordinate systems. All ~100 Jupyter notebooks are linked to from [the main NRPy+ tutorial table of contents.](NRPyPlus_Tutorial.ipynb)
# 
# This 10-minute overview however is designed to introduce only the very basic features of NRPy+, with a core focus on the idea that *NRPy+ can be used to benefit any numerical relativity code*.

# <a id='christoffel_symbols'></a>
# 
# # Part 2: Constructing 3-Christoffels $\Gamma^i_{jk}$ as symbolic expressions in terms of 3-metric $\gamma_{ij}$ and derivatives \[Back to [top](#toc)\]
# $$\label{christoffel_symbols}$$
# 
# **Problem statement**: Given a three-metric $\gamma_{ij}$, construct all 18 independent Christoffel symbols $\Gamma^i_{jk}$, which involves first derivatives of the metric. Assume that $\gamma_{ij}$ *and its derivatives* are given numerically, requiring the derivatives to be defined as symbols.
# 
# In NRPy+ we adopt a rigid syntax for tensors and indexed expressions involving Python lists, so that for example
# 
# * $\gamma_{ij}=$ `gammaDD[i][j]`
# * $\gamma_{ij,k}=$ `gammaDD_dD[i][j][k]`
# 
# Christoffel symbols (of the first kind) are defined as ([source](https://en.wikipedia.org/wiki/Christoffel_symbols#Christoffel_symbols_of_the_first_kind)):
# 
# \begin{align}
# \Gamma_{ij}^k &= \frac{1}{2} \gamma^{kl}\left(\gamma_{jl,i} + \gamma_{il,j} - \gamma_{ij,l}\right)\\
# \end{align}
# 
# So first we'll define $\gamma_{ij}$ and its inverse using NRPy+ functions

# In[1]:


import sympy as sp        # SymPy: The Python computer algebra package upon which NRPy+ depends
import indexedexp as ixp  # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support

# gammaDD is a rank-2 indexed expression symmetric in indexes 0 and 1
gammaDD = ixp.declarerank2("gammaDD", symmetry="sym01", DIM=3)

# gammaUU is the inverse of gammaDD, and
# gammaDET (unused) is the determinant of gammaDD
gammaUU, gammaDET = ixp.symm_matrix_inverter3x3(gammaDD)  # inverts a 3x3 symmetric matrix


# Below we interact a bit with these generated expressions, confirming they 
# 
# *Exercise to students*: Verify that $\gamma_{ij} \gamma^{ij} = 3$ using the provided data structures & Python loops. (*Hint 1*: Scroll down a couple of cells to see how the Christoffel symbols are implemented.  *Hint 2*: You will need to use SymPy's `simplify()` function to simplify the expression obtained.)

# In[2]:


print("Check that gammaDD[0][1] = gammaDD[1][0]:")
print("gammaDD[0][1] = ", gammaDD[0][1])
print("gammaDD[1][0] = ", gammaDD[1][0])

print("\nOutput gammaUU[1][0] = ", gammaUU[1][0])
print("\nCheck that gammaUU[1][0] - gammaUU[0][1] = 0:")
print("gammaUU[1][0] - gammaUU[0][1] = ", gammaUU[1][0] - gammaUU[0][1])


# Define Christoffel symbols in terms of the inverse metric and metric first derivatives:
# $$
# \Gamma_{ij}^k = \frac{1}{2} \gamma^{kl}\left(\gamma_{jl,i} + \gamma_{il,j} - \gamma_{ij,l}\right)
# $$

# In[3]:


# First define symbolic expressions for metric derivatives
gammaDD_dD = ixp.declarerank3("gammaDD_dD", symmetry="sym01", DIM=3)

# Initialize GammaUDD (3-Christoffel) to zero
GammaUDD = ixp.zerorank3(DIM=3)
for i in range(3):
    for j in range(3):
        for k in range(3):
            for l in range(3):
                GammaUDD[k][i][j] += sp.Rational(1, 2) * gammaUU[k][l] * (
                    gammaDD_dD[j][l][i] + gammaDD_dD[i][l][j] - gammaDD_dD[i][j][l])


# Now let's confirm that $\Gamma^{k}_{ij} = \Gamma^k_{ji}$:

# In[4]:


for i in range(3):
    for j in range(3):
        for k in range(3):
            print(GammaUDD[i][j][k] - GammaUDD[i][k][j])


# <a id='outputc'></a>
# 
# # Part 3:  `outputC()` example: Output $\Gamma^i_{jk}$ expressions as optimized C code, assuming derivatives already specified \[Back to [top](#toc)\]
# $$\label{outputc}$$
# 
# At the core of NRPy+ is the ability to convert SymPy expressions to highly optimized C code.
# 
# **Problem statement**: Output all 18 unique Christoffel symbols with 3 different levels of optimization supported by NRPy+'s core C output routine `outputC()`.
# 
# First we store all 18 unique Christoffel symbols, as well as their desired variable names in C, to Python lists.

# In[5]:


symbols_list = []
varname_list = []
for i in range(3):
    for j in range(3):
        for k in range(j, 3):
            symbols_list += [GammaUDD[i][j][k]]
            varname_list += ["ChristoffelUDD" + str(i) + str(j) + str(k)]


# Next we input these lists into NRPy+'s C/C++ code generation function `outputC()`, at three different levels of optimization.
# 
# **Optimization Level 0 (don't ever do it this way)**: Compute each Christoffel symbol independently.

# In[6]:


import outputC as outC  # NRPy+: Core C code output module
outC.outputC(symbols_list, varname_list, filename="stdout", params="CSE_enable=False,outCverbose=False")


# Notice in the above code that many expressions are recomputed time and time again. This is *incredibly inefficient*, and generally compilers won't optimize this properly. So in optimization level 1, we use common subexpression elimination.
# 
# **Optimization Level 1**: Use [common subexpression elimination (CSE)](https://en.wikipedia.org/wiki/Common_subexpression_elimination) to group common subexpressions

# In[7]:


outC.outputC(symbols_list, varname_list, filename="stdout", params="CSE_enable=True,outCverbose=False")


# **Optimization Level 2**: Use CSE and take advantage of [single instruction, multiple data (SIMD) macros](https://en.wikipedia.org/wiki/Single_instruction,_multiple_data). NRPy+ translates these macros into assembler-level-optimized code for (x86_64) CPUs. Use case: looping over data on a numerical grid; can evaluate expressions at up to 8 gridpoints simultaneously *per CPU core*. Note that `FusedMulSubSIMD(a,b,c)` and `FusedMulAddSIMD(a,b,c)` perform *two* floating-point operations per cycle (on modern CPUs), whereas for other arithmetic operations at most one FP operation can be performed per cycle.

# In[8]:


outC.outputC(symbols_list, varname_list, filename="stdout", params="CSE_enable=True,enable_SIMD=True,outCverbose=False")


# <a id='fd_outputc'></a>
# 
# # Part 4:  `FD_outputC()` example: Specify numerical derivatives within $\Gamma^i_{jk}$ expressions as finite differences \[Back to [top](#toc)\]
# $$\label{fd_outputc}$$
# 
# To emphasize the infrastructure-agnostic nature of NRPy+, in the above C codes we assumed that the derivatives were already computed (e.g., using finite-difference derivatives, discontinuous Galerkin methods, pseudospectral methods, finite-element methods, etc.)
# 
# In this part, we demonstrate NRPy+'s `FD_outputC()` code, which basically prepends the above C/C++ codes with the code needed for computing arbitrary-order finite-difference derivatives.
# 
# To do this, NRPy+ makes the standard assumption that the underlying grid is *uniform* and that derivatives are taken with respect to functions stored at each point on the numerical grid. Appropriately, we call such functions "gridfunctions".
# 
# In the following, we assume that the input, $\gamma_{ij}$, is stored at each gridpoint. Also, we wish to store each independent component of $\Gamma^i_{jk}$ at each gridpoint as output. The below code cell also introduces the NRPy+ parameter interface (accessed via `NRPy_param_funcs.py`) to set the finite-differencing order to 6.
# 
# We'll use Optimization Level 1, to make the code easier to read; change `enable_SIMD=False` in the below `FD_outputC()` function call to `enable_SIMD=True` to see the Optimization Level 2 version.

# In[9]:


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 NRPy_param_funcs as par   # NRPy+: parameter interface

# First specify the finite-differencing order to be 6th order (all even orders > 0 supported!)
par.set_parval_from_str("finite_difference::FD_CENTDERIVS_ORDER", 6)

# Next register gridfunctions for Christoffel symbols
gri.register_gridfunctions("AUX", varname_list, rank=3, is_indexed=True, DIM=3)

# Then register gamma_{ij} as gridfunctions
ixp.register_gridfunctions_for_single_rank2("EVOL", "gammaDD", "sym01", DIM=3)

# Finally output the C code using FD_outputC()
# Step 1: FD_outputC() requires that left-hand side and right-hand side of each
#     expression be specified in a named-tuple "lhrh" defined in outputC.py
outputC_lhrh = []
for i, varname in enumerate(varname_list):
    outputC_lhrh += [outC.lhrh(lhs=gri.gfaccess("out_gf", varname), rhs=symbols_list[i])]
# Step 2: call FD_outputC.
fin.FD_outputC("stdout", outputC_lhrh, params="CSE_enable=True,enable_SIMD=False,outCverbose=False")


# <a id='what_next'></a>
# 
# # Part 5: What next? Navigating the NRPy+ tutorial \[Back to [top](#toc)\]
# $$\label{what_next}$$
# 
# As mentioned previously, NRPy+ is meant to be "learn by example". To that end, there are more than 100 fully documented Jupyter notebooks covering a wide variety of topics relevant to numerical relativity research and optimized algorithms for solving PDEs numerically.
# 
# So the answer to the question "*What next?*" is, naturally, "*Whatever you like!*" To continue the journey, check out [**the main NRPy+ tutorial table of contents**.](NRPyPlus_Tutorial.ipynb).

# <a id='latex_pdf_output'></a>
# 
# # Part 6: 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-NRPyPlus_10_Minute_Overview.pdf](Tutorial-NRPyPlus_10_Minute_Overview.pdf). (Note that clicking on this link may not work; you may need to open the PDF file through another means.)

# In[10]:


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