{ "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 }