{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Start-to-Finish Example: Unit Testing `GiRaFFE_NRPy`: Boundary Conditions\n", "\n", "## Author: Patrick Nelson\n", "\n", "#### Edits by Terrence Pierre Jacques\n", "\n", "## This module Validates the Boundary Conditions routines for `GiRaFFE_NRPy`.\n", "\n", "**Notebook Status:** Validated\n", "\n", "**Validation Notes:** This module will validate the routines in [Tutorial-GiRaFFE_NRPy-BCs](Tutorial-GiRaFFE_NRPy-BCs.ipynb).\n", "\n", "### NRPy+ Source Code for this module: \n", "* [GiRaFFE_NRPy/GiRaFFE_NRPy_BCs.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_BCs.py) [\\[**tutorial**\\]](Tutorial-GiRaFFE_NRPy-BCs.ipynb) Generates the driver apply boundary conditions to the vector potential and velocity.\n", "\n", "## Introduction:\n", "\n", "This notebook validates the C code used to apply boundary conditions to components of the vector potential and valencia velocity.\n", "\n", "It is, in general, good coding practice to unit test functions individually to verify that they produce the expected and intended output. We will generate test data with arbitrarily-chosen analytic functions and calculate gridfunctions at the cell centers on a small numeric grid. We will then compute the values for the ghost zones in two ways: first with the boundary condition C code driver, then we compute them analytically.\n", "\n", "When this notebook is run, the significant digits of agreement between the approximate and exact values in the ghost zones will be evaluated. If the agreement falls below a thresold, the point, quantity, and level of agreement are reported [here](#compile_run)." ] }, { "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](#setup): Set up core functions and parameters for unit testing the BCs algorithm\n", " 1. [Step 1.a](#expressions) Write expressions for the gridfunctions we will test\n", " 1. [Step 1.b](#ccodekernels) Generate C functions to calculate the gridfunctions\n", " 1. [Step 1.c](#free_parameters) Set free parameters in the code\n", "1. [Step 2](#mainc): `BCs_unit_test.c`: The Main C Code\n", " 1. [Step 2.a](#compile_run): Compile and run the code\n", "1. [Step 3](#convergence): Code validation: Verify that relative error in numerical solution converges to zero at the expected order\n", "1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 1: Set up core functions and parameters for unit testing the BCs algorithm \\[Back to [top](#toc)\\]\n", "$$\\label{setup}$$\n", "\n", "We'll start by appending the relevant paths to `sys.path` so that we can access sympy modules in other places. Then, we'll import NRPy+ core functionality and set up a directory in which to carry out our test. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2022-11-24T14:15:17.501444Z", "iopub.status.busy": "2022-11-24T14:15:17.501240Z", "iopub.status.idle": "2022-11-24T14:15:17.696580Z", "shell.execute_reply": "2022-11-24T14:15:17.696303Z" } }, "outputs": [], "source": [ "import os, sys, shutil # Standard Python modules for multiplatform OS-level functions\n", "# First, we'll add the parent directory to the list of directories Python will check for modules.\n", "nrpy_dir_path = os.path.join(\"..\")\n", "if nrpy_dir_path not in sys.path:\n", " sys.path.append(nrpy_dir_path)\n", "nrpy_dir_path = os.path.join(\"..\", \"Deprecated\")\n", "if nrpy_dir_path not in sys.path:\n", " sys.path.append(nrpy_dir_path)\n", "\n", "from outputC import outCfunction, lhrh # NRPy+: Core C code output module\n", "import sympy as sp # SymPy: The Python computer algebra package upon which NRPy+ depends\n", "import finite_difference as fin # NRPy+: Finite difference C code generation module\n", "import NRPy_param_funcs as par # NRPy+: Parameter interface\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 reference_metric as rfm # NRPy+: Reference metric support\n", "import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n", "\n", "Ccodesdir = \"Start-to-Finish-UnitTests/BCs_UnitTest/\"\n", "\n", "# 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(Ccodesdir, ignore_errors=True)\n", "# Then create a fresh directory\n", "cmd.mkdir(Ccodesdir)\n", "\n", "outdir = os.path.join(Ccodesdir,\"output/\")\n", "cmd.mkdir(outdir)\n", "\n", "thismodule = \"Start_to_Finish_UnitTest-GiRaFFE_NRPy-BCs\"\n", "\n", "# Set the finite-differencing order to 2\n", "par.set_parval_from_str(\"finite_difference::FD_CENTDERIVS_ORDER\", 2)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 1.a: Write expressions for the gridfunctions we will test \\[Back to [top](#toc)\\]\n", "$$\\label{expressions}$$\n", "\n", "Now, we'll choose some functions with arbitrary forms to generate test data. We'll need to set seven gridfunctions, so expressions are being pulled from several previously written unit tests.\n", "\n", "\\begin{align}\n", "A_x &= dy + ez + f \\\\\n", "A_y &= mx + nz + p \\\\\n", "A_z &= sx + ty + u. \\\\\n", "\\bar{v}^x &= ax + by + cz \\\\\n", "\\bar{v}^y &= bx + cy + az \\\\\n", "\\bar{v}^z &= cx + ay + bz \\\\\n", "[\\sqrt{\\gamma} \\Phi] &= 1 - (x+2y+z) \\\\\n", "\\end{align}\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2022-11-24T14:15:17.698741Z", "iopub.status.busy": "2022-11-24T14:15:17.698636Z", "iopub.status.idle": "2022-11-24T14:15:17.734875Z", "shell.execute_reply": "2022-11-24T14:15:17.734596Z" } }, "outputs": [], "source": [ "a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t,u = par.Cparameters(\"REAL\",thismodule,[\"a\",\"b\",\"c\",\"d\",\"e\",\"f\",\"g\",\"h\",\"l\",\"m\",\"n\",\"o\",\"p\",\"q\",\"r\",\"s\",\"t\",\"u\"],1e300)\n", "M_PI = par.Cparameters(\"#define\",thismodule,[\"M_PI\"], \"\")\n", "\n", "AD = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"AD\",DIM=3)\n", "ValenciavU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"ValenciavU\",DIM=3)\n", "psi6Phi = gri.register_gridfunctions(\"EVOL\",\"psi6Phi\")\n", "\n", "par.set_parval_from_str(\"reference_metric::CoordSystem\",\"Cartesian\")\n", "rfm.reference_metric()\n", "x = rfm.xx_to_Cart[0]\n", "y = rfm.xx_to_Cart[1]\n", "z = rfm.xx_to_Cart[2]\n", "\n", "AD[0] = d*y + e*z + f\n", "AD[1] = m*x + n*z + o\n", "AD[2] = s*x + t*y + u\n", "\n", "ValenciavU[0] = a#*x + b*y + c*z\n", "ValenciavU[1] = b#*x + c*y + a*z\n", "ValenciavU[2] = c#*x + a*y + b*z\n", "\n", "psi6Phi = sp.sympify(1) - (x + sp.sympify(2)*y + z)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 1.b: Generate C functions to calculate the gridfunctions \\[Back to [top](#toc)\\]\n", "$$\\label{ccodekernels}$$\n", "\n", "Here, we will use the NRPy+ function `outCfunction()` to generate C code that will calculate our metric gridfunctions over an entire grid; note that we call the function twice, once over just the interior points, and once over all points. This will allow us to compare against exact values in the ghostzones. We will also call the function to generate the boundary conditions function we are testing. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2022-11-24T14:15:17.737178Z", "iopub.status.busy": "2022-11-24T14:15:17.737078Z", "iopub.status.idle": "2022-11-24T14:15:17.782989Z", "shell.execute_reply": "2022-11-24T14:15:17.782709Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output C function calculate_test_data() to file Start-to-Finish-UnitTests/BCs_UnitTest/calculate_test_data.h\n", "Output C function calculate_test_data_exact() to file Start-to-Finish-UnitTests/BCs_UnitTest/calculate_test_data_exact.h\n" ] } ], "source": [ "metric_gfs_to_print = [\\\n", " lhrh(lhs=gri.gfaccess(\"evol_gfs\",\"AD0\"),rhs=AD[0]),\\\n", " lhrh(lhs=gri.gfaccess(\"evol_gfs\",\"AD1\"),rhs=AD[1]),\\\n", " lhrh(lhs=gri.gfaccess(\"evol_gfs\",\"AD2\"),rhs=AD[2]),\\\n", " lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"ValenciavU0\"),rhs=ValenciavU[0]),\\\n", " lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"ValenciavU1\"),rhs=ValenciavU[1]),\\\n", " lhrh(lhs=gri.gfaccess(\"auxevol_gfs\",\"ValenciavU2\"),rhs=ValenciavU[2]),\\\n", " lhrh(lhs=gri.gfaccess(\"evol_gfs\",\"psi6Phi\"),rhs=psi6Phi),\\\n", " ]\n", "\n", "desc = \"Calculate test data on the interior grid for boundary conditions\"\n", "name = \"calculate_test_data\"\n", "outCfunction(\n", " outfile = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n", " params =\"const paramstruct *restrict params,REAL *restrict xx[3],REAL *restrict auxevol_gfs,REAL *restrict evol_gfs\",\n", " body = fin.FD_outputC(\"returnstring\",metric_gfs_to_print,params=\"outCverbose=False\"),\n", " loopopts=\"InteriorPoints,Read_xxs\")\n", "\n", "desc = \"Calculate test data at all points for comparison\"\n", "name = \"calculate_test_data_exact\"\n", "outCfunction(\n", " outfile = os.path.join(Ccodesdir,name+\".h\"), desc=desc, name=name,\n", " params =\"const paramstruct *restrict params,REAL *restrict xx[3],REAL *restrict auxevol_gfs,REAL *restrict evol_gfs\",\n", " body = fin.FD_outputC(\"returnstring\",metric_gfs_to_print,params=\"outCverbose=False\"),\n", " loopopts=\"AllPoints,Read_xxs\")\n", "\n", "import GiRaFFE_NRPy.GiRaFFE_NRPy_BCs as BC\n", "BC.GiRaFFE_NRPy_BCs(os.path.join(Ccodesdir,\"boundary_conditions\"))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 1.c: Set free parameters in the code \\[Back to [top](#toc)\\]\n", "$$\\label{free_parameters}$$\n", "\n", "We also need to create the files that interact with NRPy's C parameter interface. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2022-11-24T14:15:17.799919Z", "iopub.status.busy": "2022-11-24T14:15:17.799834Z", "iopub.status.idle": "2022-11-24T14:15:17.802985Z", "shell.execute_reply": "2022-11-24T14:15:17.802776Z" } }, "outputs": [], "source": [ "# Step 3.d.i: Generate declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n", "# par.generate_Cparameters_Ccodes(os.path.join(out_dir))\n", "\n", "# Step 3.d.ii: Set free_parameters.h\n", "with open(os.path.join(Ccodesdir,\"free_parameters.h\"),\"w\") as file:\n", " file.write(\"\"\"\n", "// Override parameter defaults with values based on command line arguments and NGHOSTS.\n", "params.Nxx0 = atoi(argv[1]);\n", "params.Nxx1 = atoi(argv[2]);\n", "params.Nxx2 = atoi(argv[3]);\n", "params.Nxx_plus_2NGHOSTS0 = params.Nxx0 + 2*NGHOSTS;\n", "params.Nxx_plus_2NGHOSTS1 = params.Nxx1 + 2*NGHOSTS;\n", "params.Nxx_plus_2NGHOSTS2 = params.Nxx2 + 2*NGHOSTS;\n", "// Step 0d: Set up space and time coordinates\n", "// Step 0d.i: Declare \\Delta x^i=dxx{0,1,2} and invdxx{0,1,2}, as well as xxmin[3] and xxmax[3]:\n", "const REAL xxmin[3] = {-1.0,-1.0,-1.0};\n", "const REAL xxmax[3] = { 1.0, 1.0, 1.0};\n", "\n", "params.dxx0 = (xxmax[0] - xxmin[0]) / ((REAL)params.Nxx_plus_2NGHOSTS0-1.0);\n", "params.dxx1 = (xxmax[1] - xxmin[1]) / ((REAL)params.Nxx_plus_2NGHOSTS1-1.0);\n", "params.dxx2 = (xxmax[2] - xxmin[2]) / ((REAL)params.Nxx_plus_2NGHOSTS2-1.0);\n", "params.invdx0 = 1.0 / params.dxx0;\n", "params.invdx1 = 1.0 / params.dxx1;\n", "params.invdx2 = 1.0 / params.dxx2;\n", "\\n\"\"\")\n", "\n", "# Generates declare_Cparameters_struct.h, set_Cparameters_default.h, and set_Cparameters[-SIMD].h\n", "import deprecated_NRPy_param_funcs as evil_par\n", "evil_par.generate_Cparameters_Ccodes(os.path.join(Ccodesdir))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: `BCs_unit_test.c`: The Main C Code \\[Back to [top](#toc)\\]\n", "$$\\label{mainc}$$\n", "\n", "Here we compare the results of our boundary conditions, `apply_bcs_potential()` and `apply_bcs_velocity()`, against the exact results, looping over the entire numerical grid." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2022-11-24T14:15:17.804213Z", "iopub.status.busy": "2022-11-24T14:15:17.804140Z", "iopub.status.idle": "2022-11-24T14:15:17.808115Z", "shell.execute_reply": "2022-11-24T14:15:17.807894Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing Start-to-Finish-UnitTests/BCs_UnitTest//BCs_unit_test.c\n" ] } ], "source": [ "%%writefile $Ccodesdir/BCs_unit_test.c\n", "\n", "// These are common packages that we are likely to need.\n", "#include \"stdio.h\"\n", "#include \"stdlib.h\"\n", "#include \"math.h\"\n", "#include \"string.h\" // Needed for strncmp, etc.\n", "#include \"stdint.h\" // Needed for Windows GCC 6.x compatibility\n", "#include // Needed to set a random seed.\n", "\n", "#define REAL double\n", "#include \"declare_Cparameters_struct.h\"\n", "\n", "const int NGHOSTS = 3;\n", "\n", "REAL a,b,c,d,e,f,g,h,l,m,n,o,p,q,r,s,t,u;\n", "\n", "// Standard NRPy+ memory access:\n", "#define IDX4S(g,i,j,k) \\\n", "( (i) + Nxx_plus_2NGHOSTS0 * ( (j) + Nxx_plus_2NGHOSTS1 * ( (k) + Nxx_plus_2NGHOSTS2 * (g) ) ) )\n", "\n", "// Standard formula to calculate significant digits of agreement:\n", "#define SDA(a,b) 1.0-log10(2.0*fabs(a-b)/(fabs(a)+fabs(b)))\n", "\n", "// Give gridfunctions their names:\n", "#define VALENCIAVU0GF 0\n", "#define VALENCIAVU1GF 1\n", "#define VALENCIAVU2GF 2\n", "#define NUM_AUXEVOL_GFS 3\n", "\n", "#define AD0GF 0\n", "#define AD1GF 1\n", "#define AD2GF 2\n", "#define STILDED0GF 3\n", "#define STILDED1GF 4\n", "#define STILDED2GF 5\n", "#define PSI6PHIGF 6\n", "#define NUM_EVOL_GFS 7\n", "\n", "#include \"calculate_test_data.h\"\n", "#include \"calculate_test_data_exact.h\"\n", "#include \"boundary_conditions/GiRaFFE_boundary_conditions.h\"\n", "\n", "int main(int argc, const char *argv[]) {\n", " paramstruct params;\n", "#include \"set_Cparameters_default.h\"\n", "\n", " // Step 0c: Set free parameters, overwriting Cparameters defaults\n", " // by hand or with command-line input, as desired.\n", "#include \"free_parameters.h\"\n", "#include \"set_Cparameters-nopointer.h\"\n", "\n", " // We'll define our grid slightly different from how we normally would. We let our outermost\n", " // ghostzones coincide with xxmin and xxmax instead of the interior of the grid. This means\n", " // that the ghostzone points will have identical positions so we can do convergence tests of them.\n", " // Step 0e: Set up cell-centered Cartesian coordinate grids\n", " REAL *xx[3];\n", " xx[0] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS0);\n", " xx[1] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS1);\n", " xx[2] = (REAL *)malloc(sizeof(REAL)*Nxx_plus_2NGHOSTS2);\n", " for(int j=0;j\n", "\n", "## Step 2.a: Compile and run the code \\[Back to [top](#toc)\\]\n", "$$\\label{compile_run}$$\n", "\n", "Now that we have our file, we can compile it and run the executable." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2022-11-24T14:15:17.810084Z", "iopub.status.busy": "2022-11-24T14:15:17.809978Z", "iopub.status.idle": "2022-11-24T14:15:18.441418Z", "shell.execute_reply": "2022-11-24T14:15:18.440988Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Now compiling, should take ~2 seconds...\n", "\n", "Compiling executable...\n", "(EXEC): Executing `gcc -std=gnu99 -Ofast -fopenmp -march=native -funroll-loops Start-to-Finish-UnitTests/BCs_UnitTest/BCs_unit_test.c -o Start-to-Finish-UnitTests/BCs_UnitTest/output/BCs_unit_test -lm`...\n", "(BENCH): Finished executing in 0.60 seconds.\n", "Finished compilation.\n", "Finished in 0.6055264472961426 seconds.\n", "\n", "\n", "Now running...\n", "\n", "(EXEC): Executing `taskset -c 1,3,5,7,9,11,13,15 ./BCs_unit_test 2 2 2`...\n", "(BENCH): Finished executing in 0.20 seconds.\n", "Finished in 0.23199677467346191 seconds.\n", "\n", "\n" ] } ], "source": [ "import time\n", "results_file = \"out.txt\"\n", "\n", "\n", "print(\"Now compiling, should take ~2 seconds...\\n\")\n", "start = time.time()\n", "cmd.C_compile(os.path.join(Ccodesdir,\"BCs_unit_test.c\"), os.path.join(outdir,\"BCs_unit_test\"))\n", "end = time.time()\n", "print(\"Finished in \"+str(end-start)+\" seconds.\\n\\n\")\n", "\n", "# Change to output directory\n", "os.chdir(outdir)\n", "\n", "print(\"Now running...\\n\")\n", "start = time.time()\n", "cmd.Execute(os.path.join(\"BCs_unit_test\"),\"2 2 2\",file_to_redirect_stdout=results_file)\n", "# To do a convergence test, we'll also need a second grid with twice the resolution.\n", "# cmd.Execute(os.path.join(\"Validation\",\"BCs_unit_test\"),\"9 9 9\",file_to_redirect_stdout=os.path.join(out_dir,results_file))\n", "end = time.time()\n", "print(\"Finished in \"+str(end-start)+\" seconds.\\n\\n\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we will interpret our output and verify that we produced the correct results." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "execution": { "iopub.execute_input": "2022-11-24T14:15:18.443448Z", "iopub.status.busy": "2022-11-24T14:15:18.443323Z", "iopub.status.idle": "2022-11-24T14:15:18.446383Z", "shell.execute_reply": "2022-11-24T14:15:18.446101Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "All quantities agree at all points!\n", "\n" ] } ], "source": [ "with open(results_file,\"r\") as file:\n", " output = file.readline()\n", " print(output)\n", " if output!=\"All quantities agree at all points!\\n\": # If this isn't the first line of this file, something went wrong!\n", " sys.exit(1)" ] }, { "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-GiRaFFE_NRPy-Metric_Face_Values.pdf](Tutorial-Start_to_Finish-GiRaFFE_NRPy-Metric_Face_Values.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": 8, "metadata": { "execution": { "iopub.execute_input": "2022-11-24T14:15:18.447858Z", "iopub.status.busy": "2022-11-24T14:15:18.447745Z", "iopub.status.idle": "2022-11-24T14:15:21.464180Z", "shell.execute_reply": "2022-11-24T14:15:21.463926Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Created Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-BCs.tex, and\n", " compiled LaTeX file to PDF file Tutorial-Start_to_Finish_UnitTest-\n", " GiRaFFE_NRPy-BCs.pdf\n" ] } ], "source": [ "import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n", "\n", "# Change to NRPy directory\n", "os.chdir(\"../../../\")\n", "\n", "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Start_to_Finish_UnitTest-GiRaFFE_NRPy-BCs\")" ] } ], "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.10.11" } }, "nbformat": 4, "nbformat_minor": 4 }