{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"# Finite-Difference Playground: Using NRPy+-Generated C Codes in a Larger Project\n",
"\n",
"## Author: Zach Etienne\n",
"### Formatting improvements courtesy Brandon Clark\n",
"\n",
"## This notebook demonstrates the use of NRPy+ generated C codes in finite-difference derivatives computation and verifies fourth-order convergence. It invites further exploration of higher-order accuracies and highlights the advantages of NRPy+'s handling of complex higher-dimensional expressions.\n",
"\n",
"## Introduction: \n",
"To illustrate how NRPy+-based codes can be used, we write a C code that makes use of the NRPy+-generated C code from the [previous module](Tutorial-Finite_Difference_Derivatives.ipynb). This is a rather silly example, as the C code generated by NRPy+ could be easily generated by hand. However, as we will see in later modules, NRPy+'s true strengths lie in its automatic handling of far more complex and generic expressions, in higher dimensions. For the time being, bear with NRPy+; its true powers will become clear soon!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"This notebook is organized as follows\n",
"\n",
"1. [Step 1](#outputc): Register the C function `finite_diff_tutorial__second_deriv()`\n",
"1. [Step 2](#fdplayground): Finite-Difference Playground: A Complete C Code for Analyzing Finite-Difference Expressions Output by NRPy+\n",
"1. [Step 3](#exercise): Exercises\n",
"1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 1: Register the C function `finite_diff_tutorial__second_deriv()` \\[Back to [top](#toc)\\]\n",
"$$\\label{outputc}$$\n",
"\n",
"We start with the NRPy+ code from the [previous module](Tutorial-Finite_Difference_Derivatives.ipynb):"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:32:17.521891Z",
"iopub.status.busy": "2021-03-07T17:32:17.520599Z",
"iopub.status.idle": "2021-03-07T17:32:17.976712Z",
"shell.execute_reply": "2021-03-07T17:32:17.975986Z"
}
},
"outputs": [],
"source": [
"# Step P1: Import needed NRPy+ core modules:\n",
"import outputC as outC # NRPy+: Core C code output module\n",
"import NRPy_param_funcs as par # NRPy+: Parameter interface\n",
"import finite_difference as fin # NRPy+: Finite difference C code generation module\n",
"import grid as gri # NRPy+: Functions having to do with numerical grids\n",
"import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n",
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"import os, shutil # Standard Python: multiplatform OS funcs\n",
"\n",
"# Set the spatial dimension to 1\n",
"par.set_paramsvals_value(\"grid::DIM = 1\")\n",
"\n",
"# Register the input gridfunction \"phi\" and the gridfunction to which data are output, \"output\":\n",
"phi, output = gri.register_gridfunctions(\"AUX\",[\"phi\",\"output\"])\n",
"\n",
"# Declare phi_dDD as a rank-2 indexed expression: phi_dDD[i][j] = \\partial_i \\partial_j phi\n",
"phi_dDD = ixp.declarerank2(\"phi_dDD\",\"nosym\")\n",
"\n",
"# Set output to \\partial_0^2 phi\n",
"output = phi_dDD[0][0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now for the C code generation. First create the output directory for C codes."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Step P2: Create C code output directory:\n",
"Ccodesrootdir = os.path.join(\"FiniteDifferencePlayground_Ccodes/\")\n",
"# P2.a: First remove C code output directory if it exists\n",
"# Courtesy https://stackoverflow.com/questions/303200/how-do-i-remove-delete-a-folder-that-is-not-empty\n",
"shutil.rmtree(Ccodesrootdir, ignore_errors=True)\n",
"# P2.b: Then create a fresh directory\n",
"# Then create a fresh directory\n",
"cmd.mkdir(Ccodesrootdir)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then declare C macros that will appear within `NRPy_basic_defines.h`, `#include`d in most C files within NRPy+ generated C codes."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"outC.outC_NRPy_basic_defines_h_dict[\"FiniteDifferencePlayground\"] = r\"\"\"\n",
"// Declare the IDX2S(gf,i) macro, which enables us to store 2-dimensions of\n",
"// data in a 1D array. In this case, consecutive values of \"i\"\n",
"// (\"gf\" held to a fixed value) are consecutive in memory, where\n",
"// consecutive values of \"gf\" (fixing \"i\") are separated by N elements in\n",
"// memory.\n",
"#define IDX2S(gf, i) ( (i) + Npts_in_stencil * (gf) )\n",
"\n",
"// Declare PHIGF and OUTPUTGF, used for IDX2S's gf input.\n",
"#define PHIGF 0\n",
"#define OUTPUTGF 1\n",
"\"\"\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, register the C function `finite_diff_tutorial__second_deriv()`, which contains the core C-code kernel for evaluating the finite-difference derivative.\n",
"\n",
"Registering C functions with NRPy+ in this way encourages developers to abide by NRPy+ standards, in which \n",
"\n",
"* core functions are located within files of the same name\n",
"* code comments describing the function are well-formatted and appear at the top of the function\n",
"* function prototypes are automatically generated from information input into `add_to_Cfunction_dict()`\n",
"* generating a fully functional `Makefile` from lists of C functions is automated and provided by `outputC`'s `construct_Makefile_from_outC_function_dict()`.\n",
"\n",
"Further, it is conventional in NRPy+ to register C functions within a Python function, called at the bottom of a notebook/workflow. Turns out C-code generation of very complex expressions can take some time, and doing it via function calls enables us to parallelize the process (via Python's [multiprocessing](https://docs.python.org/3/library/multiprocessing.html) module), resulting in much faster turnaround times.\n",
"\n",
"Doing it this way also enables functions to be automatically generated with multiple options, with each function specializing in a given option. For example, one could imagine outputting multiple such functions for different finite difference orders. When compared to C++ template metaprogramming (\"[Difficulty level: Hard](https://www.geeksforgeeks.org/template-metaprogramming-in-c/)\"), it is far simpler to gain familiarity with this more explicit approach."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def add_to_Cfunction_dict_finite_diff_tutorial__second_deriv():\n",
" # includes field: list of files to #include at the top. NRPy_basic_defines.h defines C macros used here\n",
" # for accessing gridfunction data.\n",
" includes = [\"stdio.h\", \"NRPy_basic_defines.h\"] # Add stdio.h in case we'd like to use printf() for debugging.\n",
" prefunc = \"\" # would contain other functions that are *only* called by this function, declared static.\n",
" desc = r\"\"\"Evaluate phi''(x0) the second derivative of phi, a function sampled on a uniform grid.\n",
"\n",
"The derivative is evaluated using a standard 4th-order-accurate finite-difference approach.\n",
"\n",
"See Tutorial-Finite_Difference_Derivatives.ipynb in https://github.com/zachetienne/nrpytutorial\n",
"for more details.\n",
"\"\"\"\n",
" c_type = \"void\"\n",
" name = \"finite_diff_tutorial__second_deriv\"\n",
" params = \"const double *in_gfs, const double invdx0, const int i0, const int Npts_in_stencil, double *aux_gfs\"\n",
" body = fin.FD_outputC(\"returnstring\",\n",
" outC.lhrh(lhs=gri.gfaccess(\"aux_gfs\", \"output\"), rhs=output),\n",
" params=\"preindent=1,includebraces=False\") + \"\\n\"\n",
" outC.add_to_Cfunction_dict(\n",
" includes=includes,\n",
" prefunc=prefunc,\n",
" desc=desc,\n",
" c_type=c_type, name=name, params=params,\n",
" body=body,\n",
" enableCparameters=False)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take a look at NRPy+'s generated version of this function. All we need to do is call the above function `add_to_Cfunction_dict_finite_diff_tutorial__second_deriv()` to register the function, and then it can be accessed directly from `outputC`'s `outC_function_master_list` list."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:32:17.521891Z",
"iopub.status.busy": "2021-03-07T17:32:17.520599Z",
"iopub.status.idle": "2021-03-07T17:32:17.976712Z",
"shell.execute_reply": "2021-03-07T17:32:17.975986Z"
},
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"#include \"stdio.h\"\n",
"#include \"./NRPy_basic_defines.h\"\n",
"/*\n",
" * Evaluate phi''(x0) the second derivative of phi, a function sampled on a uniform grid.\n",
" * \n",
" * The derivative is evaluated using a standard 4th-order-accurate finite-difference approach.\n",
" * \n",
" * See Tutorial-Finite_Difference_Derivatives.ipynb in https://github.com/zachetienne/nrpytutorial\n",
" * for more details.\n",
" */\n",
"void finite_diff_tutorial__second_deriv(const double *in_gfs, const double invdx0, const int i0, const int Npts_in_stencil, double *aux_gfs) {\n",
"\n",
" /*\n",
" * NRPy+ Finite Difference Code Generation, Step 1 of 2: Read from main memory and compute finite difference stencils:\n",
" */\n",
" /*\n",
" * Original SymPy expression:\n",
" * \"const double phi_dDD00 = invdx0**2*(-5*phi/2 + 4*phi_i0m1/3 - phi_i0m2/12 + 4*phi_i0p1/3 - phi_i0p2/12)\"\n",
" */\n",
" const double phi_i0m2 = aux_gfs[IDX2S(PHIGF, i0-2)];\n",
" const double phi_i0m1 = aux_gfs[IDX2S(PHIGF, i0-1)];\n",
" const double phi = aux_gfs[IDX2S(PHIGF, i0)];\n",
" const double phi_i0p1 = aux_gfs[IDX2S(PHIGF, i0+1)];\n",
" const double phi_i0p2 = aux_gfs[IDX2S(PHIGF, i0+2)];\n",
" const double FDPart1_Rational_5_2 = 5.0/2.0;\n",
" const double FDPart1_Rational_1_12 = 1.0/12.0;\n",
" const double FDPart1_Rational_4_3 = 4.0/3.0;\n",
" const double phi_dDD00 = ((invdx0)*(invdx0))*(FDPart1_Rational_1_12*(-phi_i0m2 - phi_i0p2) + FDPart1_Rational_4_3*(phi_i0m1 + phi_i0p1) - FDPart1_Rational_5_2*phi);\n",
" /*\n",
" * NRPy+ Finite Difference Code Generation, Step 2 of 2: Evaluate SymPy expressions and write to main memory:\n",
" */\n",
" /*\n",
" * Original SymPy expression:\n",
" * \"aux_gfs[IDX2S(OUTPUTGF, i0)] = phi_dDD00\"\n",
" */\n",
" aux_gfs[IDX2S(OUTPUTGF, i0)] = phi_dDD00;\n",
"\n",
"}\n",
"\n"
]
}
],
"source": [
"add_to_Cfunction_dict_finite_diff_tutorial__second_deriv()\n",
"for item in outC.outC_function_master_list:\n",
" print(outC.outC_function_dict[item.name])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 2: Finite-Difference Playground: A Complete C Code for Analyzing Finite-Difference Expressions Output by NRPy+ \\[Back to [top](#toc)\\]\n",
"$$\\label{fdplayground}$$\n",
"\n",
"NRPy+ is designed to generate C code \"kernels\" at the heart of more advanced projects. As an example of its utility, let's now write a simple C code that imports the above file `finite_diff_tutorial-second_deriv.h` to evaluate the finite-difference second derivative of\n",
"\n",
"$$f(x) = \\sin(x)$$\n",
"\n",
"at fourth-order accuracy. Let's call the finite-difference second derivative of $f$ evaluated at a point $x$ $f''(x)_{\\rm FD}$. A fourth-order-accurate $f''(x)_{\\rm FD}$ will, in the truncation-error-dominated regime, satisfy the equation\n",
"\n",
"$$f''(x)_{\\rm FD} = f''(x)_{\\rm exact} + \\mathcal{O}(\\Delta x^4).$$\n",
"\n",
"Therefore, the [relative error](https://en.wikipedia.org/wiki/Approximation_error) between the finite-difference derivative and the exact value should be given to good approximation by\n",
"\n",
"$$E_{\\rm Rel} = \\left| \\frac{f''(x)_{\\rm FD} - f''(x)_{\\rm exact}}{f''(x)_{\\rm exact}}\\right| \\propto \\Delta x^4,$$\n",
"\n",
"so that (taking the logarithm of both sides of the equation):\n",
"\n",
"$$\\log_{10} E_{\\rm Rel} = 4 \\log_{10} (\\Delta x) + \\log_{10} (k),$$\n",
"\n",
"where $k$ is the proportionality constant, divided by $f''(x)_{\\rm exact}$.\n",
"\n",
"Let's confirm this is true using our finite-difference playground code, which imports the NRPy+-generated C code generated above for evaluating $f''(x)_{\\rm FD}$ at fourth-order accuracy, and outputs $\\log_{10} (\\Delta x)$ and $\\log_{10} E_{\\rm Rel}$ in a range of $\\Delta x$ that is truncation-error dominated."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"def add_to_Cfunction_dict_main__finite_difference_playground():\n",
" includes = [\"stdio.h\", \"stdlib.h\", \"math.h\", \"NRPy_basic_defines.h\", \"NRPy_function_prototypes.h\"]\n",
" prefunc = r\"\"\"\n",
"// Define the function we wish to differentiate, as well as its exact second derivative:\n",
"double f(const double x) { return sin(x); } // f(x)\n",
"double f_dDD_exact(const double x) { return -sin(x); } // f''(x)\n",
"\n",
"// Define the uniform grid coordinate:\n",
"// x_i = (x_0 + i*Delta_x)\n",
"double x_i(const double x_0,const int i,const double Delta_x) {\n",
" return (x_0 + (double)i*Delta_x);\n",
"}\n",
"\"\"\"\n",
" desc = \"\"\"// main() function:\n",
"// Step 0: Read command-line input, set up grid structure, allocate memory for gridfunctions, set up coordinates\n",
"// Step 1: Write test data to gridfunctions\n",
"// Step 2: Overwrite all data in ghost zones with NaNs\n",
"// Step 3: Apply curvilinear boundary conditions\n",
"// Step 4: Print gridfunction data after curvilinear boundary conditions have been applied\n",
"// Step 5: Free all allocated memory\n",
"\"\"\"\n",
" c_type = \"int\"\n",
" name = \"main\"\n",
" params = \"int argc, const char *argv[]\"\n",
" body = r\"\"\" // Step 0: Read command-line arguments (TODO)\n",
"\n",
" // Step 1: Set some needed constants\n",
" const int Npts_in_stencil = 5; // Equal to the finite difference order, plus one. '+str(par.parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\"))+'+1;\n",
" const double PI = 3.14159265358979323846264338327950288; // The scale over which the sine function varies.\n",
" const double x_eval = PI/4.0; // x_0 = desired x at which we wish to compute f(x)\n",
"\n",
" // Step 2: Evaluate f''(x_eval) using the exact expression:\n",
" double EX = f_dDD_exact(x_eval);\n",
"\n",
" // Step 3: Allocate space for gridfunction input and output.\n",
" double *restrict gfs = (double *restrict)malloc(sizeof(double)*Npts_in_stencil*2);\n",
"\n",
" // Step 4: Loop over grid spacings\n",
" for(double Delta_x = 1e-3*(2*PI);Delta_x<=1.5e-1*(2*PI);Delta_x*=1.1) {\n",
"\n",
" // Step 4a: x_eval is the center point of the finite differencing stencil,\n",
" // thus x_0 = x_eval - 2*dx for fourth-order-accurate first & second finite difference derivs,\n",
" // and x_0 = x_eval - 3*dx for sixth-order-accurate first & second finite difference derivs, etc.\n",
" // In general, for the integer Npts_in_stencil, we have\n",
" // x_0 = x_eval - (double)(Npts_in_stencil/2)*Delta_x,\n",
" // where we rely upon integer arithmetic (which always rounds down) to ensure\n",
" // Npts_in_stencil/2 = 5/2 = 2 for fourth-order-accurate first & second finite difference derivs:\n",
" const double x_0 = x_eval - (double)(Npts_in_stencil/2)*Delta_x;\n",
"\n",
" // Step 4b: Set \\phi=PHIGF to be f(x) as defined in the\n",
" // f(const double x) function above, where x_i = stencil_start_x + i*Delta_x:\n",
" for(int ii=0;ii"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# from https://stackoverflow.com/questions/52386747/matplotlib-check-for-empty-plot\n",
"import numpy\n",
"x, y = numpy.loadtxt(fname='data.txt', delimiter='\\t', unpack=True)\n",
"fig, ax = plt.subplots()\n",
"ax.set_xlabel('log10(Delta_x)')\n",
"ax.set_ylabel('log10([Relative Error])')\n",
"ax.plot(x, y, color='b')\n",
"\n",
"os.chdir(\"..\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A quick glance at the above plot indicates that between $\\log_{10}(\\Delta x) \\approx -2.0$ and $\\log_{10}(\\Delta x) \\approx -1.0$, the logarithmic relative error $\\log_{10} E_{\\rm Rel}$ increases by about 4, indicating a positive slope of approximately 4. Thus we have confirmed fourth-order convergence."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 3: Exercises \\[Back to [top](#toc)\\]\n",
"$$\\label{exercise}$$\n",
"\n",
"1. Use NumPy's [`polyfit()`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html) function to evaluate the least-squares slope of the above line.\n",
"2. Explore $\\log_{10}(\\Delta x)$ outside the above (truncation-error-dominated) range. What other errors dominate outside the truncation-error-dominated regime?\n",
"3. Adjust the above NRPy+ and C codes to support 6th-order-accurate finite differencing. What should the slope of the resulting plot of $\\log_{10} E_{\\rm Rel}$ versus $\\log_{10}(\\Delta x)$ be? Explain why this case does not provide as clean a slope as the 4th-order case."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 4: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
"$$\\label{latex_pdf_output}$$\n",
"\n",
"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\n",
"[Tutorial-Start_to_Finish-Finite_Difference_Playground.pdf](Tutorial-Start_to_Finish-Finite_Difference_Playground.pdf). (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-03-07T17:32:18.967244Z",
"iopub.status.busy": "2021-03-07T17:32:18.966638Z",
"iopub.status.idle": "2021-03-07T17:32:22.469505Z",
"shell.execute_reply": "2021-03-07T17:32:22.468459Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Created Tutorial-Start_to_Finish-Finite_Difference_Playground.tex, and\n",
" compiled LaTeX file to PDF file Tutorial-Start_to_Finish-\n",
" Finite_Difference_Playground.pdf\n"
]
}
],
"source": [
"import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n",
"cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Start_to_Finish-Finite_Difference_Playground\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.1"
}
},
"nbformat": 4,
"nbformat_minor": 4
}