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