{ "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", "## This notebook provides a comprehensive tutorial on the implementation of the Piecewise Parabolic Method (PPM) in the IllinoisGRMHD code, detailing functions for variable reconstruction, slope limitation, steepening, flattening, monotonizing, and index shifting. The tutorial incorporates the usage of relevant functions, methods to evaluate shock conditions, and code validation processes.\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. Flattening\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 within 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 the 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 need 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 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.1" } }, "nbformat": 4, "nbformat_minor": 4 }