{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"# Tutorial-IllinoisGRMHD: reconstruct_set_of_prims_PPM.C\n",
"\n",
"## Authors: Leo Werneck & Zach Etienne\n",
"\n",
"**This module is currently under development**\n",
"\n",
"## In this tutorial module we explain the piecewise parabolic method (PPM) used to reconstruct the primitive variables within IllinoisGRMHD\n",
"\n",
"### Required and recommended citations:\n",
"\n",
"* **(Required)** Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. ([arxiv:1501.07276](http://arxiv.org/abs/1501.07276)).\n",
"* **(Required)** Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) ([astro-ph/0512420](https://arxiv.org/abs/astro-ph/0512420)).\n",
"* **(Recommended)** Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 ([astro-ph/0210618](https://arxiv.org/abs/astro-ph/0210618))."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Table of Contents\n",
"$$\\label{toc}$$\n",
"\n",
"This module is organized as follows\n",
"\n",
"0. [Step 0](#src_dir): **Source directory creation**\n",
"1. [Step 1](#introduction): **Introduction**\n",
" 1. [Step 1.a](#ppm): *The Piecewise Parabolic Method (PPM)*\n",
" 1. [Step 1.b](#loop_defines_reconstruction): *The `loop_defines_reconstruction.h` header file*\n",
" 1. [Step 1.c](#preamble_reconstruct_set_of_prims_ppm): *The preamble to the `reconstruct_set_of_prims_PPM.C` code file*\n",
"1. [Step 2](#reconstruct_set_of_prims_ppm): **The `reconstruct_set_of_prims_PPM()` function**\n",
" 1. [Step 2.a](#reading_the_input_gfs): *Reading the input gridfunctions*\n",
" 1. [Step 2.b](#computation_of_du): *Evaluation of $\\delta U_{i}$*\n",
" 1. [Step 2.c](#computation_of_ur_and_ul): *Computing $U_{r}$ and $U_{l}$*\n",
" 1. [Step 2.d](#steepening_rhob): *Steepening $\\rho_{b}$*\n",
" 1. [Step 2.e](#flattening_and_monotonizing): *Flattening and monotonizing*\n",
" 1. [Step 2.f](#shifting_ur_and_ul): *Shifting $U_{r}$ and $U_{l}$*\n",
"1. [Step 3](#slope_limit): **The `slope_limit()` function**\n",
"1. [Step 4](#steepen_rho): **The `steepen_rho()` function**\n",
"1. [Step 5](#monotonize): **The `monotonize()` function**\n",
"1. [Step 6](#compute_p_cold__Gamma_cold): **The `compute_P_cold__Gamma_cold()` function**\n",
"1. [Step 7](#ftilde_gf_compute): **The `ftilde_gf_compute()` function**\n",
"1. [Step 8](#ftilde_compute): **The `ftilde_compute()` function**\n",
"1. [Step 9](#code_validation): **Code validation**\n",
" 1. [Step 9.a](#loop_defines_reconstruction__h_validation): *`loop_defines_reconstruction.h`*\n",
" 1. [Step 9.b](#reconstruct_set_of_prims_ppm__c_validation): *`reconstruct_set_of_prims_PPM.C`*\n",
"1. [Step 10](#latex_pdf_output): **Output this notebook to $\\LaTeX$-formatted PDF file**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 0: Source directory creation \\[Back to [top](#toc)\\]\n",
"$$\\label{src_dir}$$\n",
"\n",
"We will now use the [cmdline_helper.py NRPy+ module](Tutorial-Tutorial-cmdline_helper.ipynb) to create the source directory within the `IllinoisGRMHD` NRPy+ directory, if it does not exist yet."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Step 0: Creation of the IllinoisGRMHD source directory\n",
"# Step 0a: 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",
"# Step 0b: Load up cmdline_helper and create the directory\n",
"import cmdline_helper as cmd\n",
"IGM_src_dir_path = os.path.join(\"..\",\"src\")\n",
"cmd.mkdir(IGM_src_dir_path)\n",
"\n",
"# Step 0c: Create the output file path\n",
"outfile_path__reconstruct_set_of_prims_PPM__C = os.path.join(IGM_src_dir_path,\"reconstruct_set_of_prims_PPM.C\")\n",
"outfile_path__loop_defines_reconstruction__h = os.path.join(IGM_src_dir_path,\"loop_defines_reconstruction.h\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 1: Introduction \\[Back to [top](#toc)\\]\n",
"$$\\label{introduction}$$\n",
"\n",
"In this tutorial module we will go through the implementation of the piecewise parabolic method (PPM), introduced by [Colella & Woodward (1984)](https://crd.lbl.gov/assets/pubs_presos/AMCS/ANAG/A141984.pdf) (which shall henceforth be our main reference), used by `IllinoisGRMHD`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 1.a: The Piecewise Parabolic Method (PPM) \\[Back to [top](#toc)\\]\n",
"$$\\label{ppm}$$\n",
"\n",
"The piecewise parabolic method (PPM) is an algorithm used to construct the values of primitive variables, $U$, at cell interfaces. The interpolation procedure alone can lead to unstable evolutions. To remedy this, we introduce three different techniques:\n",
"\n",
"1. Steepening\n",
"1. Flatenning\n",
"1. Monotonizing\n",
"\n",
"These algorithms are intended to also produce narrower profiles near the vicinity of a shock. These steps tend to reduce the third-order accuracy of the interpolation code, but only in cases where a third-order interpolation algorithm would produce *worse* results (e.g. at local extrema).\n",
"\n",
"The algorithmic flow of the code is as follows:\n",
"\n",
"1. **Read the input**: determine which primitives are to be \"reconstructed\" (i.e. interpolated)\n",
"2. **Slope-limited gradient**: this must be computed for each of the primitives that will be reconstructed. We have:\n",
"\n",
"$$\n",
"\\boxed{\\delta U^{\\rm slope-lim} \\equiv\n",
"\\left\\{\n",
"\\begin{matrix}\n",
"{\\rm sign}\\left(\\delta_{m} U_{i}\\right)\\min\\left(\\left|\\delta_{m} U_{i}\\right|,c\\left|\\delta U_{i}\\right|,c\\left|\\delta U_{i+1}\\right|\\right) & ,\\ {\\rm if}\\ dU_{i}dU_{i+1} > 0\\\\\n",
"0 &,\\ {\\rm otherwise}\n",
"\\end{matrix}\n",
"\\right.}\\ ,\n",
"$$\n",
"\n",
" where $\\delta U^{\\rm slope-lim}$ is referred to as the slope-limited gradient of $U$ and\n",
"\n",
"\\begin{align}\n",
"\\delta U_{i} &\\equiv U_{i} - U_{i-1}\\ ,\\\\\n",
"\\delta_{m} U_{i} &\\equiv \\frac{\\delta U_{i} + \\delta U_{i+1}}{2} = \\frac{U_{i+1} - U_{i-1}}{2}\\ .\n",
"\\end{align}\n",
"\n",
"3. **Perform the interpolation**: We wish to determine $U_{r}$ and $U_{l}$, the values of $U$ at the cell interfaces. From the given set of known values $\\left\\{U_{i-2},U_{i-1},U_{i},U_{i+1},U_{i+2}\\right\\}$, one *interpolates* (with third-order accuracy) the values\n",
"\n",
"$$\n",
"\\begin{matrix}\n",
"U_{i+1/2} \\equiv U_{r,i} = U_{i+0} + \\frac{1}{2}\\left(U_{i+1} - U_{i+0}\\right) + \\frac{1}{6}\\left(\\delta U^{\\rm slope-lim}_{i+0} - \\delta U^{\\rm slope-lim}_{i+1}\\right)\\\\\n",
"U_{i-1/2} \\equiv U_{l,i} = U_{i-1} + \\frac{1}{2}\\left(U_{i+0} - U_{i-1}\\right) + \\frac{1}{6}\\left(\\delta U^{\\rm slope-lim}_{i-1} - \\delta U^{\\rm slope-lim}_{i+0}\\right)\n",
"\\end{matrix}\n",
"$$\n",
"\n",
"4. **Compute $P_{\\rm cold}$ and $\\Gamma_{\\rm cold}$**: In order to decide whether or not to apply the steepening procedure, we must evaluate the contact discontinuity condition:\n",
"\n",
"$$\n",
"\\boxed{\\Gamma_{\\rm cold} K_{0}\\frac{\\left|\\rho_{i+1}-\\rho_{i-1}\\right|}{\\min\\left(\\rho_{i+1},\\rho_{i-1}\\right)} \\geq \\frac{\\left|\\left(P_{\\rm cold}\\right)_{i+1}-\\left(P_{\\rm cold}\\right)_{i-1}\\right|}{\\min\\left[\\left(P_{\\rm cold}\\right)_{i+1},\\left(P_{\\rm cold}\\right)_{i-1}\\right]}}\\ ,\n",
"$$\n",
"\n",
"with $K_{0}$ a problem dependent constant.\n",
"\n",
"5. **Steepening**: *if necessary*, apply the steepening proceedure *only* to $U = \\rho_{b}^{\\color{red}\\ddagger}$. This involves performing the replacement\n",
"\n",
"$$\n",
"\\boxed{\n",
"\\begin{matrix}\n",
"\\rho_{r}\\to \\rho_{r}(1-\\eta) + \\rho^{\\rm MC}_{r}\\eta\\\\\n",
"\\rho_{l}\\to \\rho_{l}(1-\\eta) + \\rho^{\\rm MC}_{l}\\eta\n",
"\\end{matrix}\n",
"}\\ ,\n",
"$$\n",
"\n",
"where\n",
"\n",
"\\begin{align}\n",
"\\rho^{\\rm MC}_{r,i+1} &= \\rho_{i+1} - \\frac{1}{2}\\delta\\rho^{\\rm slope-lim}_{i+1}\\ ,\\\\\n",
"\\rho^{\\rm MC}_{l,i+0} &= \\rho_{i-1} + \\frac{1}{2}\\delta\\rho^{\\rm slope-lim}_{i-1}\\ ,\\\\\n",
"\\eta_{i} &= \\max\\left\\{0,\\min\\left[\\eta_{1}\\left(\\tilde\\eta_{i}-\\eta_{2}\\right),1\\right]\\right\\}\\ ,\n",
"\\end{align}\n",
"\n",
"with $\\eta_{1}$ and $\\eta_{2}$ constants and\n",
"\n",
"$$\n",
"\\tilde\\eta_{i} = \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"0&, \\ {\\rm if}\\ \\delta\\rho_{i} = 0\\ ,\\\\\n",
"-\\frac{1}{6}\\left(\\frac{\\delta^{2}\\rho_{i+1} - \\delta^{2}\\rho_{i-1}}{2\\delta\\rho_{i}}\\right)&, \\ {\\rm otherwise}\\ ,\n",
"\\end{matrix}\n",
"\\right.\n",
"$$\n",
"\n",
"and finally\n",
"\n",
"\\begin{align}\n",
"\\delta\\rho_{i+0} &= \\frac{\\rho_{i+1} - \\rho_{i-1}}{2}\\ ,\\\\\n",
"\\delta^{2}\\rho_{i-1} &= \\rho_{i+0} - 2\\rho_{i-1} + \\rho_{i-2}\\ ,\\\\\n",
"\\delta^{2}\\rho_{i+1} &= \\rho_{i+2} - 2\\rho_{i+1} + \\rho_{i+0}\\ .\n",
"\\end{align}\n",
"\n",
"$^{\\color{red}\\ddagger}$: note that the common notation is $\\rho_{0}$, but we will use $\\rho_{b}$ to represent the baryonic matter density.\n",
"\n",
"6. **Flattening**: *if necessary*, apply the flattening procedure to *all* primitives which are to be reconstructed. In the flattening procedure we modify either $U_{r}$, $U_{l}$, or both of them, according to\n",
"\n",
"$$\n",
"\\boxed{\n",
"\\begin{matrix}\n",
"U_{r,i+0} = U_{i+0}\\tilde{f} + U_{r,i+0}\\left(1-\\tilde{f}\\right)\\\\\n",
"U_{l,i+0} = U_{i+0}\\tilde{f} + U_{l,i+0}\\left(1-\\tilde{f}\\right)\n",
"\\end{matrix}\n",
"}\\ ,\n",
"$$\n",
"\n",
"where\n",
"\n",
"\\begin{align}\n",
"\\tilde{f} &= \\min\\left[1,w\\max\\left(0,q_{1}\\right)\\right]\\ ,\\\\\n",
"w &= \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"1\\ , &\\ {\\rm if}\\ q_{2} > \\epsilon_{2}\\ {\\rm and}\\ q_{2}\\left(v^{\\rm flux\\ dirn}_{i-1}-v^{\\rm flux\\ dirn}_{i+1}\\right)>0\\ \\left({\\rm inside\\ shock}\\right)\\ ,\\\\\n",
"0\\ , &\\ {\\rm otherwise}\\ \\left({\\rm outside\\ shock}\\right)\\ ,\n",
"\\end{matrix}\n",
"\\right.\n",
"\\end{align}\n",
"\n",
"and\n",
"\n",
"\\begin{align}\n",
"q_{1} &= \\left(\\frac{\\delta P_{1}}{\\delta P_{2}}-\\omega_{1}\\right)\\omega_{2}\\ ,\\\\\n",
"q_{2} &= \\frac{\\left|\\delta P_{1}\\right|}{\\min\\left(P_{i+1},P_{i-1}\\right)}\\ ,\n",
"\\end{align}\n",
"\n",
"with $\\omega_{1}$ and $\\omega_{2}$ constants and $\\delta P_{n} \\equiv P_{i+n} - P_{i-n}$.\n",
"\n",
"7. **Monotonizing**: *if necessary*, apply the monotonizing procedure to *all* primitives which are to be reconstructed. We check three different cases, modifying $U_{r}$ and $U_{l}$ as follows:\n",
"\n",
"$$\n",
"\\boxed{\n",
"\\begin{matrix}\n",
"\\text{Case 1: if}\\ \\left(U_{r} - U\\right)\\left(U - U_{l}\\right)\\leq 0 \\ \\text{then:}\\ \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"U_{r} \\to U\\ ,\\\\\n",
"\\ U_{l}\\to U\\ .\n",
"\\end{matrix}\n",
"\\right.\\\\\n",
"\\text{Case 2: if}\\ \\delta U\\left(U - \\delta_{m}U\\right) > \\frac{\\left(\\delta U\\right)^{2}}{6} \\ \\text{then:}\\ \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"U_{r} \\to U_{r}\\ ,\\\\\n",
"\\ U_{l}\\to 3U-2U_{r}\\ .\n",
"\\end{matrix}\n",
"\\right.\\\\\n",
"\\text{Case 3: if}\\ \\delta U\\left(U - \\delta_{m}U\\right) < -\\frac{\\left(\\delta U\\right)^{2}}{6} \\ \\text{then:}\\ \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"U_{r} \\to 3U-2U_{l}\\ ,\\\\\n",
"\\ U_{l}\\to U_{l}\\ .\n",
"\\end{matrix}\n",
"\\right.\n",
"\\end{matrix}\n",
"}\\ ,\n",
"$$\n",
"\n",
"where\n",
"\n",
"\\begin{align}\n",
"\\delta U &\\equiv U_{r} - U_{l}\\ ,\\\\\n",
"\\delta_{m}U &\\equiv \\frac{U_{r} + U_{l}}{2}\\ .\n",
"\\end{align}\n",
"\n",
"8. **Index shifting**: shift the indices of $U_{r}$ and $U_{l}$. We have, at this point\n",
"\n",
"\\begin{align}\n",
"U_{r,i} = U_{i+1/2}\\ ,\\\\\n",
"U_{l,i} = U_{i-1/2}\\ .\n",
"\\end{align}\n",
"\n",
"We then perform the following shift\n",
"\n",
"$$\n",
"\\boxed{\n",
"\\begin{matrix}\n",
"U_{i-1/2+\\epsilon} = U_{l,i}^{\\rm old} = U_{r,i}^{\\rm new}\\ ,\\\\\n",
"U_{i-1/2-\\epsilon} = U_{r,i-1}^{\\rm old} = U_{l,i}^{\\rm new}\\ .\\\\\n",
"\\end{matrix}\n",
"}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 1.b: The `loop_defines_reconstruction.h` header file \\[Back to [top](#toc)\\]\n",
"$$\\label{loop_defines_reconstruction}$$\n",
"\n",
"This header file defines useful quantities to be used throughout the `reconstruct_set_of_prims_PPM.C` code. They are:\n",
"\n",
"1. LOOP_DEFINE $\\rightarrow$ sets up a loop over the $x$, $y$, and $z$ directions, including the ghostzones\n",
"1. SET_INDEX_ARRAYS $\\rightarrow$ for a given direction, $x^{j}$, and coordinate range, $i\\in\\left[i_\\min,i_\\max\\right]$, finds the appropriate array indices\n",
"1. SET_INDEX_ARRAYS_3DBLOCK $\\rightarrow$ finds the appropriate array indices in all directions of a given range"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing ../src/loop_defines_reconstruction.h\n"
]
}
],
"source": [
"%%writefile $outfile_path__loop_defines_reconstruction__h\n",
"#ifndef LOOP_DEFINES_RECONSTRUCTION_H_\n",
"#define LOOP_DEFINES_RECONSTRUCTION_H_\n",
"\n",
"#define LOOP_DEFINE(gz_shift_lo,gz_shift_hi, ext,flux_dirn, ijkgz_lo_hi,gz_lo,gz_hi) \\\n",
" for(int rr=1;rr<=3;rr++) { \\\n",
" ijkgz_lo_hi[rr][0]= gz_lo[rr]; \\\n",
" ijkgz_lo_hi[rr][1]=ext[rr-1]-gz_hi[rr]; \\\n",
" } \\\n",
" ijkgz_lo_hi[flux_dirn][0] += gz_shift_lo; \\\n",
" ijkgz_lo_hi[flux_dirn][1] -= gz_shift_hi; \\\n",
" /* The following line is valid C99 */ \\\n",
" _Pragma(\"omp parallel for private(U,dU,slope_lim_dU,Ur,Ul)\") \\\n",
" for(int k=ijkgz_lo_hi[3][0];kmax_shift) CCTK_VError(VERR_DEF_PARAMS,\"FIX MAXNUMINDICES!\"); */ \\\n",
" int index_arr[4][MAXNUMINDICES]; \\\n",
" for(int idx=IMIN;idx<=IMAX;idx++) { \\\n",
" index_arr[flux_dirn][idx+max_shift]= \\\n",
" CCTK_GFINDEX3D(cctkGH, \\\n",
" i+idx*kronecker_delta[flux_dirn][0], \\\n",
" j+idx*kronecker_delta[flux_dirn][1], \\\n",
" k+idx*kronecker_delta[flux_dirn][2]); \\\n",
" }\n",
"\n",
"#define SET_INDEX_ARRAYS_3DBLOCK(IJKLOHI) \\\n",
" int max_shift=(MAXNUMINDICES/2); \\\n",
" int index_arr_3DB[MAXNUMINDICES][MAXNUMINDICES][MAXNUMINDICES]; \\\n",
" for(int idx_k=IJKLOHI[4];idx_k<=IJKLOHI[5];idx_k++) for(int idx_j=IJKLOHI[2];idx_j<=IJKLOHI[3];idx_j++) for(int idx_i=IJKLOHI[0];idx_i<=IJKLOHI[1];idx_i++) { \\\n",
" index_arr_3DB[idx_k+max_shift][idx_j+max_shift][idx_i+max_shift]=CCTK_GFINDEX3D(cctkGH,i+idx_i,j+idx_j,k+idx_k); \\\n",
" }\n",
"\n",
"#endif /* LOOP_DEFINES_RECONSTRUCTION_H_ */\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 1.c: The preamble to the `reconstruct_set_of_prims_PPM.C` code file \\[Back to [top](#toc)\\]\n",
"$$\\label{preamble_reconstruct_set_of_prims_ppm}$$\n",
"\n",
"We then initialize the `reconstruct_set_of_prims_PPM.C` code file with a basic preamble, some simple definitions, and function headers. Notice that these functions are defined in the other Steps of this tutorial notebook. See the [table of contents](#toc) above for more information."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"/*****************************************\n",
" * PPM Reconstruction Interface.\n",
" * Zachariah B. Etienne (2013)\n",
" *\n",
" * This version of PPM implements the standard\n",
" * Colella & Woodward PPM, though modified as in GRHydro\n",
" * to have 3 ghostzones instead of 4.\n",
" *****************************************/\n",
"\n",
"#define MINUS2 0\n",
"#define MINUS1 1\n",
"#define PLUS0 2\n",
"#define PLUS1 3\n",
"#define PLUS2 4\n",
"#define MAXNUMINDICES 5\n",
"// ^^^^^^^^^^^^^ Be _sure_ to define MAXNUMINDICES appropriately!\n",
"\n",
"// You'll find the #define's for LOOP_DEFINE and SET_INDEX_ARRAYS inside:\n",
"#include \"loop_defines_reconstruction.h\"\n",
"\n",
"static inline CCTK_REAL ftilde_compute(const int flux_dirn,CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES]);\n",
"static inline CCTK_REAL slope_limit(CCTK_REAL dU,CCTK_REAL dUp1);\n",
"static inline void steepen_rho(CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL slope_lim_dU[MAXNUMVARS][MAXNUMINDICES],\n",
" CCTK_REAL Gamma_th,CCTK_REAL P_cold,CCTK_REAL Gamma_cold,\n",
" CCTK_REAL *rho_br_ppm,CCTK_REAL *rho_bl_ppm);\n",
"static inline void compute_P_cold__Gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &Gamma_cold);\n",
"static inline void monotonize(CCTK_REAL U,CCTK_REAL &Ur,CCTK_REAL &Ul);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 2: The `reconstruct_set_of_prims_PPM()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{reconstruct_set_of_prims_ppm}$$\n",
"\n",
"The `reconstruct_set_of_prims_PPM()` function receives as input the direction in which to perform the reconstruction ($\\rm flux\\_dirn$), the number of primitives which are to be reconstructed ($\\rm num\\_prims\\_to\\_reconstruct$), which primitives are to be reconstructed ($\\rm which\\_prims\\_to\\_reconstruct$), the equation of state ($\\rm eos$), and the array which stores the primitives ($\\rm in\\_prims$). The reconstructed primitive values are then stored in the output arrays $\\rm out\\_prims\\_r$ and $\\rm out\\_prims\\_l$.\n",
"\n",
"**Note**: you can find more information on the ETK specific parameters (such as ${\\rm cctkGH}$ and $\\rm cctk\\_lsh$) by looking at the **Cactus Variables** paragraph of section C1.6.2 of the [Einstein Toolkit UserGuide](https://einsteintoolkit.org/usersguide/UsersGuidech9.html#x13-81000C1.6).\n",
"\n",
"Notice that we will start a loop that runs from 0 to $\\rm num\\_prims\\_to\\_reconstruct$, meaning that the discussion here will apply to *each* of the primitives that we wish to reconstruct."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
"static void reconstruct_set_of_prims_PPM(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,const int num_prims_to_reconstruct,const int *which_prims_to_reconstruct,eos_struct &eos,\n",
" gf_and_gz_struct *in_prims,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l,CCTK_REAL *ftilde_gf, CCTK_REAL *temporary) {\n",
"\n",
" DECLARE_CCTK_PARAMETERS;\n",
"\n",
" CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES],dU[MAXNUMVARS][MAXNUMINDICES],slope_lim_dU[MAXNUMVARS][MAXNUMINDICES],\n",
" Ur[MAXNUMVARS][MAXNUMINDICES],Ul[MAXNUMVARS][MAXNUMINDICES];\n",
" int ijkgz_lo_hi[4][2];\n",
"\n",
" for(int ww=0;ww\n",
"\n",
"## Step 2.a: Reading the input gridfunctions \\[Back to [top](#toc)\\]\n",
"$$\\label{reading_the_input_gfs}$$\n",
"\n",
"We start by reading in the input primitive gridfunction and storing them to a variable $U$. Notice that for a given direction and a given point $i$, we will know: $\\left\\{U_{i-2},U_{i-1},U_{i},U_{i+1},U_{i+2}\\right\\}$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" // *** LOOP 1: Interpolate to Ur and Ul, which are face values ***\n",
" // You will find that Ur depends on U at MINUS1,PLUS0, PLUS1,PLUS2, and\n",
" // Ul depends on U at MINUS2,MINUS1,PLUS0,PLUS1.\n",
" // However, we define the below loop from MINUS2 to PLUS2. Why not split\n",
" // this up and get additional points? The reason is that later on,\n",
" // Ur and Ul depend on ftilde, which is defined from MINUS2 to PLUS2,\n",
" // so we would lose those points anyway.\n",
" LOOP_DEFINE(2,2, cctk_lsh,flux_dirn, ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {\n",
" SET_INDEX_ARRAYS(-2,2,flux_dirn);\n",
" /* *** LOOP 1a: READ INPUT *** */\n",
" // Read in a primitive at all gridpoints between m = MINUS2 & PLUS2, where m's direction is given by flux_dirn. Store to U.\n",
" for(int ii=MINUS2;ii<=PLUS2;ii++) U[whichvar][ii] = in_prims[whichvar].gf[index_arr[flux_dirn][ii]];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 2.b: Evaluation of $\\delta U_{i}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{computation_of_du}$$\n",
"\n",
"We will need $\\delta U_{i} \\equiv U_{i} - U_{i-1}$ in order to compute $U_{r}$ and $U_{l}$. We will then compute (notice the notation change $i\\to i+0$ so that it is easier to understand the C code that follows):\n",
"\n",
"\\begin{align}\n",
"\\delta U_{i-1} &= U_{i-1} - U_{i-2}\\ ,\\\\\n",
"\\delta U_{i+0} &= U_{i+0} - U_{i-1}\\ ,\\\\\n",
"\\delta U_{i+1} &= U_{i+1} - U_{i+0} \\ ,\\\\\n",
"\\delta U_{i+2} &= U_{i+2} - U_{i+1}\\ .\n",
"\\end{align}\n",
"\n",
"After evaluating the differences $\\delta U$, we compute the slope-limited $\\delta U$, $\\delta U^{\\rm slope-lim}$, using the [`slope_limit()` function](#slope_limit)."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" /* *** LOOP 1b: DO COMPUTATION *** */\n",
" /* First, compute simple dU = U(i) - U(i-1), where direction of i\n",
" * is given by flux_dirn, and U is a primitive variable:\n",
" * {rho_b,P,vx,vy,vz,Bx,By,Bz}. */\n",
" // Note that for Ur and Ul at i, we must compute dU(i-1),dU(i),dU(i+1),\n",
" // and dU(i+2)\n",
" dU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2];\n",
" dU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1];\n",
" dU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0];\n",
" dU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1];\n",
"\n",
" // Then, compute slope-limited dU, using MC slope limiter:\n",
" slope_lim_dU[whichvar][MINUS1]=slope_limit(dU[whichvar][MINUS1],dU[whichvar][PLUS0]);\n",
" slope_lim_dU[whichvar][PLUS0] =slope_limit(dU[whichvar][PLUS0], dU[whichvar][PLUS1]);\n",
" slope_lim_dU[whichvar][PLUS1] =slope_limit(dU[whichvar][PLUS1], dU[whichvar][PLUS2]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 2.c: Computing $U_{r}$ and $U_{l}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{computation_of_ur_and_ul}$$\n",
"\n",
"We now compute $U_{r}$ and $U_{l}$. Keep in mind that $U_{r,i} = U_{i+1/2}$, while $U_{l,i} = U_{i-1/2}$. The implemented equation follows eq. A1 in [Duez *et al.* (2005)](http://arxiv.org/pdf/astro-ph/0503420.pdf), but with the standardd PPM coefficient of $\\frac{1}{6}$ (i.e. eq. A1 with $\\frac{1}{8}\\to\\frac{1}{6}$). Keep in mind that we simplify the equation slightly before implementing it:\n",
"\n",
"\\begin{align}\n",
"U_{r,i+0} &= U_{i+0} + \\frac{1}{2}\\left(U_{i+1} - U_{i+0}\\right) + \\frac{1}{6}\\left(\\delta U^{\\rm slope-lim}_{i+0} - \\delta U^{\\rm slope-lim}_{i+1}\\right) \\Rightarrow \\boxed{U_{r,i+0} = \\frac{1}{2}\\left(U_{i+1} + U_{i+0}\\right) + \\frac{1}{6}\\left(\\delta U^{\\rm slope-lim}_{i+0} - \\delta U^{\\rm slope-lim}_{i+1}\\right)}\\ ,\\\\\n",
"U_{l,i+0} &= U_{i-1} + \\frac{1}{2}\\left(U_{i+0} - U_{i-1}\\right) + \\frac{1}{6}\\left(\\delta U^{\\rm slope-lim}_{i-1} - \\delta U^{\\rm slope-lim}_{i+0}\\right) \\Rightarrow \\boxed{U_{l,i+0} = \\frac{1}{2}\\left(U_{i+0} + U_{i-1}\\right) + \\frac{1}{6}\\left(\\delta U^{\\rm slope-lim}_{i-1} - \\delta U^{\\rm slope-lim}_{i+0}\\right)}\\ .\n",
"\\end{align}\n",
"\n",
"After this step, the values of $U_{r,l,i+0}$ are stored as outputs."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" // Finally, compute face values Ur and Ul based on the PPM prescription\n",
" // (Eq. A1 in http://arxiv.org/pdf/astro-ph/0503420.pdf, but using standard 1/6=(1.0/6.0) coefficient)\n",
" // Ur[PLUS0] represents U(i+1/2)\n",
" // We applied a simplification to the following line: Ur=U+0.5*(U(i+1)-U) + ... = 0.5*(U(i+1)+U) + ...\n",
" Ur[whichvar][PLUS0] = 0.5*(U[whichvar][PLUS1] + U[whichvar][PLUS0] ) + (1.0/6.0)*(slope_lim_dU[whichvar][PLUS0] - slope_lim_dU[whichvar][PLUS1]);\n",
" // Ul[PLUS0] represents U(i-1/2)\n",
" // We applied a simplification to the following line: Ul=U(i-1)+0.5*(U-U(i-1)) + ... = 0.5*(U+U(i-1)) + ...\n",
" Ul[whichvar][PLUS0] = 0.5*(U[whichvar][PLUS0] + U[whichvar][MINUS1]) + (1.0/6.0)*(slope_lim_dU[whichvar][MINUS1] - slope_lim_dU[whichvar][PLUS0]);\n",
"\n",
" /* *** LOOP 1c: WRITE OUTPUT *** */\n",
" // Store right face values to {rho_br,Pr,vxr,vyr,vzr,Bxr,Byr,Bzr},\n",
" // and left face values to {rho_bl,Pl,vxl,vyl,vzl,Bxl,Byl,Bzl}\n",
" out_prims_r[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ur[whichvar][PLUS0];\n",
" out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ul[whichvar][PLUS0];\n",
" }"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 2.d: Steepening $\\rho_{b}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{steepening_rhob}$$\n",
"\n",
"Following the procedure described in [Step 4](#steepen_rho), we will now steepen $\\rho_{b}$ using the [`steepen_rho()` function](#steepen_rho). Keep in mind that although we will loop over all primitives which are set for reconstruction, the steepening procedure is applied to $\\rho_{b}$ *only*."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" // *** LOOP 2: STEEPEN RHOB ***\n",
" // Note that this loop applies ONLY to RHOB.\n",
" if(whichvar==RHOB) {\n",
" LOOP_DEFINE(2,2, cctk_lsh,flux_dirn, ijkgz_lo_hi,in_prims[whichvar].gz_lo,in_prims[whichvar].gz_hi) {\n",
"\tSET_INDEX_ARRAYS(-2,2,flux_dirn);\n",
"\t// Set rho and P separately, since within this loop,\n",
"\t// 1) steepen_rho() depends on RHOB(MINUS2,MINUS1,PLUS0,PLUS1,PLUS2)\n",
"\n",
"\t// Read in all primitives between MINUS2 & PLUS2. Store to U.\n",
"\tfor(int ii=MINUS2;ii<=PLUS2;ii++) U[RHOB][ii] = in_prims[RHOB ].gf[index_arr[flux_dirn][ii]];\n",
"\tfor(int ii=MINUS1;ii<=PLUS1;ii++) U[PRESSURE][ii] = in_prims[PRESSURE].gf[index_arr[flux_dirn][ii]];\n",
"\tUr[RHOB][PLUS0] = out_prims_r[RHOB].gf[index_arr[flux_dirn][PLUS0]];\n",
"\tUl[RHOB][PLUS0] = out_prims_l[RHOB].gf[index_arr[flux_dirn][PLUS0]];\n",
"\n",
"\tdU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2];\n",
"\tdU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1];\n",
"\tdU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0];\n",
"\tdU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1];\n",
"\n",
"\tslope_lim_dU[whichvar][MINUS1]=slope_limit(dU[whichvar][MINUS1],dU[whichvar][PLUS0]);\n",
"\t//slope_lim_dU[whichvar][PLUS0] =slope_limit(dU[whichvar][PLUS0], dU[whichvar][PLUS1]);\n",
"\tslope_lim_dU[whichvar][PLUS1] =slope_limit(dU[whichvar][PLUS1], dU[whichvar][PLUS2]);\n",
"\n",
"\t// Steepen rho\n",
"\t// DEPENDENCIES: RHOB face values, RHOB(MINUS2,MINUS1,PLUS0,PLUS1,PLUS2), P(MINUS1,PLUS0,PLUS1), and slope_lim_dU[RHOB](MINUS1,PLUS1)\n",
"\tCCTK_REAL P_cold,Gamma_cold;\n",
"\tcompute_P_cold__Gamma_cold(U[RHOB][PLUS0],eos, P_cold,Gamma_cold);\n",
"\tsteepen_rho(U,slope_lim_dU, Gamma_th,P_cold,Gamma_cold, Ur[RHOB],Ul[RHOB]);\n",
"\n",
"\t// Output rho\n",
"\tout_prims_r[RHOB].gf[index_arr[flux_dirn][PLUS0]] = Ur[RHOB][PLUS0];\n",
"\tout_prims_l[RHOB].gf[index_arr[flux_dirn][PLUS0]] = Ul[RHOB][PLUS0];\n",
" }\n",
" }\n",
" }"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 2.e: Flattening and monotonizing \\[Back to [top](#toc)\\]\n",
"$$\\label{flattening_and_monotonizing}$$\n",
"\n",
"The flattening procedure modifies $U_{r}$ and $U_{l}$ via\n",
"\n",
"$$\n",
"\\boxed{\n",
"\\begin{matrix}\n",
"U_{r,i+0} = U_{i+0}\\tilde{f} + U_{r,i+0}\\left(1-\\tilde{f}\\right)\\\\\n",
"U_{l,i+0} = U_{i+0}\\tilde{f} + U_{l,i+0}\\left(1-\\tilde{f}\\right)\n",
"\\end{matrix}\n",
"}\\ ,\n",
"$$\n",
"\n",
"where $\\tilde{f}$ is computed by the `ftilde_compute()` function, described in [Step 8](#ftilde_compute)."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" /* ORIGINAL PPM REQUIRES AT LEAST 4 GHOSTZONES, which can add\n",
" * significantly to the size of AMR ref. boundaries.\n",
" * To reduce to 3 ghostzones, we comment the following lines out:\n",
" * if ((P[indexp1] - P[indexm1]) <= 0.0) {\n",
" * f = MAX(ftilde,ftilde_p1);\n",
" * } else {\n",
" * f = MAX(ftilde,ftilde_m1);\n",
" * }\n",
" */\n",
"\n",
" // *** LOOP 3: FLATTEN BASED ON FTILDE AND MONOTONIZE ***\n",
" for(int ww=0;ww\n",
"\n",
"## Step 2.f: Shifting $U_{r}$ and $U_{l}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{shifting_ur_and_ul}$$\n",
"\n",
"At this point, we have\n",
"\n",
"\\begin{align}\n",
"U_{r,i} &= U_{i+1/2}\\ ,\\\\\n",
"U_{l,i} &= U_{i-1/2}\\ .\n",
"\\end{align}\n",
"\n",
"To keep things consistent, we shift indices to get\n",
"\n",
"$$\n",
"\\boxed{\n",
"\\begin{matrix}\n",
"U_{i-1/2+\\epsilon} = U_{l,i}^{\\rm old} = U_{r,i}^{\\rm new}\\ ,\\\\\n",
"U_{i-1/2-\\epsilon} = U_{r,i-1}^{\\rm old} = U_{l,i}^{\\rm new}\\ ,\\\\\n",
"\\end{matrix}\n",
"}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
" // Ur depends on ftilde, which depends on points of U between MINUS2 and PLUS2\n",
" out_prims_r[whichvar].gz_lo[flux_dirn]+=2;\n",
" out_prims_r[whichvar].gz_hi[flux_dirn]+=2;\n",
" // Ul depends on ftilde, which depends on points of U between MINUS2 and PLUS2\n",
" out_prims_l[whichvar].gz_lo[flux_dirn]+=2;\n",
" out_prims_l[whichvar].gz_hi[flux_dirn]+=2;\n",
" }\n",
"\n",
" // *** LOOP 4: SHIFT Ur AND Ul ***\n",
" /* Currently face values are set so that\n",
" * a) Ur(i) represents U(i+1/2), and\n",
" * b) Ul(i) represents U(i-1/2)\n",
" * Here, we shift so that the indices are consistent:\n",
" * a) U(i-1/2+epsilon) = oldUl(i) = newUr(i)\n",
" * b) U(i-1/2-epsilon) = oldUr(i-1) = newUl(i)\n",
" * Note that this step is not strictly necessary if you keep\n",
" * track of indices when computing the flux. */\n",
" for(int ww=0;ww\n",
"\n",
"# Step 3: The `slope_limit()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{slope_limit}$$\n",
"\n",
"We will now show the definition of $\\delta U^{\\rm slope-lim}$. The reason why we introduce this slope-limited procedure is twofold. First, it leads to steeper representations of discontinuities, and second it guarantees that $U_{i+1/2}$ lies inside the range $\\left[U_{i},U_{i+1}\\right]$.\n",
"\n",
"We start by defining\n",
"\n",
"\\begin{align}\n",
"\\delta U_{i} &\\equiv U_{i} - U_{i-1}\\ ,\\\\\n",
"\\delta_{m} U_{i} &\\equiv \\frac{\\delta U_{i} + \\delta U_{i+1}}{2} = \\frac{U_{i+1} - U_{i-1}}{2}\\ .\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
"// Set SLOPE_LIMITER_COEFF = 2.0 for MC, 1 for minmod\n",
"#define SLOPE_LIMITER_COEFF 2.0\n",
"\n",
"//Eq. 60 in JOURNAL OF COMPUTATIONAL PHYSICS 123, 1-14 (1996)\n",
"// [note the factor of 2 missing in the |a_{j+1} - a_{j}| term].\n",
"// Recall that dU = U_{i} - U_{i-1}.\n",
"static inline CCTK_REAL slope_limit(CCTK_REAL dU,CCTK_REAL dUp1) {\n",
" if(dU*dUp1 > 0.0) {\n",
" //delta_m_U=0.5 * [ (u_(i+1)-u_i) + (u_i-u_(i-1)) ] = (u_(i+1) - u_(i-1))/2 <-- first derivative, second-order; this should happen most of the time (smooth flows)\n",
" CCTK_REAL delta_m_U = 0.5*(dU + dUp1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we implement the slope-limited $\\delta U$ as\n",
"\n",
"$$\n",
"\\boxed{\\delta U^{\\rm slope-lim} \\equiv\n",
"\\left\\{\n",
"\\begin{matrix}\n",
"{\\rm sign}\\left(\\delta_{m} U_{i}\\right)\\min\\left(\\left|\\delta_{m} U_{i}\\right|,c\\left|\\delta U_{i}\\right|,c\\left|\\delta U_{i+1}\\right|\\right) & ,\\ {\\rm if}\\ dU_{i}dU_{i+1} > 0\\\\\n",
"0 &,\\ {\\rm otherwise}\n",
"\\end{matrix}\n",
"\\right.}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
" // EXPLANATION OF BELOW LINE OF CODE.\n",
" // In short, sign_delta_a_j = sign(delta_m_U) = (0.0 < delta_m_U) - (delta_m_U < 0.0).\n",
" // If delta_m_U>0, then (0.0 < delta_m_U)==1, and (delta_m_U < 0.0)==0, so sign_delta_a_j=+1\n",
" // If delta_m_U<0, then (0.0 < delta_m_U)==0, and (delta_m_U < 0.0)==1, so sign_delta_a_j=-1\n",
" // If delta_m_U==0,then (0.0 < delta_m_U)==0, and (delta_m_U < 0.0)==0, so sign_delta_a_j=0\n",
" int sign_delta_m_U = (0.0 < delta_m_U) - (delta_m_U < 0.0);\n",
" //Decide whether to use 2nd order derivative or first-order derivative, limiting slope.\n",
" return sign_delta_m_U*MIN(fabs(delta_m_U),MIN(SLOPE_LIMITER_COEFF*fabs(dUp1),SLOPE_LIMITER_COEFF*fabs(dU)));\n",
" }\n",
" return 0.0;\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 4: The `steepen_rho()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{steepen_rho}$$\n",
"\n",
"The steepening procedure withing the PPM algorithm is applied only to $\\rho_{b}$. The idea here is to produce narrower profiles near the vicinity of a contact discontinuity.\n",
"\n",
"**A NOTE ON NOTATION**: in the discussion below we will refer to $\\rho$ as $\\rho$ to keep the notation a bit lighter. No confusion should arise from this since there is no other quantity $\\rho$ involved.\n",
"\n",
"We start the algorithm by computing\n",
"\n",
"\\begin{align}\n",
"\\delta\\rho_{i+0} &= \\frac{\\rho_{i+1} - \\rho_{i-1}}{2}\\ ,\\\\\n",
"\\delta^{2}\\rho_{i-1} &= \\rho_{i+0} - 2\\rho_{i-1} + \\rho_{i-2}\\ ,\\\\\n",
"\\delta^{2}\\rho_{i+1} &= \\rho_{i+2} - 2\\rho_{i+1} + \\rho_{i+0}\\ .\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
"// standard Colella-Woodward parameters:\n",
"// K0 = 0.1d0, eta1 = 20.0, eta2 = 0.05, epsilon = 0.01d0\n",
"#define K0 0.1\n",
"#define ETA1 20.0\n",
"#define ETA2 0.05\n",
"#define EPSILON 0.01\n",
"static inline void steepen_rho(CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL slope_lim_dU[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL Gamma_th,CCTK_REAL P_cold,CCTK_REAL Gamma_cold,\n",
" CCTK_REAL *rho_br_ppm,CCTK_REAL *rho_bl_ppm) {\n",
"\n",
" // Next compute centered differences d RHOB and d^2 RHOB\n",
" CCTK_REAL d1rho_b = 0.5*(U[RHOB][PLUS1] - U[RHOB][MINUS1]);\n",
" CCTK_REAL d2rho_b_m1 = U[RHOB][PLUS0] - 2.0*U[RHOB][MINUS1] + U[RHOB][MINUS2];\n",
" CCTK_REAL d2rho_b_p1 = U[RHOB][PLUS2] - 2.0*U[RHOB][PLUS1] + U[RHOB][PLUS0];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we evaluate\n",
"\n",
"$$\n",
"\\Gamma = \\left.\\left(\\frac{\\partial P}{\\partial\\rho}\\right)\\middle/\\left(\\frac{P}{\\rho}\\right)\\right. = \\Gamma_{\\rm th} + \\left(\\Gamma_{\\rm cold} - \\Gamma_{\\rm th}\\right)\\frac{P_{\\rm cold}}{P}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" // Compute effective Gamma = (partial P / partial rho0)_s /(P/rho0)\n",
" CCTK_REAL Gamma = Gamma_th + (Gamma_cold-Gamma_th)*P_cold/U[PRESSURE][PLUS0];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next the contact discontinuity condition, eq. (3.2) of [Colella & Woodward (1984)](https://crd.lbl.gov/assets/pubs_presos/AMCS/ANAG/A141984.pdf), is checked:\n",
"\n",
"$$\n",
"\\Gamma K_{0}\\frac{\\left|\\rho_{i+1}-\\rho_{i-1}\\right|}{\\min\\left(\\rho_{i+1},\\rho_{i-1}\\right)} \\geq \\frac{\\left|P_{i+1}-P_{i-1}\\right|}{\\min\\left(P_{i+1},P_{i-1}\\right)}\\ ,\n",
"$$\n",
"\n",
"where $K_{0}$ is a problem dependent constant. Keep in mind that we implement the quantity\n",
"\n",
"$$\n",
"\\boxed{{\\rm contact\\_discontinuity\\_check} \\equiv \\Gamma K_{0}\\left|\\rho_{i+1}-\\rho_{i-1}\\right|\\min\\left(P_{i+1},P_{i-1}\\right) - \\left|P_{i+1}-P_{i-1}\\right|\\min\\left(\\rho_{i+1},\\rho_{i-1}\\right)}\\ ,\n",
"$$\n",
"\n",
"and verify whether ${\\rm contact\\_discontinuity\\_check} \\geq 0$ to verify the discontinuity condition. We also define the quantities\n",
"\n",
"$$\n",
"\\boxed{{\\rm second\\_deriv\\_check} \\equiv - \\delta^{2}\\rho_{i-1}\\delta^{2}\\rho_{i+1}}\\ ,\n",
"$$\n",
"\n",
"and\n",
"\n",
"$$\n",
"\\boxed{{\\rm relative\\_change\\_check} \\equiv 2\\left|\\delta\\rho\\right| - \\epsilon\\min\\left(\\rho_{i+1},\\rho_{i-1}\\right)}\\ ,\n",
"$$\n",
"\n",
"where again $\\epsilon$ is a constant. The contact discontinuity condition is then satisfied when all three quantities inside boxes above are non-negative. When that is the case, we evaluate\n",
"\n",
"$$\n",
"\\boxed{\\eta_{i} = \\max\\left\\{0,\\min\\left[\\eta_{1}\\left(\\tilde\\eta_{i}-\\eta_{2}\\right),1\\right]\\right\\}}\\ ,\n",
"$$\n",
"\n",
"where $\\eta_{1}$ and $\\eta_{2}$ are constants and\n",
"\n",
"$$\n",
"\\boxed{\\tilde\\eta_{i} = \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"0&, \\ {\\rm if}\\ \\delta\\rho_{i} = 0\\ ,\\\\\n",
"-\\frac{1}{6}\\left(\\frac{\\delta^{2}\\rho_{i+1} - \\delta^{2}\\rho_{i-1}}{2\\delta\\rho_{i}}\\right)&, \\ {\\rm otherwise}\\ .\n",
"\\end{matrix}\n",
"\\right.}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
" CCTK_REAL contact_discontinuity_check = Gamma*K0*fabs(U[RHOB][PLUS1]-U[RHOB][MINUS1])*\n",
" MIN(U[PRESSURE][PLUS1],U[PRESSURE][MINUS1])\n",
" -fabs(U[PRESSURE][PLUS1]-U[PRESSURE][MINUS1])*MIN(U[RHOB][PLUS1],U[RHOB][MINUS1]);\n",
" CCTK_REAL second_deriv_check = -d2rho_b_p1*d2rho_b_m1;\n",
" CCTK_REAL relative_change_check = fabs(2.0*d1rho_b) - EPSILON*MIN(U[RHOB][PLUS1],U[RHOB][MINUS1]);\n",
"\n",
" if(contact_discontinuity_check >= 0.0 && second_deriv_check >= 0.0\n",
" && relative_change_check >= 0.0) {\n",
"\n",
" CCTK_REAL eta_tilde=0.0;\n",
" if (fabs(d1rho_b) > 0.0) {\n",
" eta_tilde = -(1.0/6.0)*(d2rho_b_p1-d2rho_b_m1)/(2.0*d1rho_b);\n",
" }\n",
" CCTK_REAL eta = MAX(0.0,MIN(ETA1*(eta_tilde - ETA2),1.0));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We then apply the monotonized central (MC) scheme of [van Leer](https://www.sciencedirect.com/science/article/pii/002199917790095X) (see also [Step 3](#slope_limit) for a discussion on the quantities $\\delta\\rho^{\\rm slope-lim}_{i}$ below),\n",
"\n",
"\\begin{align}\n",
"\\rho^{\\rm MC}_{r,i+1} &= \\rho_{i+1} - \\frac{1}{2}\\delta\\rho^{\\rm slope-lim}_{i+1}\\ ,\\\\\n",
"\\rho^{\\rm MC}_{l,i+0} &= \\rho_{i-1} + \\frac{1}{2}\\delta\\rho^{\\rm slope-lim}_{i-1}\\ ,\n",
"\\end{align}\n",
"\n",
"so that, finally, the steepening algorithm sets\n",
"\n",
"$$\n",
"\\begin{matrix}\n",
"\\rho_{r}\\to \\rho_{r}(1-\\eta) + \\rho^{\\rm MC}_{r}\\eta\\ ,\\\\\n",
"\\rho_{l}\\to \\rho_{l}(1-\\eta) + \\rho^{\\rm MC}_{l}\\eta\\ ,\n",
"\\end{matrix}\n",
"$$\n",
"\n",
"or, as implemented below:\n",
"\n",
"$$\n",
"\\boxed{\\begin{matrix}\n",
"\\rho_{r,i+0}\\rightarrow \\rho_{r,i+0}\\left(1-\\eta_{i+0}\\right) + \\rho^{\\rm MC}_{r,i+1}\\eta_{i+0}\\ ,\\\\\n",
"\\rho_{l,i+0}\\rightarrow \\rho_{r,i+0}\\left(1-\\eta_{i+0}\\right) + \\rho^{\\rm MC}_{l,i+0}\\eta_{i+0}\\ .\n",
"\\end{matrix}}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
" // Next compute Urp1 and Ul for RHOB, using the MC prescription:\n",
" // Ur_p1 = U_p1 - 0.5*slope_lim_dU_p1\n",
" CCTK_REAL rho_br_mc_p1 = U[RHOB][PLUS1] - 0.5*slope_lim_dU[RHOB][PLUS1];\n",
" // Ul = U_m1 + 0.5*slope_lim_dU_m1\n",
" // Based on this line of code, Ur[index] = a_j - \\delta_m a_j / 2. (cf. Eq. 65 in Marti & Muller's \"PPM Method for 1D Relativistic Hydro.\" paper)\n",
" // So: Ur[indexp1] = a_{j+1} - \\delta_m a_{j+1} / 2. This is why we have rho_br_mc[indexp1]\n",
" CCTK_REAL rho_bl_mc = U[RHOB][MINUS1] + 0.5*slope_lim_dU[RHOB][MINUS1];\n",
"\n",
" rho_bl_ppm[PLUS0] = rho_bl_ppm[PLUS0]*(1.0-eta) + rho_bl_mc*eta;\n",
" rho_br_ppm[PLUS0] = rho_br_ppm[PLUS0]*(1.0-eta) + rho_br_mc_p1*eta;\n",
"\n",
" }\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 5: The `monotonize()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{monotonize}$$\n",
"\n",
"The value $U_{i+1/2}$ will be assigned to $U_{l,i}$ and $U_{r,i-1}$ for most values of $i$, but in some cases this would lead to incorrect interpolation results. Near discontinuities, the value of either $U_{l}$, $U_{r}$, or both needs to be adjusted.\n",
"\n",
"Consider, then, the following quantities:\n",
"\n",
"\\begin{align}\n",
"\\delta U &\\equiv U_{r} - U_{l}\\ ,\\\\\n",
"\\delta_{m}U &\\equiv \\frac{U_{r} + U_{l}}{2}\\ .\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
"static inline void monotonize(CCTK_REAL U,CCTK_REAL &Ur,CCTK_REAL &Ul) {\n",
" CCTK_REAL dU = Ur - Ul;\n",
" CCTK_REAL mU = 0.5*(Ur+Ul);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then, following eq. (1.10) of [Colella & Woodward (1984)](https://crd.lbl.gov/assets/pubs_presos/AMCS/ANAG/A141984.pdf), we will check the following three cases:\n",
"\n",
"$$\n",
"\\boxed{\\text{Case 1: if}\\ \\left(U_{r} - U\\right)\\left(U - U_{l}\\right)\\leq 0 \\ \\text{then:}\\ \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"U_{r} \\to U\\ ,\\\\\n",
"\\ U_{l}\\to U\\ .\n",
"\\end{matrix}\n",
"\\right.}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" if ( (Ur-U)*(U-Ul) <= 0.0) {\n",
" Ur = U;\n",
" Ul = U;\n",
" return;\n",
" }"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\boxed{\\text{Case 2: if}\\ \\delta U\\left(U - \\delta_{m}U\\right) > \\frac{\\left(\\delta U\\right)^{2}}{6} \\ \\text{then:}\\ \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"U_{r} \\to U_{r}\\ ,\\\\\n",
"\\ U_{l}\\to 3U-2U_{r}\\ .\n",
"\\end{matrix}\n",
"\\right.}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
" if ( dU*(U-mU) > (1.0/6.0)*SQR(dU)) {\n",
" Ul = 3.0*U - 2.0*Ur;\n",
" return;\n",
" }"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\boxed{\\text{Case 3: if}\\ \\delta U\\left(U - \\delta_{m}U\\right) < -\\frac{\\left(\\delta U\\right)^{2}}{6} \\ \\text{then:}\\ \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"U_{r} \\to 3U-2U_{l}\\ ,\\\\\n",
"\\ U_{l}\\to U_{l}\\ .\n",
"\\end{matrix}\n",
"\\right.}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
" if ( dU*(U-mU) < -(1.0/6.0)*SQR(dU)) {\n",
" Ur = 3.0*U - 2.0*Ul;\n",
" return;\n",
" }\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 6: The `compute_P_cold__Gamma_cold()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{compute_p_cold__Gamma_cold}$$\n",
"\n",
"This part of the code evaluates $P_{\\rm cold}$ and $\\Gamma_{\\rm cold}$ for the equations of state (EOS) presented in eqs. 13-16 of [Stephens *et al.* (2008)](http://arxiv.org/pdf/0802.0200.pdf).\n",
"\n",
"First, if $\\rho_{b} = 0$, then $P_{\\rm cold} = 0$ and $\\Gamma_{\\rm cold}$ simply receives its tabulated value."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
"static inline void compute_P_cold__Gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &Gamma_cold) {\n",
" // This code handles equations of state of the form defined\n",
" // in Eqs 13-16 in http://arxiv.org/pdf/0802.0200.pdf\n",
"\n",
" // Default in case rho_b == 0.0\n",
" if(rho_b==0.0) { P_cold = 0.0; Gamma_cold = eos.Gamma_ppoly_tab[0]; return; }"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we consider the case where the EOS is given by a single-polytrope\n",
"\n",
"$$\n",
"\\boxed{P_{\\rm cold} = \\kappa \\rho_{b}^{\\Gamma_{\\rm cold}}}\\ ,\n",
"$$\n",
"\n",
"and also the piecewise polytrope EOS\n",
"\n",
"$$\n",
"\\boxed{\n",
"P_{\\rm cold} =\n",
"\\left\\{\n",
"\\begin{matrix}\n",
"K_{0}\\rho^{\\Gamma_{0}} & , & \\rho \\leq \\rho_{0}\\\\\n",
"K_{1}\\rho^{\\Gamma_{1}} & , & \\rho_{0} \\leq \\rho \\leq \\rho_{1}\\\\\n",
"\\vdots & & \\vdots\\\\\n",
"K_{j}\\rho^{\\Gamma_{j}} & , & \\rho_{j-1} \\leq \\rho \\leq \\rho_{j}\\\\\n",
"\\vdots & & \\vdots\\\\\n",
"K_{N-1}\\rho^{\\Gamma_{N-1}} & , & \\rho_{N-2} \\leq \\rho \\leq \\rho_{N-1}\\\\\n",
"K_{N}\\rho^{\\Gamma_{N}} & , & \\rho \\geq \\rho_{N-1}\n",
"\\end{matrix}\n",
"\\right.\n",
"}\\ .\n",
"$$\n",
"\n",
"Notice that we left the fact that $\\Gamma_{i} \\equiv \\Gamma_{{\\rm cold},i}$ implicit above."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
" /***********************************\n",
" * Piecewise Polytropic EOS Patch *\n",
" * Computing P_cold and Gamma_cold *\n",
" ***********************************/\n",
" int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b);\n",
" Gamma_cold = eos.Gamma_ppoly_tab[polytropic_index];\n",
" P_cold = eos.K_ppoly_tab[polytropic_index]*pow(rho_b,Gamma_cold);\n",
"\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 7: The `ftilde_gf_compute()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{ftilde_gf_compute}$$\n",
"\n",
"This is the driver function of the `ftilde_compute()` function, setting up the dependencies needed to compute $\\tilde{f}$. Please refer to the next step to see how $\\tilde{f}$ is computed."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
"#define OMEGA1 0.75\n",
"#define OMEGA2 10.0\n",
"#define EPSILON2 0.33\n",
"static void ftilde_gf_compute(const cGH *cctkGH,const int *cctk_lsh,const int flux_dirn,gf_and_gz_struct *in_prims,CCTK_REAL *ftilde_gf) {\n",
" int ijkgz_lo_hi[4][2];\n",
" CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES];\n",
" /*Remove gcc unused variable warning/error Re: Pragma statement in loop define:*/\n",
" CCTK_REAL dU,slope_lim_dU,Ur,Ul; dU=slope_lim_dU=Ur=Ul=0.0; dU*=0;\n",
" // Compute ftilde, which is used for flattening left and right face values\n",
" LOOP_DEFINE(2,2, cctk_lsh,flux_dirn, ijkgz_lo_hi,in_prims[VX+(flux_dirn-1)].gz_lo,in_prims[VX+(flux_dirn-1)].gz_hi) {\n",
" SET_INDEX_ARRAYS(-2,2,flux_dirn);\n",
" for(int ii=MINUS2;ii<=PLUS2;ii++) U[PRESSURE][ii] = in_prims[PRESSURE].gf[index_arr[flux_dirn][ii]];\n",
" U[VX+(flux_dirn-1)][MINUS1] = in_prims[VX+(flux_dirn-1)].gf[index_arr[flux_dirn][MINUS1]];\n",
" U[VX+(flux_dirn-1)][PLUS1] = in_prims[VX+(flux_dirn-1)].gf[index_arr[flux_dirn][PLUS1]];\n",
"\n",
" // Compute ftilde, which is used for flattening left and right face values\n",
" // DEPENDENCIES: P(MINUS2,MINUS1,PLUS1,PLUS2) and v^m(MINUS1,PLUS1), where m=flux_dirn={1,2,3}={x,y,z}.\n",
" ftilde_gf[index_arr[flux_dirn][PLUS0]] = ftilde_compute(flux_dirn,U);\n",
" }\n",
"}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 8: The `ftilde_compute()` function \\[Back to [top](#toc)\\]\n",
"$$\\label{ftilde_compute}$$\n",
"\n",
"We start by evaluating\n",
"\n",
"\\begin{align}\n",
"\\delta P_{1} &\\equiv P_{i+1} - P_{i-1}\\ ,\\\\\n",
"\\delta P_{2} &\\equiv P_{i+2} - P_{i-2}\\ .\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
"static inline CCTK_REAL ftilde_compute(const int flux_dirn,CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES]) {\n",
" CCTK_REAL dP1 = U[PRESSURE][PLUS1] - U[PRESSURE][MINUS1];\n",
" CCTK_REAL dP2 = U[PRESSURE][PLUS2] - U[PRESSURE][MINUS2];"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we modify the standard PPM algorithm slightly by introducing the following conditions:\n",
"\n",
"\\begin{align}\n",
"{\\rm if}\\ \\left|\\frac{\\delta P_{1}}{\\delta_{m}P_{1}}\\right| = 0\\ {\\rm then\\ set}\\ \\delta P_{1}=0\\ ,\\\\\n",
"{\\rm if}\\ \\left|\\frac{\\delta P_{2}}{\\delta_{m}P_{2}}\\right| = 0\\ {\\rm then\\ set}\\ \\delta P_{2}=0\\ ,\n",
"\\end{align}\n",
"\n",
"where\n",
"\n",
"\\begin{align}\n",
"\\delta_{m} P_{1} &\\frac{\\equiv P_{i+1} + P_{i-1}}{2}\\ ,\\\\\n",
"\\delta_{m} P_{2} &\\frac{\\equiv P_{i+2} + P_{i-2}}{2}\\ .\n",
"\\end{align}\n",
"\n",
"Note that if the first condition above is satisfied then we are *not* inside a shock, while if the second condition is triggered *alone* there *may* be a shock."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" // MODIFICATION TO STANDARD PPM:\n",
" // Cure roundoff error issues when dP1==0 or dP2==0 to 15 or more significant digits.\n",
" CCTK_REAL avg1=0.5*(U[PRESSURE][PLUS1] + U[PRESSURE][MINUS1]);\n",
" CCTK_REAL avg2=0.5*(U[PRESSURE][PLUS2] + U[PRESSURE][MINUS2]);\n",
" if(fabs(dP1)/avg1<1e-15) dP1=0.0; /* If this is triggered, there is NO shock */\n",
" if(fabs(dP2)/avg2<1e-15) dP2=0.0; /* If this is triggered alone, there may be a shock. Otherwise if triggered with above, NO shock. */"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we set\n",
"\n",
"$$\n",
"{\\rm dP1\\_over\\_dP2} = \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"\\frac{\\delta P_{1}}{\\delta P_{2}} &,\\ {\\rm if}\\ \\delta P_{2} \\neq 0\\\\\n",
"1 &,\\ {\\rm otherwise}\n",
"\\end{matrix}\n",
"\\right.\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" CCTK_REAL dP1_over_dP2=1.0;\n",
" if (dP2 != 0.0) dP1_over_dP2 = dP1/dP2;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We then construct\n",
"\n",
"\\begin{align}\n",
"q_{1} &= \\left({\\rm dP1\\_over\\_dP2}-\\omega_{1}\\right)\\omega_{2}\\ ,\\\\\n",
"q_{2} &= \\frac{\\left|\\delta P_{1}\\right|}{\\min\\left(P_{i+1},P_{i-1}\\right)}\\ .\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" CCTK_REAL q1 = (dP1_over_dP2-OMEGA1)*OMEGA2;\n",
" CCTK_REAL q2 = fabs(dP1)/MIN(U[PRESSURE][PLUS1],U[PRESSURE][MINUS1]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We then initialize a new variable, $w$, to\n",
"\n",
"$$\n",
"w = \n",
"\\left\\{\n",
"\\begin{matrix}\n",
"1\\ , &\\ {\\rm if}\\ q_{2} > \\epsilon_{2}\\ {\\rm and}\\ q_{2}\\left(v^{\\rm flux\\ dirn}_{i-1}-v^{\\rm flux\\ dirn}_{i+1}\\right)>0\\ \\left({\\rm inside\\ shock}\\right)\\ ,\\\\\n",
"0\\ , &\\ {\\rm otherwise}\\ \\left({\\rm outside\\ shock}\\right)\\ ,\n",
"\\end{matrix}\n",
"\\right.\n",
"$$\n",
"\n",
"where $v^{\\rm flux\\ dirn}$ represents either $v^{x}$, $v^{y}$, or $v^{z}$, dependding on the flux direction. Finally, we get\n",
"\n",
"$$\n",
"\\boxed{\\tilde{f} = \\min\\left[1,w\\max\\left(0,q_{1}\\right)\\right]}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/reconstruct_set_of_prims_PPM.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"\n",
" // w==0 -> NOT inside a shock\n",
" CCTK_REAL w=0.0;\n",
"\n",
" // w==1 -> inside a shock\n",
" if (q2 > EPSILON2 && q2*( (U[VX+(flux_dirn-1)][MINUS1]) - (U[VX+(flux_dirn-1)][PLUS1]) ) > 0.0) w = 1.0;\n",
"\n",
" return MIN(1.0, w*MAX(0.0,q1));\n",
"}\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 9: Code validation \\[Back to [top](#toc)\\]\n",
"$$\\label{code_validation}$$\n",
"\n",
"First we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 9.a: `loop_defines_reconstruction.h` \\[Back to [top](#toc)\\]\n",
"$$\\label{loop_defines_reconstruction__h_validation}$$"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Validation test for loop_defines_reconstruction.h: FAILED!\n",
"Diff:\n",
"17c17\n",
"< // This define only sets indices. \n",
"---\n",
"> // This define only sets indices.\n",
"39a40\n",
"> \n"
]
}
],
"source": [
"# Verify if the code generated by this tutorial module\n",
"# matches the original IllinoisGRMHD source code\n",
"\n",
"# First download the original IllinoisGRMHD source code\n",
"import urllib\n",
"from os import path\n",
"\n",
"original_IGM_file_url = \"https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/loop_defines_reconstruction.h\"\n",
"original_IGM_file_name = \"loop_defines_reconstruction-original.h\"\n",
"original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)\n",
"\n",
"# Then download the original IllinoisGRMHD source code\n",
"# We try it here in a couple of ways in an attempt to keep\n",
"# the code more portable\n",
"try:\n",
" original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
" # Write down the file the original IllinoisGRMHD source code\n",
" with open(original_IGM_file_path,\"w\") as file:\n",
" file.write(original_IGM_file_code)\n",
"except:\n",
" try:\n",
" original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
" # Write down the file the original IllinoisGRMHD source code\n",
" with open(original_IGM_file_path,\"w\") as file:\n",
" file.write(original_IGM_file_code)\n",
" except:\n",
" # If all else fails, hope wget does the job\n",
" !wget -O $original_IGM_file_path $original_IGM_file_url\n",
"\n",
"# Perform validation\n",
"Validation__loop_defines_reconstruction__h = !diff $original_IGM_file_path $outfile_path__loop_defines_reconstruction__h\n",
"\n",
"if Validation__loop_defines_reconstruction__h == []:\n",
" # If the validation passes, we do not need to store the original IGM source code file\n",
" !rm $original_IGM_file_path\n",
" print(\"Validation test for loop_defines_reconstruction.h: PASSED!\")\n",
"else:\n",
" # If the validation fails, we keep the original IGM source code file\n",
" print(\"Validation test for loop_defines_reconstruction.h: FAILED!\")\n",
" # We also print out the difference between the code generated\n",
" # in this tutorial module and the original IGM source code\n",
" print(\"Diff:\")\n",
" for diff_line in Validation__loop_defines_reconstruction__h:\n",
" print(diff_line)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"## Step 9.b: `reconstruct_set_of_prims_PPM.C` \\[Back to [top](#toc)\\]\n",
"$$\\label{reconstruct_set_of_prims_ppm__c_validation}$$"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Validation test for reconstruct_set_of_prims_PPM.C: FAILED!\n",
"Diff:\n",
"5c5\n",
"< * This version of PPM implements the standard \n",
"---\n",
"> * This version of PPM implements the standard\n",
"7c7\n",
"< * to have 3 ghostzones instead of 4. \n",
"---\n",
"> * to have 3 ghostzones instead of 4.\n",
"24c24\n",
"< CCTK_REAL gamma_th,CCTK_REAL P_cold,CCTK_REAL gamma_cold,\n",
"---\n",
"> CCTK_REAL Gamma_th,CCTK_REAL P_cold,CCTK_REAL Gamma_cold,\n",
"26c26\n",
"< static inline void compute_P_cold__gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &gamma_cold);\n",
"---\n",
"> static inline void compute_P_cold__Gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &Gamma_cold);\n",
"28a29\n",
"> \n",
"31a33,34\n",
"> DECLARE_CCTK_PARAMETERS;\n",
"> \n",
"44c47,48\n",
"< \n",
"---\n",
"> \n",
"> \n",
"55c59\n",
"< // Read in a primitive at all gridpoints between m = MINUS2 & PLUS2, where m's direction is given by flux_dirn. Store to U. \n",
"---\n",
"> // Read in a primitive at all gridpoints between m = MINUS2 & PLUS2, where m's direction is given by flux_dirn. Store to U.\n",
"57c61,62\n",
"< \n",
"---\n",
"> \n",
"> \n",
"59,60c64,65\n",
"< /* First, compute simple dU = U(i) - U(i-1), where direction of i \n",
"< * is given by flux_dirn, and U is a primitive variable: \n",
"---\n",
"> /* First, compute simple dU = U(i) - U(i-1), where direction of i\n",
"> * is given by flux_dirn, and U is a primitive variable:\n",
"62c67\n",
"< // Note that for Ur and Ul at i, we must compute dU(i-1),dU(i),dU(i+1), \n",
"---\n",
"> // Note that for Ur and Ul at i, we must compute dU(i-1),dU(i),dU(i+1),\n",
"64,67c69,72\n",
"< dU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2]; \n",
"< dU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1]; \n",
"< dU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0]; \n",
"< dU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1]; \n",
"---\n",
"> dU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2];\n",
"> dU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1];\n",
"> dU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0];\n",
"> dU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1];\n",
"74c79,80\n",
"< // Finally, compute face values Ur and Ul based on the PPM prescription \n",
"---\n",
"> \n",
"> // Finally, compute face values Ur and Ul based on the PPM prescription\n",
"87c93\n",
"< out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ul[whichvar][PLUS0]; \n",
"---\n",
"> out_prims_l[whichvar].gf[index_arr[flux_dirn][PLUS0]] = Ul[whichvar][PLUS0];\n",
"89a96\n",
"> \n",
"104,107c111,114\n",
"< \tdU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2]; \n",
"< \tdU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1]; \n",
"< \tdU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0]; \n",
"< \tdU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1]; \n",
"---\n",
"> \tdU[whichvar][MINUS1] = U[whichvar][MINUS1]- U[whichvar][MINUS2];\n",
"> \tdU[whichvar][PLUS0] = U[whichvar][PLUS0] - U[whichvar][MINUS1];\n",
"> \tdU[whichvar][PLUS1] = U[whichvar][PLUS1] - U[whichvar][PLUS0];\n",
"> \tdU[whichvar][PLUS2] = U[whichvar][PLUS2] - U[whichvar][PLUS1];\n",
"115,117c122,124\n",
"< \tCCTK_REAL P_cold=0,gamma_cold=0;\n",
"< \tcompute_P_cold__gamma_cold(U[RHOB][PLUS0],eos, P_cold,gamma_cold);\n",
"< \tsteepen_rho(U,slope_lim_dU, eos.gamma_th,P_cold,gamma_cold, Ur[RHOB],Ul[RHOB]);\n",
"---\n",
"> \tCCTK_REAL P_cold,Gamma_cold;\n",
"> \tcompute_P_cold__Gamma_cold(U[RHOB][PLUS0],eos, P_cold,Gamma_cold);\n",
"> \tsteepen_rho(U,slope_lim_dU, Gamma_th,P_cold,Gamma_cold, Ur[RHOB],Ul[RHOB]);\n",
"125a133\n",
"> \n",
"142c150\n",
"< \n",
"---\n",
"> \n",
"152a161\n",
"> \n",
"158a168\n",
"> \n",
"160c170\n",
"< out_prims_r[whichvar].gz_lo[flux_dirn]+=2; \n",
"---\n",
"> out_prims_r[whichvar].gz_lo[flux_dirn]+=2;\n",
"185c195\n",
"< // Then shift so that Ur represents the gridpoint at i-1/2+epsilon, \n",
"---\n",
"> // Then shift so that Ur represents the gridpoint at i-1/2+epsilon,\n",
"195c205\n",
"< // As for Ur, we didn't need to get rid of another ghostzone, \n",
"---\n",
"> // As for Ur, we didn't need to get rid of another ghostzone,\n",
"202a213\n",
"> \n",
"206,207c217,218\n",
"< //Eq. 60 in JOURNAL OF COMPUTATIONAL PHYSICS 123, 1-14 (1996) \n",
"< // [note the factor of 2 missing in the |a_{j+1} - a_{j}| term]. \n",
"---\n",
"> //Eq. 60 in JOURNAL OF COMPUTATIONAL PHYSICS 123, 1-14 (1996)\n",
"> // [note the factor of 2 missing in the |a_{j+1} - a_{j}| term].\n",
"212a224\n",
"> \n",
"224a237\n",
"> \n",
"231c244\n",
"< static inline void steepen_rho(CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL slope_lim_dU[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL gamma_th,CCTK_REAL P_cold,CCTK_REAL gamma_cold,\n",
"---\n",
"> static inline void steepen_rho(CCTK_REAL U[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL slope_lim_dU[MAXNUMVARS][MAXNUMINDICES],CCTK_REAL Gamma_th,CCTK_REAL P_cold,CCTK_REAL Gamma_cold,\n",
"238a252\n",
"> \n",
"240,242c254,257\n",
"< CCTK_REAL Gamma = gamma_th + (gamma_cold-gamma_th)*P_cold/U[PRESSURE][PLUS0];\n",
"< CCTK_REAL contact_discontinuity_check = Gamma*K0*fabs(U[RHOB][PLUS1]-U[RHOB][MINUS1])* \n",
"< MIN(U[PRESSURE][PLUS1],U[PRESSURE][MINUS1]) \n",
"---\n",
"> CCTK_REAL Gamma = Gamma_th + (Gamma_cold-Gamma_th)*P_cold/U[PRESSURE][PLUS0];\n",
"> \n",
"> CCTK_REAL contact_discontinuity_check = Gamma*K0*fabs(U[RHOB][PLUS1]-U[RHOB][MINUS1])*\n",
"> MIN(U[PRESSURE][PLUS1],U[PRESSURE][MINUS1])\n",
"247c262\n",
"< if(contact_discontinuity_check >= 0.0 && second_deriv_check >= 0.0 \n",
"---\n",
"> if(contact_discontinuity_check >= 0.0 && second_deriv_check >= 0.0\n",
"254a270\n",
"> \n",
"263c279\n",
"< rho_bl_ppm[PLUS0] = rho_bl_ppm[PLUS0]*(1.0-eta) + rho_bl_mc*eta; \n",
"---\n",
"> rho_bl_ppm[PLUS0] = rho_bl_ppm[PLUS0]*(1.0-eta) + rho_bl_mc*eta;\n",
"268a285\n",
"> \n",
"272,273c289,291\n",
"< \n",
"< if ( (Ur-U)*(U-Ul) <= 0.0) { \n",
"---\n",
"> \n",
"> \n",
"> if ( (Ur-U)*(U-Ul) <= 0.0) {\n",
"278c296,297\n",
"< if ( dU*(U-mU) > (1.0/6.0)*SQR(dU)) { \n",
"---\n",
"> \n",
"> if ( dU*(U-mU) > (1.0/6.0)*SQR(dU)) {\n",
"281a301\n",
"> \n",
"288c308,309\n",
"< static inline void compute_P_cold__gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &gamma_cold) {\n",
"---\n",
"> \n",
"> static inline void compute_P_cold__Gamma_cold(CCTK_REAL rho_b,eos_struct &eos, CCTK_REAL &P_cold,CCTK_REAL &Gamma_cold) {\n",
"293c314,322\n",
"< if(rho_b==0.0) { P_cold = 0.0; gamma_cold = eos.gamma_tab[0]; return; }\n",
"---\n",
"> if(rho_b==0.0) { P_cold = 0.0; Gamma_cold = eos.Gamma_ppoly_tab[0]; return; }\n",
"> \n",
"> /***********************************\n",
"> * Piecewise Polytropic EOS Patch *\n",
"> * Computing P_cold and Gamma_cold *\n",
"> ***********************************/\n",
"> int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b);\n",
"> Gamma_cold = eos.Gamma_ppoly_tab[polytropic_index];\n",
"> P_cold = eos.K_ppoly_tab[polytropic_index]*pow(rho_b,Gamma_cold);\n",
"295,311d323\n",
"< if(eos.neos==1) {\n",
"< // Eq. 14 of http://arxiv.org/pdf/0802.0200.pdf :\n",
"< // P_{cold} = K_i rho_i^{\\Gamma_i}\n",
"< P_cold = eos.k_tab[0]*fasterpow_ppm_reconstruct(rho_b,eos.gamma_tab[0]);\n",
"< gamma_cold = eos.gamma_tab[0];\n",
"< return;\n",
"< }\n",
"< for(int nn=1;nn eos.rho_tab[nn-1]) {\n",
"< P_cold = eos.k_tab[nn]*fasterpow_ppm_reconstruct(rho_b,eos.gamma_tab[nn]);\n",
"< gamma_cold = eos.gamma_tab[nn];\n",
"< }\n",
"< }\n",
"< if (rho_b > eos.rho_tab[eos.neos-1]) {\n",
"< P_cold = eos.k_tab[eos.neos]*fasterpow_ppm_reconstruct(rho_b,eos.gamma_tab[eos.neos]);\n",
"< gamma_cold = eos.gamma_tab[eos.neos];\n",
"< }\n",
"313a326\n",
"> \n",
"334a348\n",
"> \n",
"338a353\n",
"> \n",
"345a361\n",
"> \n",
"348a365\n",
"> \n",
"351a369\n",
"> \n",
"359a378\n",
"> \n"
]
}
],
"source": [
"# Verify if the code generated by this tutorial module\n",
"# matches the original IllinoisGRMHD source code\n",
"\n",
"# First download the original IllinoisGRMHD source code\n",
"import urllib\n",
"from os import path\n",
"\n",
"original_IGM_file_url = \"https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/reconstruct_set_of_prims_PPM.C\"\n",
"original_IGM_file_name = \"reconstruct_set_of_prims_PPM-original.C\"\n",
"original_IGM_file_path = os.path.join(IGM_src_dir_path,original_IGM_file_name)\n",
"\n",
"# Then download the original IllinoisGRMHD source code\n",
"# We try it here in a couple of ways in an attempt to keep\n",
"# the code more portable\n",
"try:\n",
" original_IGM_file_code = urllib.request.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
" # Write down the file the original IllinoisGRMHD source code\n",
" with open(original_IGM_file_path,\"w\") as file:\n",
" file.write(original_IGM_file_code)\n",
"except:\n",
" try:\n",
" original_IGM_file_code = urllib.urlopen(original_IGM_file_url).read().decode(\"utf-8\")\n",
" # Write down the file the original IllinoisGRMHD source code\n",
" with open(original_IGM_file_path,\"w\") as file:\n",
" file.write(original_IGM_file_code)\n",
" except:\n",
" # If all else fails, hope wget does the job\n",
" !wget -O $original_IGM_file_path $original_IGM_file_url\n",
"\n",
"# Perform validation\n",
"Validation__reconstruct_set_of_prims_PPM__C = !diff $original_IGM_file_path $outfile_path__reconstruct_set_of_prims_PPM__C\n",
"\n",
"if Validation__reconstruct_set_of_prims_PPM__C == []:\n",
" # If the validation passes, we do not need to store the original IGM source code file\n",
" !rm $original_IGM_file_path\n",
" print(\"Validation test for reconstruct_set_of_prims_PPM.C: PASSED!\")\n",
"else:\n",
" # If the validation fails, we keep the original IGM source code file\n",
" print(\"Validation test for reconstruct_set_of_prims_PPM.C: FAILED!\")\n",
" # We also print out the difference between the code generated\n",
" # in this tutorial module and the original IGM source code\n",
" print(\"Diff:\")\n",
" for diff_line in Validation__reconstruct_set_of_prims_PPM__C:\n",
" print(diff_line)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"# Step 10: 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-IllinoisGRMHD_reconstruct_set_of_prims_PPM.pdf](Tutorial-IllinoisGRMHD_reconstruct_set_of_prims_PPM.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": 32,
"metadata": {},
"outputs": [],
"source": [
"latex_nrpy_style_path = os.path.join(nrpy_dir_path,\"latex_nrpy_style.tplx\")\n",
"#!jupyter nbconvert --to latex --template $latex_nrpy_style_path --log-level='WARN' Tutorial-IllinoisGRMHD__reconstruct_set_of_prims_PPM.ipynb\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD_reconstruct__set_of_prims_PPM.tex\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD_reconstruct__set_of_prims_PPM.tex\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD_reconstruct__set_of_prims_PPM.tex\n",
"!rm -f Tut*.out Tut*.aux Tut*.log"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}