{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# `GiRaFFE_NRPy`: Source Terms\n", "\n", "## Author: Patrick Nelson\n", "\n", "\n", "\n", "**Notebook Status:** Validated \n", "\n", "**Validation Notes:** This code produces the expected results for generated functions.\n", "\n", "## This module presents the functionality of [GiRaFFE_NRPy_Source_Terms.py](../../edit/in_progress/GiRaFFE_NRPy/GiRaFFE_NRPy_Source_Terms.py).\n", "\n", "## Introduction: \n", "This writes and documents the C code that `GiRaFFE_NRPy` uses to compute the source terms for the right-hand sides of the evolution equations for the unstaggered prescription.\n", "\n", "The equations themselves are already coded up in other functions; however, for the $\\tilde{S}_i$ source term, we will need derivatives of the metric. It will be most efficient and accurate to take them using the interpolated metric values that we will have calculated anyway; however, we will need to write our derivatives in a nonstandard way within NRPy+ in order to take advantage of this, writing our own code for memory access." ] }, { "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](#stilde_source): The $\\tilde{S}_i$ source term\n", "1. [Step 2](#code_validation): Code Validation against original C code\n", "1. [Step 3](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Step 0: Add NRPy's directory to the path\n", "# https://stackoverflow.com/questions/16780014/import-file-from-parent-directory\n", "import os,sys\n", "nrpy_dir_path = os.path.join(\"..\")\n", "if nrpy_dir_path not in sys.path:\n", " sys.path.append(nrpy_dir_path)\n", "\n", "import cmdline_helper as cmd\n", "outdir = os.path.join(\"GiRaFFE_NRPy\",\"GiRaFFE_Ccode_validation\",\"RHSs\")\n", "cmd.mkdir(outdir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 1: The $\\tilde{S}_i$ source term \\[Back to [top](#toc)\\]\n", "$$\\label{stilde_source}$$\n", "\n", "We start in the usual way - import the modules we need. We will also import the Levi-Civita symbol from `indexedexp.py` and use it to set the Levi-Civita tensor $\\epsilon^{ijk} = [ijk]/\\sqrt{\\gamma}$." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Step 1: The StildeD RHS *source* term\n", "from outputC import outputC, outCfunction, add_to_Cfunction_dict # NRPy+: Core C code output module\n", "import indexedexp as ixp # NRPy+: Symbolic indexed expression (e.g., tensors, vectors, etc.) support\n", "import GRHD.equations as GRHD # NRPy+: Generate general relativistic hydrodynamics equations\n", "import GRFFE.equations as GRFFE # NRPy+: Generate general relativistic force-free electrodynamics equations\n", "\n", "thismodule = \"GiRaFFE_NRPy_Source_Terms\"\n", "\n", "def generate_memory_access_code():\n", " # There are several pieces of C code that we will write ourselves because we need to do things\n", " # a little bit outside of what NRPy+ is built for.\n", " # First, we will write general memory access. We will read in values from memory at a given point\n", " # for each quantity we care about.\n", " global general_access\n", " general_access = \"\"\n", " for var in [\"GAMMADD00\", \"GAMMADD01\", \"GAMMADD02\",\n", " \"GAMMADD11\", \"GAMMADD12\", \"GAMMADD22\",\n", " \"BETAU0\", \"BETAU1\", \"BETAU2\",\"ALPHA\",\n", " \"BU0\",\"BU1\",\"BU2\",\n", " \"VALENCIAVU0\",\"VALENCIAVU1\",\"VALENCIAVU2\"]:\n", " lhsvar = var.lower().replace(\"dd\",\"DD\").replace(\"u\",\"U\").replace(\"bU\",\"BU\").replace(\"valencia\",\"Valencia\")\n", " # e.g.,\n", " # const REAL gammaDD00dD0 = auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1,i2)];\n", " general_access += \"const REAL \"+lhsvar+\" = auxevol_gfs[IDX4S(\"+var+\"GF,i0,i1,i2)];\\n\"\n", "\n", " # This quick function returns a nearby point for memory access. We need this because derivatives are not local operations.\n", " def idxp1(dirn):\n", " if dirn==0:\n", " return \"i0+1,i1,i2\"\n", " if dirn==1:\n", " return \"i0,i1+1,i2\"\n", " if dirn==2:\n", " return \"i0,i1,i2+1\"\n", "\n", " # Next we evaluate needed derivatives of the metric, based on their values at cell faces\n", " global metric_deriv_access\n", " metric_deriv_access = []\n", " for dirn in range(3):\n", " metric_deriv_access.append(\"\")\n", " for var in [\"GAMMA_FACEDDdD00\", \"GAMMA_FACEDDdD01\", \"GAMMA_FACEDDdD02\",\n", " \"GAMMA_FACEDDdD11\", \"GAMMA_FACEDDdD12\", \"GAMMA_FACEDDdD22\",\n", " \"BETA_FACEUdD0\", \"BETA_FACEUdD1\", \"BETA_FACEUdD2\",\"ALPHA_FACEdD\"]:\n", " lhsvar = var.lower().replace(\"dddd\",\"DDdD\").replace(\"udd\",\"UdD\").replace(\"dd\",\"dD\").replace(\"u\",\"U\").replace(\"_face\",\"\")\n", " rhsvar = var.replace(\"dD\",\"\")\n", " # e.g.,\n", " # const REAL gammaDDdD000 = (auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0+1,i1,i2)]-auxevol_gfs[IDX4S(GAMMA_FACEDD00GF,i0,i1,i2)])/dxx0;\n", " metric_deriv_access[dirn] += \"const REAL \"+lhsvar+str(dirn)+\" = (auxevol_gfs[IDX4S(\"+rhsvar+\"GF,\"+idxp1(dirn)+\")]-auxevol_gfs[IDX4S(\"+rhsvar+\"GF,i0,i1,i2)])/dxx\"+str(dirn)+\";\\n\"\n", " metric_deriv_access[dirn] += \"REAL Stilde_rhsD\"+str(dirn)+\";\\n\"\n", "\n", " # This creates the C code that writes to the Stilde_rhs direction specified.\n", " global write_final_quantity\n", " write_final_quantity = []\n", " for dirn in range(3):\n", " write_final_quantity.append(\"\")\n", " write_final_quantity[dirn] += \"rhs_gfs[IDX4S(STILDED\"+str(dirn)+\"GF,i0,i1,i2)] += Stilde_rhsD\"+str(dirn)+\";\"\n", "\n", "def write_out_functions_for_StildeD_source_term(outdir,outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi):\n", " generate_memory_access_code()\n", " # First, we declare some dummy tensors that we will use for the codegen.\n", " gammaDDdD = ixp.declarerank3(\"gammaDDdD\",\"sym01\",DIM=3)\n", " betaUdD = ixp.declarerank2(\"betaUdD\",\"nosym\",DIM=3)\n", " alphadD = ixp.declarerank1(\"alphadD\",DIM=3)\n", "\n", " # We need to rerun a few of these functions with the reset lists to make sure these functions\n", " # don't cheat by using analytic expressions\n", " GRHD.compute_sqrtgammaDET(gammaDD)\n", " GRHD.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha, betaU, gammaDD, ValenciavU)\n", " GRFFE.compute_smallb4U(gammaDD, betaU, alpha, GRHD.u4U_ito_ValenciavU, BU, sqrt4pi)\n", " GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)\n", " GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4U, GRFFE.smallbsquared,GRHD.u4U_ito_ValenciavU)\n", " GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDDdD,betaUdD,alphadD)\n", " GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD, GRFFE.TEM4UU)\n", " for i in range(3):\n", " desc = \"Adds the source term to StildeD\"+str(i)+\".\"\n", " name = \"calculate_StildeD\"+str(i)+\"_source_term\"\n", " outCfunction(\n", " outfile = os.path.join(outdir,name+\".h\"), desc=desc, name=name,\n", " params =\"const paramstruct *params,const REAL *auxevol_gfs, REAL *rhs_gfs\",\n", " body = general_access \\\n", " +metric_deriv_access[i]\\\n", " +outputC(GRHD.S_tilde_source_termD[i],\"Stilde_rhsD\"+str(i),\"returnstring\",params=outCparams)\\\n", " +write_final_quantity[i],\n", " loopopts =\"InteriorPoints\",\n", " rel_path_to_Cparams=os.path.join(\"../\"))\n", "\n", "def add_to_Cfunction_dict__functions_for_StildeD_source_term(outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi,\n", " includes=None, rel_path_to_Cparams=os.path.join(\"../\"),\n", " path_from_rootsrcdir_to_this_Cfunc=os.path.join(\"RHSs/\")):\n", " generate_memory_access_code()\n", " # First, we declare some dummy tensors that we will use for the codegen.\n", " gammaDDdD = ixp.declarerank3(\"gammaDDdD\",\"sym01\",DIM=3)\n", " betaUdD = ixp.declarerank2(\"betaUdD\",\"nosym\",DIM=3)\n", " alphadD = ixp.declarerank1(\"alphadD\",DIM=3)\n", "\n", " # We need to rerun a few of these functions with the reset lists to make sure these functions\n", " # don't cheat by using analytic expressions\n", " GRHD.compute_sqrtgammaDET(gammaDD)\n", " GRHD.u4U_in_terms_of_ValenciavU__rescale_ValenciavU_by_applying_speed_limit(alpha, betaU, gammaDD, ValenciavU)\n", " GRFFE.compute_smallb4U(gammaDD, betaU, alpha, GRHD.u4U_ito_ValenciavU, BU, sqrt4pi)\n", " GRFFE.compute_smallbsquared(gammaDD, betaU, alpha, GRFFE.smallb4U)\n", " GRFFE.compute_TEM4UU(gammaDD,betaU,alpha, GRFFE.smallb4U, GRFFE.smallbsquared,GRHD.u4U_ito_ValenciavU)\n", " GRHD.compute_g4DD_zerotimederiv_dD(gammaDD,betaU,alpha, gammaDDdD,betaUdD,alphadD)\n", " GRHD.compute_S_tilde_source_termD(alpha, GRHD.sqrtgammaDET,GRHD.g4DD_zerotimederiv_dD, GRFFE.TEM4UU)\n", " for i in range(3):\n", " desc = \"Adds the source term to StildeD\"+str(i)+\".\"\n", " name = \"calculate_StildeD\"+str(i)+\"_source_term\"\n", " params =\"const paramstruct *params,const REAL *auxevol_gfs, REAL *rhs_gfs\"\n", " body = general_access \\\n", " +metric_deriv_access[i]\\\n", " +outputC(GRHD.S_tilde_source_termD[i],\"Stilde_rhsD\"+str(i),\"returnstring\",params=outCparams)\\\n", " +write_final_quantity[i]\n", " loopopts =\"InteriorPoints\"\n", " add_to_Cfunction_dict(\n", " includes=includes,\n", " desc=desc,\n", " name=name, params=params,\n", " body=body, loopopts=loopopts,\n", " path_from_rootsrcdir_to_this_Cfunc=path_from_rootsrcdir_to_this_Cfunc,\n", " rel_path_to_Cparams=rel_path_to_Cparams)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: Code Validation against original C code \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation}$$\n", "\n", "To validate the code in this tutorial we check for agreement between the files\n", "\n", "1. that were written in this tutorial and\n", "1. those that are stored in `GiRaFFE_NRPy/GiRaFFE_Ccode_library` or generated by `GiRaFFE_NRPy_A2B.py`\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Output C function calculate_StildeD0_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_validation\\RHSs\\calculate_StildeD0_source_term.h\n", "Output C function calculate_StildeD1_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_validation\\RHSs\\calculate_StildeD1_source_term.h\n", "Output C function calculate_StildeD2_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_validation\\RHSs\\calculate_StildeD2_source_term.h\n", "Output C function calculate_StildeD0_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_library\\RHSs\\calculate_StildeD0_source_term.h\n", "Output C function calculate_StildeD1_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_library\\RHSs\\calculate_StildeD1_source_term.h\n", "Output C function calculate_StildeD2_source_term() to file GiRaFFE_NRPy\\GiRaFFE_Ccode_library\\RHSs\\calculate_StildeD2_source_term.h\n", "Printing difference between original C code and this code...\n", "Checking file calculate_StildeD0_source_term.h\n", "No difference. TEST PASSED!\n", "Checking file calculate_StildeD1_source_term.h\n", "No difference. TEST PASSED!\n", "Checking file calculate_StildeD2_source_term.h\n", "No difference. TEST PASSED!\n" ] } ], "source": [ "# Declare gridfunctions necessary to generate the C code:\n", "gammaDD = ixp.register_gridfunctions_for_single_rank2(\"AUXEVOL\",\"gammaDD\",\"sym01\",DIM=3)\n", "betaU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"betaU\",DIM=3)\n", "alpha = gri.register_gridfunctions(\"AUXEVOL\",\"alpha\",DIM=3)\n", "BU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"BU\",DIM=3)\n", "ValenciavU = ixp.register_gridfunctions_for_single_rank1(\"AUXEVOL\",\"ValenciavU\",DIM=3)\n", "StildeD = ixp.register_gridfunctions_for_single_rank1(\"EVOL\",\"StildeD\",DIM=3)\n", "# Declare this symbol:\n", "sqrt4pi = par.Cparameters(\"REAL\",thismodule,\"sqrt4pi\",\"sqrt(4.0*M_PI)\")\n", "\n", "# First, we generate the file using the functions written in this notebook:\n", "outCparams = \"outCverbose=False\"\n", "write_out_functions_for_StildeD_source_term(outdir,outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi)\n", "\n", "# Define the directory that we wish to validate against:\n", "valdir = os.path.join(\"GiRaFFE_NRPy\",\"GiRaFFE_Ccode_library\",\"RHSs\")\n", "cmd.mkdir(valdir)\n", "\n", "import GiRaFFE_NRPy.GiRaFFE_NRPy_Source_Terms as source\n", "source.write_out_functions_for_StildeD_source_term(valdir,outCparams,gammaDD,betaU,alpha,ValenciavU,BU,sqrt4pi)\n", "\n", "import difflib\n", "import sys\n", "\n", "print(\"Printing difference between original C code and this code...\")\n", "# Open the files to compare\n", "files = [\"calculate_StildeD0_source_term.h\",\"calculate_StildeD1_source_term.h\",\"calculate_StildeD2_source_term.h\"]\n", "\n", "for file in files:\n", " print(\"Checking file \" + file)\n", " with open(os.path.join(valdir,file)) as file1, open(os.path.join(outdir,file)) as file2:\n", " # Read the lines of each file\n", " file1_lines = file1.readlines()\n", " file2_lines = file2.readlines()\n", " num_diffs = 0\n", " for line in difflib.unified_diff(file1_lines, file2_lines, fromfile=os.path.join(valdir+file), tofile=os.path.join(outdir+file)):\n", " sys.stdout.writelines(line)\n", " num_diffs = num_diffs + 1\n", " if num_diffs == 0:\n", " print(\"No difference. TEST PASSED!\")\n", " else:\n", " print(\"ERROR: Disagreement found with .py file. See differences above.\")\n", " sys.exit(1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 3: 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-GiRaFFE_NRPy_C_code_library-Source_Terms](TTutorial-GiRaFFE_NRPy_C_code_library-Source_Terms.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": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Notebook output to PDF is only supported on Linux systems, with pdflatex installed.\n" ] } ], "source": [ "import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n", "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-GiRaFFE_NRPy-Source_Terms\")" ] } ], "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.4" } }, "nbformat": 4, "nbformat_minor": 2 }