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