{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Tutorial-IllinoisGRMHD: The conservative to primitive algorithm\n", "\n", "## Authors: Leo Werneck & Zach Etienne\n", "\n", "**This module is currently under development**\n", "\n", "## In this tutorial module we explain the algorithm used to get the primitive variables out of the conservative ones\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 2](#driver_conserv_to_prims): **`driver_conserv_to_prims.C`**\n", " 1. [Step 2.a](#adm_to_bssn__enforcing_detgammabar_equal_one): *Converting ADM quantities to BSSN quantities and enforcing $\\bar\\gamma=1$*\n", " 1. [Step 2.b](#equatorial_symmetry): *Applying equatorial symmetry*\n", " 1. [Step 2.c](#variable_setup): *Setting up the variables needed by `HARM`*\n", " 1. [Step 2.c.i](#variable_setup__bssn): BSSN quantities\n", " 1. [Step 2.c.ii](#variable_setup__prims): Primitives\n", " 1. [Step 2.c.iii](#variable_setup__conservs): Conservatives\n", " 1. [Step 2.c.iv](#variable_setup__lapse_and_psi): Lapse function and conformal factor\n", " 1. [Step 2.c.v](#variable_setup__phys_metric): Physical spatial metric\n", " 1. [Step 2.c.vi](#variable_setup__betadown_and_beta2): $\\beta_{i}$ and $\\beta^{2} \\equiv \\beta_{i}\\beta^{i}$\n", " 1. [Step 2.c.vii](#variable_setup__adm_4metric): The ADM 4-metric, $g_{\\mu\\nu}$, and its inverse, $g^{\\mu\\nu}$\n", " 1. [Step 2.c.viii](#variable_setup__temp_conservs): Temporary storage for current values of the conservative variables\n", " 1. [Step 2.d](#conserv_to_prim__driver): *Determining the primitives variables from the conservatives variables*\n", " 1. [Step 2.e](#enforce_limits_on_primitives_and_recompute_conservs): *Enforcing physical limits on primitives and recomputing the conservatives variables*\n", " 1. [Step 2.f](#updating_conservs_and_prims_gfs): *Updating conservative and primitive gridfunctions*\n", " 1. [Step 2.g](#diagnostics_and_debugging_tools): *Diagnostics and debugging tools*\n", "1. [Step 3](#harm_primitives_lowlevel): **`harm_primitives_lowlevel.C`**\n", " 1. [Step 3.a](#variables_needed_by_harm): *Setting up the variables needed by `HARM`*\n", " 1. [Step 3.a.i](#variables_needed_by_harm__detg): ${\\rm detg}$\n", " 1. [Step 3.a.ii](#variables_needed_by_harm__bi_harm): $B^{i}_{\\rm HARM}$\n", " 1. [Step 3.a.iii](#variables_needed_by_harm__init_rhob_pressure_vi): Initializing $\\rho_{b}$, $P$, and $v^{i}$\n", " 1. [Step 3.a.iv](#variables_needed_by_harm__original_conserv): Storing the original values of the conservative variables\n", " 1. [Step 3.a.v](#variables_needed_by_harm__guessing_rhob_pressure_vi): Guessing $\\rho_{b}$, $P$, and $v^{i}$\n", " 1. [Step 3.a.vi](#variables_needed_by_harm__conservs): Writing $\\boldsymbol{C}_{\\rm HARM}$ in terms of $\\boldsymbol{C}_{\\rm IGM}$\n", " 1. [Step 3.a.vii](#variables_needed_by_harm__prims): Writing $\\boldsymbol{P}_{\\rm HARM}$ in terms of $\\boldsymbol{P}_{\\rm IGM}$\n", " 1. [Step 3.b](#calling_harm_conservs_to_prims_solver): *Calling the `HARM` conservative-to-primitive solver*\n", " 1. [Step 3.c](#font_fix) *Applying the Font *et al.* fix, if the inversion fails*\n", " 1. [Step 3.d](#compute_utconi) *Compute $\\tilde{u}^{i}$*\n", " 1. [Step 3.e](#limiting_velocities) *Limiting velocities*\n", " 1. [Step 3.f](#primitives) *Setting the primitives*\n", "1. [Step 4](#font_fix_hybrid_eos): **`font_fix_hybrid_EOS.C`**\n", " 1. [Step 4.a](#font_fix_hybrid_eos__basic_quantities): Computing the basic quantities needed by the algorithm\n", " 1. [Step 4.b](#font_fix_hybrid_eos__sdots): Basic check: looking at $\\tilde{S}^{2}$\n", " 1. [Step 4.c](#font_fix_hybrid_eos__initial_guesses): Initial guesses for $W$, $S_{{\\rm fluid}}^{2}$, and $\\rho$\n", " 1. [Step 4.d](#font_fix_hybrid_eos__main_loop): The main loop\n", " 1. [Step 4.e](#font_fix_hybrid_eos__outputs): Output $\\rho_{b}$ and $u_{i}$\n", "1. [Step 5](#harm_primitives_headers): **`harm_primitives_headers.h`**\n", "1. [Step 6](#code_validation): **Code validation**\n", " 1. [Step 6.a](#driver_conserv_to_prims_validation): *`driver_conserv_to_prims.C`*\n", " 1. [Step 6.b](#harm_primitives_lowlevel_validation): *`harm_primitives_lowlevel.C`*\n", " 1. [Step 6.c](#font_fix_gamma_law_validation): *`font_fix_gamma_law.C`*\n", " 1. [Step 6.d](#harm_primitives_headers_validation): *`harm_primitives_headers.h`*\n", "1. [Step 7](#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__driver_conserv_to_prims__C = os.path.join(IGM_src_dir_path,\"driver_conserv_to_prims.C\")\n", "outfile_path__harm_primitives_lowlevel__C = os.path.join(IGM_src_dir_path,\"harm_primitives_lowlevel.C\")\n", "outfile_path__font_fix_hybrid_EOS__C = os.path.join(IGM_src_dir_path,\"font_fix_hybrid_EOS.C\")\n", "outfile_path__harm_primitives_headers__h = os.path.join(IGM_src_dir_path,\"harm_primitives_headers.h\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 1: Introduction \\[Back to [top](#toc)\\]\n", "$$\\label{introduction}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: `driver_conserv_to_prims.C` \\[Back to [top](#toc)\\]\n", "$$\\label{driver_conserv_to_prims}$$\n", "\n", "We start here by creating the `driver_conserv_to_prims.C` file and loading all files used by it. Note that of the files loaded, we have the following `IllinoisGRMHD` files:\n", "\n", "1. `harm_primitives_headers.h`: we will discuss this file [in step 4 of this tutorial module](#harm_primitives_headers).\n", "1. `inlined_functions.C`: this file is discussed in the [inlined_functions NRPy tutorial module](Tutorial-IllinoisGRMHD__inlined_functions.ipynb)\n", "1. `apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C`: this file is discussed in the [apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs NRPy tutorial module](Tutorial-IllinoisGRMHD__apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.ipynb)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile $outfile_path__driver_conserv_to_prims__C\n", "/* We evolve forward in time a set of functions called the\n", " * \"conservative variables\", and any time the conserv's\n", " * are updated, we must solve for the primitive variables\n", " * (rho, pressure, velocities) using a Newton-Raphson\n", " * technique, before reconstructing & evaluating the RHSs\n", " * of the MHD equations again.\n", " *\n", " * This file contains the driver routine for this Newton-\n", " * Raphson solver. Truncation errors in conservative\n", " * variables can lead to no physical solutions in\n", " * primitive variables. We correct for these errors here\n", " * through a number of tricks described in the appendices\n", " * of http://arxiv.org/pdf/1112.0568.pdf.\n", " *\n", " * This is a wrapper for the 2d solver of Noble et al. See\n", " * harm_utoprim_2d.c for references and copyright notice\n", " * for that solver. This wrapper was primarily written by\n", " * Zachariah Etienne & Yuk Tung Liu, in 2011-2013.\n", " *\n", " * For optimal compatibility, this wrapper is licensed under\n", " * the GPL v2 or any later version.\n", " *\n", " * Note that this code assumes a simple gamma law for the\n", " * moment, though it would be easy to extend to a piecewise\n", " * polytrope. */\n", "\n", "// Standard #include's\n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "\n", "\n", "#ifdef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", "#include \"standalone_conserv_to_prims_main_function.h\"\n", "#else\n", "#include \"cctk.h\"\n", "#include \"cctk_Arguments.h\"\n", "#include \"cctk_Parameters.h\"\n", "#include \"Symmetry.h\"\n", "\n", "#include \"IllinoisGRMHD_headers.h\"\n", "#include \"harm_primitives_headers.h\"\n", "#include \"harm_u2p_util.c\"\n", "#include \"inlined_functions.C\"\n", "#include \"apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.C\"\n", "\n", "extern \"C\" void IllinoisGRMHD_conserv_to_prims(CCTK_ARGUMENTS) {\n", " DECLARE_CCTK_ARGUMENTS;\n", " DECLARE_CCTK_PARAMETERS;\n", "\n", " // We use proper C++ here, for file I/O later.\n", " using namespace std;\n", "#endif\n", "\n", " /**********************************\n", " * Piecewise Polytropic EOS Patch *\n", " * Setting up the EOS struct *\n", " **********************************/\n", " /*\n", " * The short piece of code below takes care\n", " * of initializing the EOS parameters.\n", " * Please refer to the \"inlined_functions.C\"\n", " * source file for the documentation on the\n", " * function.\n", " */\n", " eos_struct eos;\n", " initialize_EOS_struct_from_input(eos);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.a: Converting ADM quantities to BSSN quantities and enforcing $\\bar\\gamma=1$ \\[Back to [top](#toc)\\]\n", "$$\\label{adm_to_bssn__enforcing_detgammabar_equal_one}$$\n", "\n", "Here we compute the conformal metric, $\\bar\\gamma_{ij}$, and its inverse, $\\bar\\gamma^{ij}$, from the physical metric $\\gamma_{ij}$. We also compute $\\phi$, the conformal factor, and $\\psi\\equiv e^{\\phi}$. Finally, we enforce the constraint $\\bar\\gamma = \\det\\left(\\bar\\gamma_{ij}\\right) = 1$. The entire procedure is explained in detail in the [convert ADM to BSSN and enforce determinant constraint NRPy tutorial module](Tutorial-IllinoisGRMHD__convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.ipynb)." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " // These BSSN-based variables are not evolved, and so are not defined anywhere that the grid has moved.\n", " // Here we convert ADM variables (from ADMBase) to the BSSN-based variables expected by this routine.\n", " IllinoisGRMHD_convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(cctkGH,cctk_lsh, gxx,gxy,gxz,gyy,gyz,gzz,alp,\n", " gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,\n", " gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,\n", " phi_bssn,psi_bssn,lapm1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.b: Applying equatorial symmetry \\[Back to [top](#toc)\\]\n", "$$\\label{equatorial_symmetry}$$\n", "\n", "We then use the [CardGrid3D ETK thorn](https://einsteintoolkit.org/thornguide/CactusBase/CartGrid3D/documentation.html) to apply equatorial symmetry to our problem." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", "#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", " if(CCTK_EQUALS(Symmetry,\"equatorial\")) {\n", " // SET SYMMETRY GHOSTZONES ON ALL CONSERVATIVE VARIABLES!\n", " int ierr=0;\n", " ierr+=CartSymGN(cctkGH,\"IllinoisGRMHD::grmhd_conservatives\");\n", " // FIXME: UGLY. Filling metric ghostzones is needed for, e.g., Cowling runs.\n", " ierr+=CartSymGN(cctkGH,\"lapse::lapse_vars\");\n", " ierr+=CartSymGN(cctkGH,\"bssn::BSSN_vars\");\n", " ierr+=CartSymGN(cctkGH,\"bssn::BSSN_AH\");\n", " ierr+=CartSymGN(cctkGH,\"shift::shift_vars\");\n", " if(ierr!=0) CCTK_VError(VERR_DEF_PARAMS,\"IllinoisGRMHD ERROR (grep for it, foo!) :(\");\n", " }\n", "#endif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.c: Setting up the variables needed by `HARM` \\[Back to [top](#toc)\\]\n", "$$\\label{variable_setup}$$\n", "\n", "We will now set up all the necessary variables to start the conservative to primitive algorithm. We begin by declaring useful debugging variables." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " //Start the timer, so we can benchmark the primitives solver during evolution.\n", " // Slower solver -> harder to find roots -> things may be going crazy!\n", " //FIXME: Replace this timing benchmark with something more meaningful, like the avg # of Newton-Raphson iterations per gridpoint!\n", " /*\n", " struct timeval start, end;\n", " long mtime, seconds, useconds;\n", " gettimeofday(&start, NULL);\n", " */\n", "\n", " int failures=0,font_fixes=0,vel_limited_ptcount=0;\n", " int pointcount=0;\n", " int failures_inhoriz=0;\n", " int pointcount_inhoriz=0;\n", "\n", " int pressure_cap_hit=0;\n", "\n", " CCTK_REAL error_int_numer=0,error_int_denom=0;\n", "\n", " int imin=0,jmin=0,kmin=0;\n", " int imax=cctk_lsh[0],jmax=cctk_lsh[1],kmax=cctk_lsh[2];\n", "\n", " int rho_star_fix_applied=0;\n", " long n_iter=0;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 2.c.i: BSSN quantities \\[Back to [top](#toc)\\]\n", "$$\\label{variable_setup__bssn}$$\n", "\n", "We start by loading the BSSN variables $\\left\\{\\phi, \\bar\\gamma_{ij}, \\alpha-1, \\beta^{i}, \\bar\\gamma^{ij}\\right\\}$ into a new array called $\\rm METRIC$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "#pragma omp parallel for reduction(+:failures,vel_limited_ptcount,font_fixes,pointcount,failures_inhoriz,pointcount_inhoriz,error_int_numer,error_int_denom,pressure_cap_hit,rho_star_fix_applied,n_iter) schedule(static)\n", " for(int k=kmin;k\n", "\n", "### Step 2.c.ii: Primitives \\[Back to [top](#toc)\\]\n", "$$\\label{variable_setup__prims}$$\n", "\n", "Then we load our current known values for the primitive variables $\\left\\{\\rho_{b}, P, v^{i}, B^{i}\\right\\}$ into a new array called $\\rm PRIMS$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " CCTK_REAL PRIMS[MAXNUMVARS];\n", " ww=0;\n", " PRIMS[ww] = rho_b[index]; ww++;\n", " PRIMS[ww] = P[index]; ww++;\n", " PRIMS[ww] = vx[index]; ww++;\n", " PRIMS[ww] = vy[index]; ww++;\n", " PRIMS[ww] = vz[index]; ww++;\n", " PRIMS[ww] = Bx[index]; ww++;\n", " PRIMS[ww] = By[index]; ww++;\n", " PRIMS[ww] = Bz[index]; ww++;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 2.c.iii: Conservatives \\[Back to [top](#toc)\\]\n", "$$\\label{variable_setup__conservs}$$\n", "\n", "Then we load our current known values for the conservative variables $\\left\\{\\rho_{\\star}, \\tilde{S}_{i}, \\tilde{\\tau}\\right\\}$ into a new array called $\\rm CONSERVS$." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " CCTK_REAL CONSERVS[NUM_CONSERVS] = {rho_star[index], mhd_st_x[index],mhd_st_y[index],mhd_st_z[index],tau[index]};" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 2.c.iv: Lapse function and conformal factor \\[Back to [top](#toc)\\]\n", "$$\\label{variable_setup__lapse_and_psi}$$\n", "\n", "Then we load the lapse function, $\\alpha$, and $\\psi$ related variables into the $\\rm METRIC\\_LAP\\_PSI4$ array. Notice that this is done using the ${\\rm SET\\_LAPSE\\_PSI4}()$ \"function\", which is defined in the `IllinoisGRMHD_headers.h` file from `IllinoisGRMHD`. The \"function\" itself is quite simple:\n", "\n", "```c\n", "#define SET_LAPSE_PSI4(array_name,METRIC) { \\\n", " array_name[LAPSE] = METRIC[LAPM1]+1.0; \\\n", " array_name[PSI2] = exp(2.0*METRIC[PHI]); \\\n", " array_name[PSI4] = SQR(array_name[PSI2]); \\\n", " array_name[PSI6] = array_name[PSI4]*array_name[PSI2]; \\\n", " array_name[PSIM4] = 1.0/array_name[PSI4]; \\\n", " array_name[LAPSEINV] = 1.0/array_name[LAPSE]; \\\n", " }\n", "```\n", "\n", "defining the quantities $\\left\\{\\alpha, \\psi^{2}, \\psi^{4}, \\psi^{6}, \\psi^{-4},\\alpha^{-1}\\right\\}$, where $\\psi=e^{\\phi}$." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " CCTK_REAL METRIC_LAP_PSI4[NUMVARS_METRIC_AUX];\n", " SET_LAPSE_PSI4(METRIC_LAP_PSI4,METRIC);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 2.c.v: Physical spatial metric \\[Back to [top](#toc)\\]\n", "$$\\label{variable_setup__phys_metric}$$\n", "\n", "Then we set up the physical spatial metric and its inverse through the relations\n", "\n", "$$\n", "\\boxed{\n", "\\begin{align}\n", "\\gamma_{ij} &= \\psi^{4} \\bar\\gamma_{ij}\\\\\n", "\\gamma^{ij} &= \\psi^{-4}\\bar\\gamma^{ij}\n", "\\end{align}\n", "}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " CCTK_REAL METRIC_PHYS[NUMVARS_FOR_METRIC];\n", " METRIC_PHYS[GXX] = METRIC[GXX]*METRIC_LAP_PSI4[PSI4];\n", " METRIC_PHYS[GXY] = METRIC[GXY]*METRIC_LAP_PSI4[PSI4];\n", " METRIC_PHYS[GXZ] = METRIC[GXZ]*METRIC_LAP_PSI4[PSI4];\n", " METRIC_PHYS[GYY] = METRIC[GYY]*METRIC_LAP_PSI4[PSI4];\n", " METRIC_PHYS[GYZ] = METRIC[GYZ]*METRIC_LAP_PSI4[PSI4];\n", " METRIC_PHYS[GZZ] = METRIC[GZZ]*METRIC_LAP_PSI4[PSI4];\n", " METRIC_PHYS[GUPXX] = METRIC[GUPXX]*METRIC_LAP_PSI4[PSIM4];\n", " METRIC_PHYS[GUPXY] = METRIC[GUPXY]*METRIC_LAP_PSI4[PSIM4];\n", " METRIC_PHYS[GUPXZ] = METRIC[GUPXZ]*METRIC_LAP_PSI4[PSIM4];\n", " METRIC_PHYS[GUPYY] = METRIC[GUPYY]*METRIC_LAP_PSI4[PSIM4];\n", " METRIC_PHYS[GUPYZ] = METRIC[GUPYZ]*METRIC_LAP_PSI4[PSIM4];\n", " METRIC_PHYS[GUPZZ] = METRIC[GUPZZ]*METRIC_LAP_PSI4[PSIM4];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 2.c.vi: $\\beta_{i}$ and $\\beta^{2} \\equiv \\beta_{i}\\beta^{i}$ \\[Back to [top](#toc)\\]\n", "$$\\label{variable_setup__betadown_and_beta2}$$\n", "\n", "We then evaluate\n", "\n", "$$\n", "\\beta_{i} = \\gamma_{ij}\\beta^{j} \\implies\n", "\\boxed{\n", "\\left\\{\n", "\\begin{align}\n", "\\beta_{x} &= \\gamma_{xx}\\beta^{x} + \\gamma_{xy}\\beta^{y} + \\gamma_{xz}\\beta^{z}\\\\\n", "\\beta_{y} &= \\gamma_{yx}\\beta^{x} + \\gamma_{yy}\\beta^{y} + \\gamma_{yz}\\beta^{z}\\\\\n", "\\beta_{z} &= \\gamma_{zx}\\beta^{x} + \\gamma_{zy}\\beta^{y} + \\gamma_{zz}\\beta^{z}\n", "\\end{align}\n", "\\right.\n", "}\\ ,\n", "$$\n", "\n", "and\n", "\n", "$$\n", "\\boxed{\\beta^{2} \\equiv \\beta_{i}\\beta^{i} = \\beta_{x}\\beta^{x} + \\beta_{y}\\beta^{y} + \\beta_{z}\\beta^{z}}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " CCTK_REAL TUPMUNU[10],TDNMUNU[10];\n", "\n", " CCTK_REAL shift_xL = METRIC_PHYS[GXX]*METRIC[SHIFTX] + METRIC_PHYS[GXY]*METRIC[SHIFTY] + METRIC_PHYS[GXZ]*METRIC[SHIFTZ];\n", " CCTK_REAL shift_yL = METRIC_PHYS[GXY]*METRIC[SHIFTX] + METRIC_PHYS[GYY]*METRIC[SHIFTY] + METRIC_PHYS[GYZ]*METRIC[SHIFTZ];\n", " CCTK_REAL shift_zL = METRIC_PHYS[GXZ]*METRIC[SHIFTX] + METRIC_PHYS[GYZ]*METRIC[SHIFTY] + METRIC_PHYS[GZZ]*METRIC[SHIFTZ];\n", " CCTK_REAL beta2L = shift_xL*METRIC[SHIFTX] + shift_yL*METRIC[SHIFTY] + shift_zL*METRIC[SHIFTZ];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 2.c.vii: The ADM 4-metric, $g_{\\mu\\nu}$, and its inverse, $g^{\\mu\\nu}$ \\[Back to [top](#toc)\\]\n", "$$\\label{variable_setup__adm_4metric}$$\n", "\n", "We then setup the ADM 4-metric and its inverse. We refer the reader to eqs. (2.119) and (2.122) from [Baumgarte & Shapiro's Numerical Relativity (2010)](https://books.google.com/books/about/Numerical_Relativity.html?id=dxU1OEinvRUC), which are repeated here for the sake of the reader (in reverse order)\n", "\n", "$$\n", "\\boxed{g_{\\mu\\nu} =\n", "\\begin{pmatrix}\n", "-\\alpha^{2} + \\beta_{\\ell}\\beta^{\\ell} & \\beta_{i}\\\\\n", "\\beta_{j} & \\gamma_{ij}\n", "\\end{pmatrix}}\\ .\n", "$$\n", "\n", "and\n", "\n", "$$\n", "\\boxed{\n", "g^{\\mu\\nu} =\n", "\\begin{pmatrix}\n", "-\\alpha^{-2} & \\alpha^{-2}\\beta^{i}\\\\\n", "\\alpha^{-2}\\beta^{j} & \\gamma^{ij} - \\alpha^{-2}\\beta^{i}\\beta^{j}\n", "\\end{pmatrix}\n", "}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " // Compute 4-metric, both g_{\\mu \\nu} and g^{\\mu \\nu}.\n", " // This is for computing T_{\\mu \\nu} and T^{\\mu \\nu}. Also the HARM con2prim lowlevel function requires them.\n", " CCTK_REAL g4dn[4][4],g4up[4][4];\n", " g4dn[0][0] = -SQR(METRIC_LAP_PSI4[LAPSE]) + beta2L;\n", " g4dn[0][1] = g4dn[1][0] = shift_xL;\n", " g4dn[0][2] = g4dn[2][0] = shift_yL;\n", " g4dn[0][3] = g4dn[3][0] = shift_zL;\n", " g4dn[1][1] = METRIC_PHYS[GXX];\n", " g4dn[1][2] = g4dn[2][1] = METRIC_PHYS[GXY];\n", " g4dn[1][3] = g4dn[3][1] = METRIC_PHYS[GXZ];\n", " g4dn[2][2] = METRIC_PHYS[GYY];\n", " g4dn[2][3] = g4dn[3][2] = METRIC_PHYS[GYZ];\n", " g4dn[3][3] = METRIC_PHYS[GZZ];\n", "\n", " CCTK_REAL alpha_inv_squared=SQR(METRIC_LAP_PSI4[LAPSEINV]);\n", " g4up[0][0] = -1.0*alpha_inv_squared;\n", " g4up[0][1] = g4up[1][0] = METRIC[SHIFTX]*alpha_inv_squared;\n", " g4up[0][2] = g4up[2][0] = METRIC[SHIFTY]*alpha_inv_squared;\n", " g4up[0][3] = g4up[3][0] = METRIC[SHIFTZ]*alpha_inv_squared;\n", " g4up[1][1] = METRIC_PHYS[GUPXX] - METRIC[SHIFTX]*METRIC[SHIFTX]*alpha_inv_squared;\n", " g4up[1][2] = g4up[2][1] = METRIC_PHYS[GUPXY] - METRIC[SHIFTX]*METRIC[SHIFTY]*alpha_inv_squared;\n", " g4up[1][3] = g4up[3][1] = METRIC_PHYS[GUPXZ] - METRIC[SHIFTX]*METRIC[SHIFTZ]*alpha_inv_squared;\n", " g4up[2][2] = METRIC_PHYS[GUPYY] - METRIC[SHIFTY]*METRIC[SHIFTY]*alpha_inv_squared;\n", " g4up[2][3] = g4up[3][2] = METRIC_PHYS[GUPYZ] - METRIC[SHIFTY]*METRIC[SHIFTZ]*alpha_inv_squared;\n", " g4up[3][3] = METRIC_PHYS[GUPZZ] - METRIC[SHIFTZ]*METRIC[SHIFTZ]*alpha_inv_squared;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 2.c.viii: Temporary storage for current values of the conservative variables \\[Back to [top](#toc)\\]\n", "$$\\label{variable_setup__temp_conservs}$$\n", "\n", "Instead of declaring new variables to store the currently known values of the conservative variables $\\left\\{\\rho_{\\star}, \\tilde{S}_{i}, \\tilde{\\tau}\\right\\}$, we will simply use the flux variables $\\left\\{\\rho_{\\star}^{\\rm flux}, \\tilde{S}_{i}^{\\rm flux}, \\tilde{\\tau}^{\\rm flux}\\right\\}$ which are used by `IllinoisGRMHD` as temporary storage. This is done for debugging purposes.\n", "\n", "We also store the original values in the variables $\\left\\{\\rho_{\\star}^{\\rm orig}, \\tilde{S}_{i}^{\\rm orig}, \\tilde{\\tau}^{\\rm orig}\\right\\}$." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " //FIXME: might slow down the code.\n", " if(isnan(CONSERVS[RHOSTAR]*CONSERVS[STILDEX]*CONSERVS[STILDEY]*CONSERVS[STILDEZ]*CONSERVS[TAUENERGY]*PRIMS[BX_CENTER]*PRIMS[BY_CENTER]*PRIMS[BZ_CENTER])) {\n", " CCTK_VInfo(CCTK_THORNSTRING,\"NAN FOUND: i,j,k = %d %d %d, x,y,z = %e %e %e , index=%d st_i = %e %e %e, rhostar = %e, tau = %e, Bi = %e %e %e, gij = %e %e %e %e %e %e, Psi6 = %e\",\n", " i,j,k,x[index],y[index],z[index],index,\n", " CONSERVS[STILDEX],CONSERVS[STILDEY],CONSERVS[STILDEZ],CONSERVS[RHOSTAR],CONSERVS[TAUENERGY],\n", " PRIMS[BX_CENTER],PRIMS[BY_CENTER],PRIMS[BZ_CENTER],METRIC_PHYS[GXX],METRIC_PHYS[GXY],METRIC_PHYS[GXZ],METRIC_PHYS[GYY],METRIC_PHYS[GYZ],METRIC_PHYS[GZZ],METRIC_LAP_PSI4[PSI6]);\n", " }\n", "\n", " // Here we use _flux variables as temp storage for original values of conservative variables.. This is used for debugging purposes only.\n", " rho_star_flux[index] = CONSERVS[RHOSTAR];\n", " st_x_flux[index] = CONSERVS[STILDEX];\n", " st_y_flux[index] = CONSERVS[STILDEY];\n", " st_z_flux[index] = CONSERVS[STILDEZ];\n", " tau_flux[index] = CONSERVS[TAUENERGY];\n", "\n", " CCTK_REAL rho_star_orig = CONSERVS[RHOSTAR];\n", " CCTK_REAL mhd_st_x_orig = CONSERVS[STILDEX];\n", " CCTK_REAL mhd_st_y_orig = CONSERVS[STILDEY];\n", " CCTK_REAL mhd_st_z_orig = CONSERVS[STILDEZ];\n", " CCTK_REAL tau_orig = CONSERVS[TAUENERGY];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.d: Determining the primitives variables from the conservatives variables \\[Back to [top](#toc)\\]\n", "$$\\label{conserv_to_prim__driver}$$\n", "\n", "In this part of the code we determine the primitives variables from the conservatives variables. Please note that this is only the driver function. The algorithm is discussed on [Step 3](#harm_primitives_lowlevel) below.\n", "\n", "We start by calling the `apply_tau_floor()` function to ensure that the value of $\\tilde\\tau$ is not out of physical range. This function is discussed in the [apply $\\tilde\\tau$ floor and enforce limits on primitives NRPy+ tutorial module](Tutorial-IllinoisGRMHD__apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.ipynb).\n", "\n", "Notice that this is only performed when $\\rho_{\\star}>0$. If this is not the case, we fix the issue by setting $\\rho_{b} = \\rho_{b}^{\\rm atm}$ and the conservative to primitive algorithm is skipped altogether.\n", "\n", "When $\\rho_{\\star}>0$, we call upon the primitive to conservatives algorithm through the `harm_primitives_gammalaw_lowlevel()` function, described [below](#fixme)." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " int check=0;\n", " struct output_stats stats;\n", " stats.n_iter=0;\n", " stats.vel_limited=0;\n", " stats.failure_checker=0;\n", "\n", " if(CONSERVS[RHOSTAR]>0.0) {\n", " // Apply the tau floor\n", " apply_tau_floor(index,tau_atm,rho_b_atm,Psi6threshold,PRIMS,METRIC,METRIC_PHYS,METRIC_LAP_PSI4,stats,eos, CONSERVS);\n", "\n", " stats.font_fixed=0;\n", " for(int ii=0;ii<3;ii++) {\n", " check = harm_primitives_gammalaw_lowlevel(index,i,j,k,x,y,z,METRIC,METRIC_PHYS,METRIC_LAP_PSI4,\n", " CONSERVS,PRIMS, g4dn,g4up, stats,eos);\n", " if(check==0) ii=4;\n", " else stats.failure_checker+=100000;\n", " }\n", " } else {\n", " stats.failure_checker+=1;\n", " // Set to atmosphere if rho_star<0.\n", " //FIXME: FOR GAMMA=2 ONLY:\n", " PRIMS[RHOB] = rho_b_atm;\n", "\n", " /* Set P = P_cold */\n", " int polytropic_index = find_polytropic_K_and_Gamma_index(eos, rho_b_atm);\n", " CCTK_REAL K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];\n", " CCTK_REAL Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];\n", " PRIMS[PRESSURE] = K_ppoly_tab*pow(rho_b_atm,Gamma_ppoly_tab);\n", "\n", " PRIMS[VX] =-METRIC[SHIFTX];\n", " PRIMS[VY] =-METRIC[SHIFTY];\n", " PRIMS[VZ] =-METRIC[SHIFTZ];\n", "\n", " rho_star_fix_applied++;\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.e: Enforcing physical limits on primitives and recomputing the conservatives variables \\[Back to [top](#toc)\\]\n", "$$\\label{enforce_limits_on_primitives_and_recompute_conservs}$$\n", "\n", "We now enforce that the values of the primitives variables are physically meaningful. This is done by calling the `IllinoisGRMHD_enforce_limits_on_primitives_and_recompute_conservs()` function, which is described in the [enforce physical limits on primitives and recomputive conservatives NRPy+ tutorial module](Tutorial-IllinoisGRMHD__apply_tau_floor__enforce_limits_on_primitives_and_recompute_conservs.ipynb)." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " // Enforce limits on primitive variables and recompute conservatives.\n", " static const int already_computed_physical_metric_and_inverse=1;\n", " IllinoisGRMHD_enforce_limits_on_primitives_and_recompute_conservs(already_computed_physical_metric_and_inverse,PRIMS,stats,eos,METRIC,g4dn,g4up, TUPMUNU,TDNMUNU,CONSERVS);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.f: Updating conservative and primitive gridfunctions \\[Back to [top](#toc)\\]\n", "$$\\label{updating_conservs_and_prims_gfs}$$\n", "\n", "Then we update the corresponding conservative and primitive gridfunctions with the updated values we just computed." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " rho_star[index] = CONSERVS[RHOSTAR];\n", " mhd_st_x[index] = CONSERVS[STILDEX];\n", " mhd_st_y[index] = CONSERVS[STILDEY];\n", " mhd_st_z[index] = CONSERVS[STILDEZ];\n", " tau[index] = CONSERVS[TAUENERGY];\n", "\n", " // Set primitives, and/or provide a better guess.\n", " rho_b[index] = PRIMS[RHOB];\n", " P[index] = PRIMS[PRESSURE];\n", " vx[index] = PRIMS[VX];\n", " vy[index] = PRIMS[VY];\n", " vz[index] = PRIMS[VZ];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.g: Diagnostics and debugging tools \\[Back to [top](#toc)\\]\n", "$$\\label{diagnostics_and_debugging_tools}$$\n", "\n", "Now we simply append to the file useful diagnostics and debugging tools for our code." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/driver_conserv_to_prims.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__driver_conserv_to_prims__C\n", "\n", "\n", " if(update_Tmunu) {\n", " ww=0;\n", " eTtt[index] = TDNMUNU[ww]; ww++;\n", " eTtx[index] = TDNMUNU[ww]; ww++;\n", " eTty[index] = TDNMUNU[ww]; ww++;\n", " eTtz[index] = TDNMUNU[ww]; ww++;\n", " eTxx[index] = TDNMUNU[ww]; ww++;\n", " eTxy[index] = TDNMUNU[ww]; ww++;\n", " eTxz[index] = TDNMUNU[ww]; ww++;\n", " eTyy[index] = TDNMUNU[ww]; ww++;\n", " eTyz[index] = TDNMUNU[ww]; ww++;\n", " eTzz[index] = TDNMUNU[ww];\n", " }\n", "\n", " //Finally, we set h, the enthalpy:\n", " //CCTK_REAL eps = P[index]/rho_b[index]/(GAMMA-1.0);\n", " //h[index] = 1.0 + P[index]/rho_b[index] + eps;\n", "\n", " /***************************************************************************************************************************/\n", " // DIAGNOSTICS:\n", " //Pressure cap hit?\n", " /* FIXME\n", " CCTK_REAL P_cold = rho_b[index]*rho_b[index];\n", " if(P[index]/P_cold > 0.99*1e3 && rho_b[index]>100.0*rho_b_atm) {\n", " if(exp(phi[index]*6.0) <= Psi6threshold) pressure_cap_hit++;\n", " }\n", " */\n", "\n", " //Now we compute the difference between original & new conservatives, for diagnostic purposes:\n", " error_int_numer += fabs(tau[index] - tau_orig) + fabs(rho_star[index] - rho_star_orig) +\n", " fabs(mhd_st_x[index] - mhd_st_x_orig) + fabs(mhd_st_y[index] - mhd_st_y_orig) + fabs(mhd_st_z[index] - mhd_st_z_orig);\n", " error_int_denom += tau_orig + rho_star_orig + fabs(mhd_st_x_orig) + fabs(mhd_st_y_orig) + fabs(mhd_st_z_orig);\n", "\n", " if(stats.font_fixed==1) font_fixes++;\n", " vel_limited_ptcount+=stats.vel_limited;\n", " if(check!=0) {\n", " failures++;\n", " if(exp(METRIC[PHI]*6.0)>Psi6threshold) {\n", " failures_inhoriz++;\n", " pointcount_inhoriz++;\n", " }\n", " }\n", " pointcount++;\n", " /***************************************************************************************************************************/\n", " failure_checker[index] = stats.failure_checker;\n", " n_iter += stats.n_iter;\n", " }\n", "\n", " /*\n", " gettimeofday(&end, NULL);\n", "\n", " seconds = end.tv_sec - start.tv_sec;\n", " useconds = end.tv_usec - start.tv_usec;\n", "\n", " mtime = ((seconds) * 1000 + useconds/1000.0) + 0.999; // We add 0.999 since mtime is a long int; this rounds up the result before setting the value. Here, rounding down is incorrect.\n", " solutions per second: cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2] / ((CCTK_REAL)mtime/1000.0),\n", " */\n", " if(CCTK_Equals(verbose, \"essential\") || CCTK_Equals(verbose, \"essential+iteration output\")) {\n", " CCTK_VInfo(CCTK_THORNSTRING,\"C2P: Lev: %d NumPts= %d | Fixes: Font= %d VL= %d rho*= %d | Failures: %d InHoriz= %d / %d | Error: %.3e, ErrDenom: %.3e | %.2f iters/gridpt\",\n", " (int)GetRefinementLevel(cctkGH),\n", " pointcount,font_fixes,vel_limited_ptcount,rho_star_fix_applied,\n", " failures,\n", " failures_inhoriz,pointcount_inhoriz,\n", " error_int_numer/error_int_denom,error_int_denom,\n", " (double)n_iter/( (double)(cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2]) ));\n", " }\n", "\n", " if(pressure_cap_hit!=0) {\n", " //CCTK_VInfo(CCTK_THORNSTRING,\"PRESSURE CAP HIT %d TIMES! Outputting debug file!\",pressure_cap_hit);\n", " }\n", "\n", " // Very useful con2prim debugger. If the primitives (con2prim) solver fails, this will output all data needed to\n", " // debug where and why the solver failed. Strongly suggested for experimenting with new fixes.\n", " if(conserv_to_prims_debug==1 && error_int_numer/error_int_denom > 0.05) {\n", "\n", " ofstream myfile;\n", " char filename[100];\n", " srand(time(NULL));\n", " sprintf(filename,\"primitives_debug-%e.dat\",error_int_numer/error_int_denom);\n", " //Alternative, for debugging purposes as well:\n", " //srand(time(NULL));\n", " //sprintf(filename,\"primitives_debug-%d.dat\",rand());\n", " myfile.open (filename, ios::out | ios::binary);\n", " //myfile.open (\"data.bin\", ios::out | ios::binary);\n", " myfile.write((char*)cctk_lsh, 3*sizeof(int));\n", "\n", " myfile.write((char*)&GAMMA_SPEED_LIMIT, 1*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char*)&rho_b_max, 1*sizeof(CCTK_REAL));\n", " myfile.write((char*)&rho_b_atm, 1*sizeof(CCTK_REAL));\n", " myfile.write((char*)&tau_atm, 1*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char*)&Psi6threshold, 1*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char*)&update_Tmunu, 1*sizeof(int));\n", "\n", " myfile.write((char*)&neos, 1*sizeof(int));\n", " myfile.write((char*)&Gamma_th, 1*sizeof(CCTK_REAL));\n", " myfile.write((char*)&K_ppoly_tab0, 1*sizeof(CCTK_REAL));\n", " myfile.write((char*)Gamma_ppoly_tab_in, neos*sizeof(CCTK_REAL));\n", " myfile.write((char*)rho_ppoly_tab_in, (neos-1)*sizeof(CCTK_REAL));\n", "\n", " int fullsize=cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2];\n", " myfile.write((char*)x, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)y, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)z, (fullsize)*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char *)failure_checker, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)eTtt, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)eTtx, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)eTty, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)eTtz, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)eTxx, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)eTxy, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)eTxz, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)eTyy, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)eTyz, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)eTzz, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)alp, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)gxx, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)gxy, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)gxz, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)gyy, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)gyz, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)gzz, fullsize*sizeof(CCTK_REAL));\n", " myfile.write((char *)psi_bssn, fullsize*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char*)phi_bssn, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtxx, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtxy, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtxz, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtyy, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtyz, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtzz, (fullsize)*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char*)gtupxx, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtupxy, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtupxz, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtupyy, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtupyz, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)gtupzz, (fullsize)*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char*)betax, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)betay, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)betaz, (fullsize)*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char*)lapm1, (fullsize)*sizeof(CCTK_REAL));\n", "\n", " // HERE WE USE _flux variables as temp storage for original values of conservative variables.. This is used for debugging purposes only.\n", " myfile.write((char*)tau_flux, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)st_x_flux, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)st_y_flux, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)st_z_flux, (fullsize)*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char*)rho_star_flux, (fullsize)*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char*)Bx, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)By, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)Bz, (fullsize)*sizeof(CCTK_REAL));\n", "\n", " myfile.write((char*)vx, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)vy, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)vz, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)P, (fullsize)*sizeof(CCTK_REAL));\n", " myfile.write((char*)rho_b,(fullsize)*sizeof(CCTK_REAL));\n", "\n", " int checker=1063; myfile.write((char*)&checker,sizeof(int));\n", "\n", " myfile.close();\n", " CCTK_VInfo(CCTK_THORNSTRING,\"Finished writing %s\",filename);\n", " }\n", "\n", "#ifdef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", " return 0; // int main() requires an integer be returned\n", "#endif\n", "\n", "}\n", "\n", "#include \"harm_primitives_lowlevel.C\"\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 3: `harm_primitives_lowlevel.C` \\[Back to [top](#toc)\\]\n", "$$\\label{harm_primitives_lowlevel}$$\n", "\n", "We will now begin documenting the `harm_primitives_lowlevel.C` code. We will start by declaring the `harm_primitives_gammalaw_lowlevel()` function and loading up the EOS parameters:\n", "\n", "$$\n", "\\begin{align}\n", "{\\rm kpoly} &:= \\kappa\\ ,\\\\\n", "{\\rm gamma} &:= \\Gamma\\ ,\n", "\\end{align}\n", "$$\n", "\n", "where we are currently assuming the single polytrope EOS\n", "\n", "$$\n", "P = \\kappa \\rho_{b}^{\\Gamma}\\ .\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.a: Setting up the variables needed by `HARM` \\[Back to [top](#toc)\\]\n", "$$\\label{variables_needed_by_harm}$$\n", "\n", "In this section we will describe the final set of manipulations that are required *before* calling the conservative to primitive algorithm. These manipulations are necessary because we make use of the [`HARM` software](https://arxiv.org/abs/astro-ph/0301509), which uses a different set of conservative/primitives variables than `IllinoisGRMHD`.\n", "\n", "Notice that `IllinoisGRMHD` uses the set of conservative variables $\\boldsymbol{C}_{\\rm IGM} = \\left\\{\\rho_{\\star},\\tilde{\\tau},\\tilde{S}_{i},\\tilde{B}^{i}\\right\\}$ and the primitive variables $\\boldsymbol{P}_{\\rm IGM} = \\left\\{\\rho_{b},P,v^{i},B^{i}\\right\\}$, while `HARM` makes use of the conservative variables $\\boldsymbol{C}_{\\rm HARM} = \\left\\{\\sqrt{-g}\\rho_{b}u^{0},\\sqrt{-g}\\left(T^{0}_{\\ 0}+\\rho_{b}u^{0}\\right),\\sqrt{-g}T^{0}_{\\ i},\\sqrt{-g}B^{i}_{\\rm HARM}\\right\\}$ and primitive variables $\\boldsymbol{P}_{\\rm HARM} = \\left\\{\\rho_{b},u,\\tilde{u}^{i},B^{i}_{\\rm HARM}\\right\\}$. The key step here is then writing `HARM`'s variables in terms of `IllinoisGRMHD`'s variables.\n", "\n", "\n", "\n", "### Step 3.a.i: ${\\rm detg}$ \\[Back to [top](#toc)\\]\n", "$$\\label{variables_needed_by_harm__detg}$$\n", "\n", "The first variable we will need to compute is ${\\rm detg} := \\sqrt{-g}$, where $g=\\det(g_{\\mu\\nu})$, with $g_{\\mu\\nu}$ the [physical ADM 4-metric evaluated above](#variable_setup__adm_4metric). To relate $\\sqrt{-g}$ with `IllinoisGRMHD`'s variables, we will make use of the well known relation (see eq. 2.124 in [Baumgarte & Shapiro's Numerical Relativity](https://books.google.com/books/about/Numerical_Relativity.html?id=dxU1OEinvRUC))\n", "\n", "$$\n", "\\boxed{\\sqrt{-g} = \\alpha\\sqrt{\\gamma} = \\alpha\\psi^{6}}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile $outfile_path__harm_primitives_lowlevel__C\n", "inline int harm_primitives_gammalaw_lowlevel(const int index,const int i,const int j,const int k,CCTK_REAL *X,CCTK_REAL *Y,CCTK_REAL *Z,\n", " CCTK_REAL *METRIC,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,\n", " CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,\n", " CCTK_REAL g4dn[NDIM][NDIM],CCTK_REAL g4up[NDIM][NDIM],\n", " struct output_stats &stats, eos_struct &eos) {\n", "#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", " DECLARE_CCTK_PARAMETERS;\n", "#endif" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " // declare some variables for HARM.\n", " CCTK_REAL U[NPR];\n", " CCTK_REAL prim[NPR];\n", " CCTK_REAL detg = METRIC_LAP_PSI4[LAPSE]*METRIC_LAP_PSI4[PSI6]; // == alpha sqrt{gamma} = alpha Psi^6\n", "\n", " // Check to see if the metric is positive-definite.\n", " // Note that this will slow down the code, and if the metric doesn't obey this, the run is probably too far gone to save,\n", " // though if it happens deep in the horizon, it might resurrect the run.\n", " /*\n", " CCTK_REAL lam1,lam2,lam3;\n", " CCTK_REAL M11 = METRIC[GXX], M12=METRIC[GXY], M13=METRIC[GXZ], M22=METRIC[GYY], M23=METRIC[GYZ], M33=METRIC[GZZ];\n", " eigenvalues_3by3_real_sym_matrix(lam1, lam2, lam3,M11, M12, M13, M22, M23, M33);\n", " if (lam1 < 0.0 || lam2 < 0.0 || lam3 < 0.0) {\n", " // Metric is not positive-defitive, reset the physical metric to be conformally-flat.\n", " METRIC_PHYS[GXX] = METRIC_LAP_PSI4[PSI4];\n", " METRIC_PHYS[GXY] = 0.0;\n", " METRIC_PHYS[GXZ] = 0.0;\n", " METRIC_PHYS[GYY] = METRIC_LAP_PSI4[PSI4];\n", " METRIC_PHYS[GYZ] = 0.0;\n", " METRIC_PHYS[GZZ] = METRIC_LAP_PSI4[PSI4];\n", " METRIC_PHYS[GUPXX] = METRIC_LAP_PSI4[PSIM4];\n", " METRIC_PHYS[GUPXY] = 0.0;\n", " METRIC_PHYS[GUPXZ] = 0.0;\n", " METRIC_PHYS[GUPYY] = METRIC_LAP_PSI4[PSIM4];\n", " METRIC_PHYS[GUPYZ] = 0.0;\n", " METRIC_PHYS[GUPZZ] = METRIC_LAP_PSI4[PSIM4];\n", " }\n", " */" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.a.ii: $B^{i}_{\\rm HARM}$ \\[Back to [top](#toc)\\]\n", "$$\\label{variables_needed_by_harm__bi_harm}$$\n", "\n", "Now we must relate `HARM`'s $B^{i}_{\\rm HARM}$ with the variables used by `IllinoisGRMHD`. We can start by looking at eqs. (23), (24), and (31) in [Duez *et al.* (2005)](https://arxiv.org/pdf/astro-ph/0503420.pdf) to find the following useful relations\n", "\n", "$$\n", "{\\rm IllinoisGRMHD:}\\ \n", "\\left\\{\n", "\\begin{align}\n", "\\sqrt{4\\pi}b^{0} &= \\frac{u_{i}B^{i}}{\\alpha}\\ ,\\\\\n", "\\sqrt{4\\pi}b^{i} &= \\frac{B^{i}/\\alpha + \\sqrt{4\\pi}b^{0}u^{i}}{u^{0}}\\ .\n", "\\end{align}\n", "\\right.\n", "$$\n", "\n", "Now, if we look at eqs. (16) and (17) in [Gammie & McKinney (2003)](https://arxiv.org/pdf/astro-ph/0301509.pdf) we find the following relations\n", "\n", "$$\n", "{\\rm HARM:}\\ \n", "\\left\\{\n", "\\begin{align}\n", "b^{0} &= u_{i}B^{i}_{\\rm HARM}\\ ,\\\\\n", "b^{i} &= \\frac{B^{i}_{\\rm HARM} + b^{0}u^{i}}{u^{0}}\\ ,\n", "\\end{align}\n", "\\right.\n", "$$\n", "\n", "from which we can then find the relation\n", "\n", "$$\n", "\\boxed{B^{i}_{\\rm HARM} = \\frac{B^{i}}{\\alpha\\sqrt{4\\pi}}}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " // Note that ONE_OVER_SQRT_4PI gets us to the object\n", " // referred to as B^i in the Noble et al paper (and\n", " // apparently also in the comments to their code).\n", " // This is NOT the \\mathcal{B}^i, which differs by\n", " // a factor of the lapse.\n", " CCTK_REAL BxL_over_alpha_sqrt_fourpi = PRIMS[BX_CENTER]*METRIC_LAP_PSI4[LAPSEINV]*ONE_OVER_SQRT_4PI;\n", " CCTK_REAL ByL_over_alpha_sqrt_fourpi = PRIMS[BY_CENTER]*METRIC_LAP_PSI4[LAPSEINV]*ONE_OVER_SQRT_4PI;\n", " CCTK_REAL BzL_over_alpha_sqrt_fourpi = PRIMS[BZ_CENTER]*METRIC_LAP_PSI4[LAPSEINV]*ONE_OVER_SQRT_4PI;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.a.iii: Initializing $\\rho_{b}$, $P$, and $v^{i}$ \\[Back to [top](#toc)\\]\n", "$$\\label{variables_needed_by_harm__init_rhob_pressure_vi}$$\n", "\n", "Here we simply initialize $\\rho_{b}$, $P$, and $v^{i}$." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " CCTK_REAL rho_b_oldL = PRIMS[RHOB];\n", " CCTK_REAL P_oldL = PRIMS[PRESSURE];\n", " CCTK_REAL vxL = PRIMS[VX];\n", " CCTK_REAL vyL = PRIMS[VY];\n", " CCTK_REAL vzL = PRIMS[VZ];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.a.iv: Storing the original values of the conservative variables \\[Back to [top](#toc)\\]\n", "$$\\label{variables_needed_by_harm__original_conserv}$$\n", "\n", "Here we store the currently known values of $\\left\\{\\rho_{\\star},\\tilde{\\tau},\\tilde{S}_{i}\\right\\}$." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " /*\n", " -- Driver for new prim. var. solver. The driver just translates\n", " between the two sets of definitions for U and P. The user may\n", " wish to alter the translation as they see fit.\n", "\n", "\n", " // / rho u^t \\ //\n", " // U = | T^t_t + rho u^t | * sqrt(-det(g_{\\mu\\nu})) //\n", " // | T^t_i | //\n", " // \\ B^i / //\n", " // //\n", " // / rho \\ //\n", " // P = | uu | //\n", " // | \\tilde{u}^i | //\n", " // \\ B^i / //\n", "\n", " (above equations have been fixed by Yuk Tung & Zach)\n", " */\n", "\n", " // U[NPR] = conserved variables (current values on input/output);\n", " // g4dn[NDIM][NDIM] = covariant form of the 4-metric ;\n", " // g4up[NDIM][NDIM] = contravariant form of the 4-metric ;\n", " // gdet = sqrt( - determinant of the 4-metric) ;\n", " // prim[NPR] = primitive variables (guess on input, calculated values on\n", " // output if there are no problems);\n", "\n", " // U[1] =\n", " // U[2-4] = stildei + rhostar\n", "\n", " CCTK_REAL rho_star_orig = CONSERVS[RHOSTAR];\n", " CCTK_REAL mhd_st_x_orig = CONSERVS[STILDEX];\n", " CCTK_REAL mhd_st_y_orig = CONSERVS[STILDEY];\n", " CCTK_REAL mhd_st_z_orig = CONSERVS[STILDEZ];\n", " CCTK_REAL tau_orig = CONSERVS[TAUENERGY];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.a.v: Guessing $\\rho_{b}$, $P$, and $v^{i}$ \\[Back to [top](#toc)\\]\n", "$$\\label{variables_needed_by_harm__guessing_rhob_pressure_vi}$$\n", "\n", "We will now start preparing the primitives for the conservative-to-primitive algorithm, which employs the [Newton-Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method).\n", "\n", "We offer the algorithm the initial guess\n", "\n", "$$\n", "\\boxed{\n", "{\\rm Guess\\ \\#1:}\\ \n", "\\left\\{\n", "\\begin{align}\n", "\\rho_{b}^{\\rm guess} &= \\frac{\\rho_{\\star}}{\\psi^{6}}\\\\\n", "P_{\\rm guess} &= \\kappa \\left(\\rho_{b}^{\\rm guess}\\right)^{\\Gamma}\\\\\n", "u0 &= \\frac{1}{\\alpha}\\\\\n", "v^{i} &= -\\beta^{i}\n", "\\end{align}\n", "\\right.\n", "}\\ .\n", "$$\n", "\n", "If this initial guess causes the Newton-Raphson method to fail, we try a second guess\n", "\n", "$$\n", "\\boxed{\n", "{\\rm Guess\\ \\#2:}\\ \n", "\\left\\{\n", "\\begin{align}\n", "\\rho_{b}^{\\rm guess} &= 100\\rho_{b}^{\\rm atm}\\\\\n", "P_{\\rm guess} &= \\kappa \\left(\\rho_{b}^{\\rm guess}\\right)^{\\Gamma}\\\\\n", "u0 &= \\frac{1}{\\alpha}\\\\\n", "v^{i} &= -\\beta^{i}\n", "\\end{align}\n", "\\right.\n", "}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " // Other ideas for setting the gamma speed limit\n", " //CCTK_REAL GAMMA_SPEED_LIMIT = 100.0;\n", " //if(METRIC_LAP_PSI4[PSI6]>Psi6threshold) GAMMA_SPEED_LIMIT=500.0;\n", " //if(METRIC_LAP_PSI4[PSI6]>Psi6threshold) GAMMA_SPEED_LIMIT=100.0;\n", "\n", " //FIXME: Only works if poisoning is turned on. Otherwise will access unknown memory. This trick alone speeds up the whole code (Cowling) by 2%.\n", " //int startguess=0;\n", " //if(std::isnan(PRIMS[VX])) startguess=1;\n", " int startguess=1;\n", "\n", " CCTK_REAL u0L=1.0;\n", " CCTK_REAL K_ppoly_tab,Gamma_ppoly_tab;\n", "\n", " for(int which_guess=startguess;which_guess<3;which_guess++) {\n", " int check;\n", "\n", " if(which_guess==1) {\n", " //Use a different initial guess:\n", " rho_b_oldL = CONSERVS[RHOSTAR]/METRIC_LAP_PSI4[PSI6];\n", "\n", " /**********************************\n", " * Piecewise Polytropic EOS Patch *\n", " * Finding Gamma_ppoly_tab and K_ppoly_tab *\n", " **********************************/\n", " /* Here we use our newly implemented\n", " * find_polytropic_K_and_Gamma() function\n", " * to determine the relevant polytropic\n", " * Gamma and K parameters to be used\n", " * within this function.\n", " */\n", " int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_oldL);\n", " K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];\n", " Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];\n", "\n", " // After that, we compute P_cold\n", " P_oldL = K_ppoly_tab*pow(rho_b_oldL,Gamma_ppoly_tab);\n", "\n", " u0L = METRIC_LAP_PSI4[LAPSEINV];\n", " vxL = -METRIC[SHIFTX];\n", " vyL = -METRIC[SHIFTY];\n", " vzL = -METRIC[SHIFTZ];\n", " }\n", "\n", " if(which_guess==2) {\n", " //Use atmosphere as initial guess:\n", " rho_b_oldL = 100.0*rho_b_atm;\n", "\n", " /**********************************\n", " * Piecewise Polytropic EOS Patch *\n", " * Finding Gamma_ppoly_tab and K_ppoly_tab *\n", " **********************************/\n", " /* Here we use our newly implemented\n", " * find_polytropic_K_and_Gamma() function\n", " * to determine the relevant polytropic\n", " * Gamma and K parameters to be used\n", " * within this function.\n", " */\n", " int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_oldL);\n", " K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];\n", " Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];\n", "\n", " // After that, we compute P_cold\n", " P_oldL = K_ppoly_tab*pow(rho_b_oldL,Gamma_ppoly_tab);\n", "\n", " u0L = METRIC_LAP_PSI4[LAPSEINV];\n", " vxL = -METRIC[SHIFTX];\n", " vyL = -METRIC[SHIFTY];\n", " vzL = -METRIC[SHIFTZ];\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.a.vi: Writing $\\boldsymbol{C}_{\\rm HARM}$ in terms of $\\boldsymbol{C}_{\\rm IGM}$ \\[Back to [top](#toc)\\]\n", "$$\\label{variables_needed_by_harm__conservs}$$\n", "\n", "We will now relate $\\boldsymbol{C}_{\\rm HARM}$ to $\\boldsymbol{C}_{\\rm IGM}$. As previously mentioned, we have\n", "\n", "$$\n", "\\begin{align}\n", "\\boldsymbol{C}_{\\rm IGM} &= \\left\\{\\rho_{\\star},\\tilde{\\tau},\\tilde{S}_{i},\\tilde{B}^{i}\\right\\}\\ ,\\\\\n", "\\boldsymbol{C}_{\\rm HARM} &= \\left\\{\\sqrt{-g}\\rho_{b}u^{0},\\sqrt{-g}\\left(T^{0}_{\\ 0}+\\rho_{b}u^{0}\\right),\\sqrt{-g}T^{0}_{\\ i},\\sqrt{-g}B^{i}_{\\rm HARM}\\right\\} \\equiv \\left\\{\\rho_{\\rm HARM},\\tilde{\\tau}_{\\rm HARM},\\tilde{S}_{i}^{\\rm HARM},\\tilde{B}^{i}_{\\rm HARM}\\right\\}\\ .\n", "\\end{align}\n", "$$\n", "\n", "The following relations immediately hold\n", "\n", "$$\n", "\\boxed{\n", "\\begin{align}\n", "\\rho^{\\rm HARM}_{\\star} &= \\rho_{\\star}\\\\\n", "\\tilde{S}_{i}^{\\rm HARM} &= \\tilde{S}_{i}\n", "\\end{align}\n", "}\\ .\n", "$$\n", "\n", "As previously derived in [Step 3.a.ii](#variables_needed_by_harm__bi_harm), we then have\n", "\n", "$$\n", "\\boxed{\\tilde{B}^{i}_{\\rm HARM} = \\underbrace{\\left(\\alpha\\sqrt{\\gamma}\\right)}_{\\rm detg}\\underbrace{\\left(\\frac{B^{i}}{\\alpha\\sqrt{4\\pi}}\\right)}_{\\rm BiL\\_over\\_alpha\\_sqrt\\_fourpi}}\\ .\n", "$$\n", "\n", "Finally, we must relate $\\tilde{\\tau}_{\\rm HARM} = \\sqrt{-g}\\left(T^{0}_{\\ 0}+\\rho_{b}u^{0}\\right)$ to the `IllinoisGRMHD` variables. First, consider the identity\n", "\n", "$$\n", "\\begin{align}\n", "S^{i} &= -\\gamma^{i}_{\\ \\mu}n_{\\nu}T^{\\mu\\nu}\\\\\n", " &= -\\left(\\delta^{i}_{\\mu} + n^{i}n_{\\mu}\\right)n_{\\nu}T^{\\mu\\nu}\\\\\n", " &= -\\delta^{i}_{\\mu}n_{\\nu}T^{\\mu\\nu} - n^{i}n_{\\mu}n_{\\nu}T^{\\mu\\nu}\\\\\n", " &= -n_{0}T^{i0} - n^{i}n_{0}n_{0}T^{00}\\\\\n", " &= \\alpha T^{i0} - \\left(-\\frac{\\beta^{i}}{\\alpha}\\right)\\alpha^{2}T^{00}\\\\\n", "\\implies S^{i}&= \\alpha\\left[\\beta^{i}T^{00} + T^{i0}\\right]\\ .\n", "\\end{align}\n", "$$\n", "\n", "Consider, also, the identity\n", "\n", "$$\n", "\\tilde{\\tau} = \\alpha^{2}\\sqrt{\\gamma}T^{00} - \\rho_{\\star} \\implies T^{00} = \\frac{\\left(\\tilde\\tau + \\rho_{\\star}\\right)}{\\alpha^{2}\\sqrt{\\gamma}}\\ .\n", "$$\n", "\n", "Then\n", "\n", "$$\n", "\\begin{align}\n", "T^{0}_{0} &= g_{0\\mu}T^{0\\mu}\\\\\n", " &= g_{00}T^{00} + g_{0i}T^{0i}\\\\\n", " &= \\left(-\\alpha^{2}+\\beta_{\\ell}\\beta^{\\ell}\\right)T^{00} + \\beta_{i}T^{0i}\\\\\n", " &= -\\alpha^{2}T^{00} + \\beta_{i}\\left[\\beta^{i}T^{00} + T^{i0}\\right]\\\\\n", " &= -\\alpha^{2}\\left[\\frac{\\left(\\tilde\\tau + \\rho_{\\star}\\right)}{\\alpha^{2}\\sqrt{\\gamma}}\\right] + \\frac{\\beta^{i}S_{i}}{\\alpha}\\\\\n", "\\implies T^{0}_{0} &= -\\frac{\\tilde\\tau+\\rho_{\\star}}{\\sqrt{\\gamma}} + \\frac{\\beta^{i}S_{i}}{\\alpha}\\ .\n", "\\end{align}\n", "$$\n", "\n", "Finally, we can compute\n", "\n", "$$\n", "\\begin{align}\n", "\\tilde\\tau_{\\rm HARM} &= \\sqrt{-g}\\left(T^{0}_{0} + \\rho_{b}u^{0}\\right)\\\\\n", " &= \\alpha\\sqrt{\\gamma}T^{0}_{0} + \\alpha\\sqrt{\\gamma}\\rho_{b}u^{0}\\\\\n", " &= \\alpha\\sqrt{\\gamma}\\left[-\\frac{\\tilde\\tau+\\rho_{\\star}}{\\sqrt{\\gamma}} + \\frac{\\beta^{i}S_{i}}{\\alpha}\\right] + \\rho_{\\star}\\\\\n", " &= -\\alpha\\tilde{\\tau} -\\alpha\\rho_{\\star} + \\beta^{i}\\left(\\sqrt{\\gamma}S_{i}\\right) + \\rho_{\\star}\\\\\n", "\\implies &\\boxed{\\tilde\\tau_{\\rm HARM} = -\\alpha\\tilde{\\tau} - \\left(\\alpha-1\\right)\\rho_{\\star} + \\beta^{i}\\tilde{S}_{i}}\\ .\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " // Fill the array of conserved variables according to the wishes of Utoprim_2d.\n", " U[RHO] = CONSERVS[RHOSTAR];\n", " U[UU] = -CONSERVS[TAUENERGY]*METRIC_LAP_PSI4[LAPSE] - (METRIC_LAP_PSI4[LAPSE]-1.0)*CONSERVS[RHOSTAR] +\n", " METRIC[SHIFTX]*CONSERVS[STILDEX] + METRIC[SHIFTY]*CONSERVS[STILDEY] + METRIC[SHIFTZ]*CONSERVS[STILDEZ] ; // note the minus sign on tau\n", " U[UTCON1] = CONSERVS[STILDEX];\n", " U[UTCON2] = CONSERVS[STILDEY];\n", " U[UTCON3] = CONSERVS[STILDEZ];\n", " U[BCON1] = detg*BxL_over_alpha_sqrt_fourpi;\n", " U[BCON2] = detg*ByL_over_alpha_sqrt_fourpi;\n", " U[BCON3] = detg*BzL_over_alpha_sqrt_fourpi;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.a.vii: Writing $\\boldsymbol{P}_{\\rm HARM}$ in terms of $\\boldsymbol{P}_{\\rm IGM}$ \\[Back to [top](#toc)\\]\n", "$$\\label{variables_needed_by_harm__prims}$$\n", "\n", "Keep in mind that for the primitive variables, we are not intereted in providing exact relations. Instead, what we are interested in doing is providing an educated guess of what it should look like in hopes that the Newton-Raphson method then converges to the correct values. With that in mind, the primitive variables are then initialized to:\n", "\n", "$$\n", "\\boxed{\n", "\\begin{align}\n", "\\rho_{\\rm HARM} &= \\rho_{b}\\\\\n", "u &= \\frac{P}{\\Gamma_{\\rm poly} - 1}\\\\\n", "\\tilde{u}^{i} &= u^{0}\\left(v^{i}+\\beta^{i}\\right)\\\\\n", "B^{i}_{\\rm HARM} &= \\frac{B^{i}}{\\alpha\\sqrt{4\\pi}}\n", "\\end{align}\\ ,\n", "}\n", "$$\n", "\n", "where $\\Gamma_{\\rm poly}$ stands for the *local* polytropic $\\Gamma$-factor." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " CCTK_REAL uL = P_oldL/(Gamma_ppoly_tab - 1.0);\n", " CCTK_REAL utxL = u0L*(vxL + METRIC[SHIFTX]);\n", " CCTK_REAL utyL = u0L*(vyL + METRIC[SHIFTY]);\n", " CCTK_REAL utzL = u0L*(vzL + METRIC[SHIFTZ]);\n", "\n", " prim[RHO] = rho_b_oldL;\n", " prim[UU] = uL;\n", " prim[UTCON1] = utxL;\n", " prim[UTCON2] = utyL;\n", " prim[UTCON3] = utzL;\n", " prim[BCON1] = BxL_over_alpha_sqrt_fourpi;\n", " prim[BCON2] = ByL_over_alpha_sqrt_fourpi;\n", " prim[BCON3] = BzL_over_alpha_sqrt_fourpi;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.b: Calling the `HARM` conservative-to-primitive solver \\[Back to [top](#toc)\\]\n", "$$\\label{calling_harm_conservs_to_prims_solver}$$" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " /*************************************************************/\n", " // CALL HARM PRIMITIVES SOLVER:\n", " check = Utoprim_2d(eos, U, g4dn, g4up, detg, prim,stats.n_iter);\n", " // Note that we have modified this solver, so that nearly 100%\n", " // of the time it yields either a good root, or a root with\n", " // negative epsilon (i.e., pressure).\n", " /*************************************************************/" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.c: Applying the Font *et al.* fix, if the inversion fails \\[Back to [top](#toc)\\]\n", "$$\\label{font_fix}$$\n", "\n", "If the conservative-to-primitive solver fails to converge, we apply the procedure suggested by [Font *et al.* (1998)](https://arxiv.org/abs/gr-qc/9811015). The algorithm can be summarized as the requirement that\n", "\n", "$$\n", "P=P_{\\rm cold} = \\kappa \\rho^{\\Gamma_{\\rm cold}}_{b}\\ ,\n", "$$\n", "\n", "and then recomputing the velocities $u_{i}$. We will describe this procedure in detail in [Step 4](#font_fix_gamma_law__c) below." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " // Use the new Font fix subroutine\n", " int font_fix_applied=0;\n", " if(check!=0) {\n", " font_fix_applied=1;\n", " CCTK_REAL u_xl=1e100, u_yl=1e100, u_zl=1e100; // Set to insane values to ensure they are overwritten.\n", " /************************\n", " * New Font fix routine *\n", " ************************/\n", " check = font_fix__hybrid_EOS(u_xl,u_yl,u_zl, CONSERVS,PRIMS,METRIC_PHYS,METRIC_LAP_PSI4, eos);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.d: Compute $\\tilde{u}^{i}$ \\[Back to [top](#toc)\\]\n", "$$\\label{compute_utconi}$$\n", "\n", "Now we evaluate\n", "\n", "$$\n", "\\tilde{u}^{i} = \\gamma^{ij}u_{i} \\implies\n", "\\boxed{\n", "\\left\\{\n", "\\begin{align}\n", "\\tilde{u}^{x} &= \\gamma^{xx}u_{x} + \\gamma^{xy}u_{y} + \\gamma^{xz}u_{z}\\\\\n", "\\tilde{u}^{y} &= \\gamma^{yx}u_{x} + \\gamma^{yy}u_{y} + \\gamma^{yz}u_{z}\\\\\n", "\\tilde{u}^{z} &= \\gamma^{zx}u_{x} + \\gamma^{zy}u_{y} + \\gamma^{zz}u_{z}\n", "\\end{align}\n", "\\right.\n", "}\n", "$$" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", " //Translate to HARM primitive now:\n", " prim[UTCON1] = METRIC_PHYS[GUPXX]*u_xl + METRIC_PHYS[GUPXY]*u_yl + METRIC_PHYS[GUPXZ]*u_zl;\n", " prim[UTCON2] = METRIC_PHYS[GUPXY]*u_xl + METRIC_PHYS[GUPYY]*u_yl + METRIC_PHYS[GUPYZ]*u_zl;\n", " prim[UTCON3] = METRIC_PHYS[GUPXZ]*u_xl + METRIC_PHYS[GUPYZ]*u_yl + METRIC_PHYS[GUPZZ]*u_zl;\n", " if (check==1) {\n", " CCTK_VInfo(CCTK_THORNSTRING,\"Font fix failed!\");\n", " CCTK_VInfo(CCTK_THORNSTRING,\"i,j,k = %d %d %d, stats.failure_checker = %d x,y,z = %e %e %e , index=%d st_i = %e %e %e, rhostar = %e, Bi = %e %e %e, gij = %e %e %e %e %e %e, Psi6 = %e\",i,j,k,stats.failure_checker,X[index],Y[index],Z[index],index,mhd_st_x_orig,mhd_st_y_orig,mhd_st_z_orig,rho_star_orig,PRIMS[BX_CENTER],PRIMS[BY_CENTER],PRIMS[BZ_CENTER],METRIC_PHYS[GXX],METRIC_PHYS[GXY],METRIC_PHYS[GXZ],METRIC_PHYS[GYY],METRIC_PHYS[GYZ],METRIC_PHYS[GZZ],METRIC_LAP_PSI4[PSI6]);\n", " }\n", " }\n", " stats.failure_checker+=font_fix_applied*10000;\n", " stats.font_fixed=font_fix_applied;\n", " /*************************************************************/" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.e: Limiting velocities \\[Back to [top](#toc)\\]\n", "$$\\label{limiting_velocities}$$\n", "\n", "The procedure we follow here is similar to the one discussed in the [inlined functions tutorial module](#Tutorial-IllinoisGRMHD__inlined_functions.ipynb)." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " if(check==0) {\n", " //Now that we have found some solution, we first limit velocity:\n", " //FIXME: Probably want to use exactly the same velocity limiter function here as in mhdflux.C\n", " CCTK_REAL utx_new = prim[UTCON1];\n", " CCTK_REAL uty_new = prim[UTCON2];\n", " CCTK_REAL utz_new = prim[UTCON3];\n", "\n", " //Velocity limiter:\n", " CCTK_REAL gijuiuj = METRIC_PHYS[GXX]*SQR(utx_new ) +\n", " 2.0*METRIC_PHYS[GXY]*utx_new*uty_new + 2.0*METRIC_PHYS[GXZ]*utx_new*utz_new +\n", " METRIC_PHYS[GYY]*SQR(uty_new) + 2.0*METRIC_PHYS[GYZ]*uty_new*utz_new +\n", " METRIC_PHYS[GZZ]*SQR(utz_new);\n", " CCTK_REAL au0m1 = gijuiuj/( 1.0+sqrt(1.0+gijuiuj) );\n", " u0L = (au0m1+1.0)*METRIC_LAP_PSI4[LAPSEINV];\n", "\n", " // *** Limit velocity\n", " stats.vel_limited=0;\n", " if (au0m1 > 0.9999999*(GAMMA_SPEED_LIMIT-1.0)) {\n", " CCTK_REAL fac = sqrt((SQR(GAMMA_SPEED_LIMIT)-1.0)/(SQR(1.0+au0m1) - 1.0));\n", " utx_new *= fac;\n", " uty_new *= fac;\n", " utz_new *= fac;\n", " gijuiuj = gijuiuj * SQR(fac);\n", " au0m1 = gijuiuj/( 1.0+sqrt(1.0+gijuiuj) );\n", " // Reset rho_b and u0\n", " u0L = (au0m1+1.0)*METRIC_LAP_PSI4[LAPSEINV];\n", " prim[RHO] = rho_star_orig/(METRIC_LAP_PSI4[LAPSE]*u0L*METRIC_LAP_PSI4[PSI6]);\n", " stats.vel_limited=1;\n", " stats.failure_checker+=1000;\n", " } //Finished limiting velocity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.f: Setting the primitives \\[Back to [top](#toc)\\]\n", "$$\\label{primitives}$$\n", "\n", "Finally, we update the primitives,\n", "\n", "$$\n", "\\boxed{\n", "\\begin{align}\n", "\\rho_{b} &= \\rho\\\\\n", "P &= \\left(\\Gamma - 1\\right)u = P_{\\rm cold}\\\\\n", "v^{i} &= \\frac{\\tilde{u}^{i}}{u^{0}} - \\beta^{i}\n", "\\end{align}\n", "}\n", "$$" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_primitives_lowlevel.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_primitives_lowlevel__C\n", "\n", "\n", " //The Font fix only sets the velocities. Here we set the pressure & density HARM primitives.\n", " if(font_fix_applied==1) {\n", " prim[RHO] = rho_star_orig/(METRIC_LAP_PSI4[LAPSE]*u0L*METRIC_LAP_PSI4[PSI6]);\n", " //Next set P = P_cold:\n", " CCTK_REAL P_cold;\n", "\n", " /**********************************\n", " * Piecewise Polytropic EOS Patch *\n", " * Finding Gamma_ppoly_tab and K_ppoly_tab *\n", " **********************************/\n", " /* Here we use our newly implemented\n", " * find_polytropic_K_and_Gamma() function\n", " * to determine the relevant polytropic\n", " * Gamma and K parameters to be used\n", " * within this function.\n", " */\n", " int polytropic_index = find_polytropic_K_and_Gamma_index(eos,prim[RHO]);\n", " K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];\n", " Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];\n", "\n", " // After that, we compute P_cold\n", " P_cold = K_ppoly_tab*pow(prim[RHO],Gamma_ppoly_tab);\n", "\n", " prim[UU] = P_cold/(Gamma_ppoly_tab-1.0);\n", " } //Finished setting remaining primitives if there was a Font fix.\n", "\n", " /* Set rho_b */\n", " PRIMS[RHOB] = prim[RHO];\n", "\n", " /***************\n", " * PPEOS Patch *\n", " * Hybrid EOS *\n", " ***************\n", " */\n", " /* We now compute the pressure as a function\n", " * of rhob, P_cold, eps_cold, and u = rhob*eps,\n", " * using the function pressure_rho0_u(), which\n", " * implements the equation:\n", " * .-------------------------------------------------------------.\n", " * | p(rho_b,u) = P_cold + (Gamma_th - 1)*(u - rho_b * eps_cold) |\n", " * .-------------------------------------------------------------.\n", " */\n", " PRIMS[PRESSURE] = pressure_rho0_u(eos, prim[RHO],prim[UU]);\n", "\n", " /* Already set u0L. */\n", " PRIMS[VX] = utx_new/u0L - METRIC[SHIFTX];\n", " PRIMS[VY] = uty_new/u0L - METRIC[SHIFTY];\n", " PRIMS[VZ] = utz_new/u0L - METRIC[SHIFTZ];\n", "\n", " return 0;\n", " } else {\n", " //If we didn't find a root, then try again with a different guess.\n", " }\n", " }\n", " CCTK_VInfo(CCTK_THORNSTRING,\"Couldn't find root from: %e %e %e %e %e, rhob approx=%e, rho_b_atm=%e, Bx=%e, By=%e, Bz=%e, gij_phys=%e %e %e %e %e %e, alpha=%e\",\n", "\t tau_orig,rho_star_orig,mhd_st_x_orig,mhd_st_y_orig,mhd_st_z_orig,rho_star_orig/METRIC_LAP_PSI4[PSI6],rho_b_atm,PRIMS[BX_CENTER],PRIMS[BY_CENTER],PRIMS[BZ_CENTER],METRIC_PHYS[GXX],METRIC_PHYS[GXY],METRIC_PHYS[GXZ],METRIC_PHYS[GYY],METRIC_PHYS[GYZ],METRIC_PHYS[GZZ],METRIC_LAP_PSI4[LAPSE]);\n", " return 1;\n", "}\n", "\n", "//#include \"harm_u2p_util.c\"\n", "#include \"harm_utoprim_2d.c\"\n", "#include \"eigen.C\"\n", "#include \"font_fix_hybrid_EOS.C\"\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 4: `font_fix_hybrid_EOS.C` \\[Back to [top](#toc)\\]\n", "$$\\label{font_fix_hybrid_eos}$$\n", "\n", "### Polytropic EOSs\n", "\n", "The [Font *et al.*](https://arxiv.org/pdf/gr-qc/9811015.pdf) algorithm (henceforth Font Fix algorithm) can be summarized as follows. Font fixes occur at the atmospheric region, so we start by assuming that $P$ is given only by its cold part, i.e.\n", "\n", "$$\n", "P = P_{\\rm cold} = K_{\\rm atm} \\rho_{b}^{\\Gamma_{\\rm atm}}\\ ,\n", "$$\n", "\n", "where $K_{\\rm atm}$ and $\\Gamma_{\\rm atm}$ are the constants used in the atmosphere. Then, the specific internal energy is given by\n", "\n", "$$\n", "\\begin{align}\n", "\\epsilon &= \\epsilon_{\\rm cold}\\\\\n", " &= \\int d\\rho \\frac{P_{\\rm cold}}{\\rho^{2}}\\\\\n", " &= K_{\\rm atm}\\int d\\rho \\rho^{\\Gamma_{\\rm atm}-2}\\\\\n", " &= \\frac{K_{\\rm atm}\\rho^{\\Gamma_{\\rm atm}-1}}{\\Gamma_{\\rm atm}-1}\\ .\n", "\\end{align}\n", "$$\n", "\n", "Having computed $P$ and $\\epsilon$, we can compute the enthalpy, $h$, giving\n", "\n", "$$\n", "\\begin{align}\n", "h &= 1 + \\epsilon + \\frac{P}{\\rho}\\\\\n", " &= 1 + \\frac{K_{\\rm atm}\\rho^{\\Gamma_{\\rm atm}-1}}{\\Gamma_{\\rm atm}-1} + K_{\\rm atm}\\rho^{\\Gamma_{\\rm atm} - 1}\\\\\n", " &= 1 + \\left(\\frac{1}{\\Gamma_{\\rm atm}-1}+1\\right)K_{\\rm atm}\\rho^{\\Gamma_{\\rm atm} - 1}\\\\\n", "\\implies &\\boxed{ h = 1 + \\left(\\frac{\\Gamma_{\\rm atm}}{\\Gamma_{\\rm atm}-1}\\right)K_{\\rm atm} \\rho^{\\Gamma_{\\rm atm}-1} }\\ .\n", "\\end{align}\n", "$$\n", "\n", "We then run an iterative process that updates $\\rho$ in a consistent way, based on the value of $h$. We now describe this process and its implementation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 4.a: Computing the basic quantities needed by the algorithm \\[Back to [top](#toc)\\]\n", "$$\\label{font_fix_hybrid_eos__basic_quantities}$$\n", "\n", "We start by computing all basic quantities needed by the Font Fix algorithm:\n", "\n", "$$\n", "\\boxed{\n", "\\begin{align}\n", "\\bar{B}^{i} &= \\frac{B^i}{\\sqrt{4\\pi}}\\\\\n", "\\bar{B}_{i} &= \\gamma_{ij}\\bar{B}^{j}\\\\\n", "\\bar{B}^{2} &= \\bar{B}_{i}\\bar{B}^{i}\\\\\n", "\\bar{B} &= \\sqrt{\\bar{B}^{2}}\\\\\n", "\\bar{B}\\cdot\\tilde{S} &= \\bar{B}^{i}\\tilde{S}_{i}\\\\\n", " (\\bar{B}&\\cdot\\tilde{S})^{2}\\\\\n", "\\hat{\\bar{B}}\\cdot\\tilde{S} &= \\hat{\\bar{B}}^{i}\\tilde{S}_{i} \\equiv \\left(\\frac{\\bar{B}^{i}}{\\bar{B}}\\right)\\tilde{S}_{i}\\\\\n", "\\tilde{S}\\cdot\\tilde{S} &= \\gamma^{ij}\\tilde{S}_{i}\\tilde{S}_{j}\n", "\\end{align}\n", "}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing ../src/font_fix_hybrid_EOS.C\n" ] } ], "source": [ "%%writefile $outfile_path__font_fix_hybrid_EOS__C\n", "\n", "\n", "/**********************************\n", " * Piecewise Polytropic EOS Patch *\n", " * Font fix: function call *\n", " **********************************/\n", "inline int font_fix__hybrid_EOS(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4, eos_struct eos) {\n", "\n", " CCTK_REAL Bxbar = PRIMS[BX_CENTER]*ONE_OVER_SQRT_4PI;\n", " CCTK_REAL Bybar = PRIMS[BY_CENTER]*ONE_OVER_SQRT_4PI;\n", " CCTK_REAL Bzbar = PRIMS[BZ_CENTER]*ONE_OVER_SQRT_4PI;\n", " CCTK_REAL Bbar_x = METRIC_PHYS[GXX]*Bxbar + METRIC_PHYS[GXY]*Bybar + METRIC_PHYS[GXZ]*Bzbar;\n", " CCTK_REAL Bbar_y = METRIC_PHYS[GXY]*Bxbar + METRIC_PHYS[GYY]*Bybar + METRIC_PHYS[GYZ]*Bzbar;\n", " CCTK_REAL Bbar_z = METRIC_PHYS[GXZ]*Bxbar + METRIC_PHYS[GYZ]*Bybar + METRIC_PHYS[GZZ]*Bzbar;\n", " CCTK_REAL B2bar = Bxbar*Bbar_x + Bybar*Bbar_y + Bzbar*Bbar_z;\n", " CCTK_REAL Bbar = sqrt(B2bar);\n", "\n", " CCTK_REAL check_B_small = fabs(Bxbar)+fabs(Bybar)+fabs(Bzbar);\n", " if (check_B_small>0 && check_B_small<1.e-150) {\n", " // need to compute B2bar specially to prevent floating-point underflow\n", " CCTK_REAL Bmax = fabs(Bxbar);\n", " if (Bmax < fabs(Bybar)) Bmax=fabs(Bybar);\n", " if (Bmax < fabs(Bzbar)) Bmax=fabs(Bzbar);\n", " CCTK_REAL Bxtmp=Bxbar/Bmax, Bytemp=Bybar/Bmax, Bztemp=Bzbar/Bmax;\n", " CCTK_REAL B_xtemp=Bbar_x/Bmax, B_ytemp=Bbar_y/Bmax, B_ztemp=Bbar_z/Bmax;\n", " Bbar = sqrt(Bxtmp*B_xtemp + Bytemp*B_ytemp + Bztemp*B_ztemp)*Bmax;\n", " }\n", " CCTK_REAL BbardotS = Bxbar*CONSERVS[STILDEX] + Bybar*CONSERVS[STILDEY] + Bzbar*CONSERVS[STILDEZ];\n", " CCTK_REAL BbardotS2 = BbardotS*BbardotS;\n", " CCTK_REAL hatBbardotS = BbardotS/Bbar;\n", " if (Bbar<1.e-300) hatBbardotS = 0.0;\n", " CCTK_REAL Psim6 = 1.0/METRIC_LAP_PSI4[PSI6];\n", "\n", " // Limit hatBbardotS\n", " //CCTK_REAL max_gammav = 100.0;\n", " //CCTK_REAL rhob_max = CONSERVS[RHOSTAR]*Psim6;\n", " //CCTK_REAL hmax = 1.0 + gam_gamm1_kpoly*pow(rhob_max,gam1);\n", " //CCTK_REAL abs_hatBbardotS_max = sqrt(SQR(max_gammav)-1.0)*CONSERVS[RHOSTAR]*hmax;\n", " //if (fabs(hatBbardotS) > abs_hatBbardotS_max) {\n", " // CCTK_REAL fac_reduce = abs_hatBbardotS_max/fabs(hatBbardotS);\n", " // CCTK_REAL hatBbardotS_max = hatBbardotS*fac_reduce;\n", " // CCTK_REAL Bbar_inv = 1.0/Bbar;\n", " // CCTK_REAL hat_Bbar_x = Bbar_x*Bbar_inv;\n", " // CCTK_REAL hat_Bbar_y = Bbar_y*Bbar_inv;\n", " // CCTK_REAL hat_Bbar_z = Bbar_z*Bbar_inv;\n", " // CCTK_REAL sub_fact = hatBbardotS_max - hatBbardotS;\n", " // CONSERVS[STILDEX] += sub_fact*hat_Bbar_x;\n", " // CONSERVS[STILDEY] += sub_fact*hat_Bbar_y;\n", " // CONSERVS[STILDEZ] += sub_fact*hat_Bbar_z;\n", " // hatBbardotS = hatBbardotS_max;\n", " // BbardotS *= fac_reduce;\n", " // BbardotS2 = BbardotS*BbardotS;\n", " //}\n", "\n", " CCTK_REAL sdots = METRIC_PHYS[GUPXX]*SQR(CONSERVS[STILDEX]) + METRIC_PHYS[GUPYY]*SQR(CONSERVS[STILDEY]) + METRIC_PHYS[GUPZZ]*SQR(CONSERVS[STILDEZ])\n", " + 2.0*( METRIC_PHYS[GUPXY]*CONSERVS[STILDEX]*CONSERVS[STILDEY] + METRIC_PHYS[GUPXZ]*CONSERVS[STILDEX]*CONSERVS[STILDEZ]\n", " + METRIC_PHYS[GUPYZ]*CONSERVS[STILDEY]*CONSERVS[STILDEZ]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 4.b: Basic check: looking at $\\tilde{S}^{2}$ \\[Back to [top](#toc)\\]\n", "$$\\label{font_fix_hybrid_eos__sdots}$$\n", "\n", "We start by looking at the dot product $\\tilde{S}^{2}$. Recall that\n", "\n", "$$\n", "\\tilde{S}_{i} = \\left(\\rho_{\\star}h + \\alpha\\sqrt{\\gamma}\\, u^{0}b^{2}\\right)u_{i}-\\alpha\\sqrt{\\gamma}\\, b^{0}b_{i}\\ .\n", "$$\n", "\n", "If $\\tilde{S}^{2} = 0$, then we must be in a region where $u_{i} = 0 = b_{i}$. In this case, we return\n", "\n", "$$\n", "\\begin{align}\n", "\\rho_{b} &= \\psi^{-6}\\rho_{\\star}\\ ,\\\\\n", "u^{i} &= 0\\ ,\n", "\\end{align}\n", "$$\n", "\n", "and terminate the function call." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/font_fix_hybrid_EOS.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__font_fix_hybrid_EOS__C\n", "\n", "\n", " CCTK_REAL rhob;\n", " if (sdots<1.e-300) {\n", " rhob = CONSERVS[RHOSTAR]*Psim6;\n", " u_x=0.0; u_y=0.0; u_z=0.0;\n", " return 0;\n", " }\n", " /* This test has some problem.\n", " if (fabs(BbardotS2 - sdots*B2bar) > 1e-8) {\n", " CCTK_VInfo(CCTK_THORNSTRING,\"(Bbar dot S)^2, Bbar^2 * sdotS, %e %e\",SQR(BbardotS),sdots*B2bar);\n", " CCTK_VInfo(CCTK_THORNSTRING,\"Cauchy-Schwartz inequality is violated!\");\n", " }\n", " */" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 4.c: Initial guesses for $W$, $S_{{\\rm fluid}}^{2}$, and $\\rho$ \\[Back to [top](#toc)\\]\n", "$$\\label{font_fix_hybrid_eos__initial_guesses}$$\n", "\n", "If $\\tilde{S}^{2} \\neq 0$, then we move on to the iterative procedure previously mentioned. We start by setting the initial data based on eqs. (A52), (A53), and (A59) found in [Appendix A of Zachariah *et al.* (2012)](https://arxiv.org/pdf/1112.0568.pdf):\n", "\n", "$$\n", "\\boxed{\n", "\\begin{align}\n", "W_{0} &= \\psi^{-6}\\sqrt{\\left(\\hat{\\bar{B}}\\cdot\\tilde{S}\\right)^{2} + \\rho_{\\star}^{2}}\\\\\n", "S_{{\\rm fluid},0}^{2} &=\\frac{W_{0}^{2}\\left(\\tilde{S}\\cdot\\tilde{S}\\right)+\\left(\\bar{B}\\cdot\\tilde{S}\\right)^{2}\\left(\\bar{B}^{2} + 2W_{0}\\right)}{\\left(W_{0} + \\bar{B}^{2}\\right)^{2}}\\\\\n", "\\rho_{0} &= \\frac{\\psi^{-6}\\rho_{\\star}}{\\sqrt{1+\\frac{S_{{\\rm fluid},0}^{2}}{\\rho_{\\star}^{2}}}}\n", "\\end{align}\n", "}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/font_fix_hybrid_EOS.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__font_fix_hybrid_EOS__C\n", "\n", "\n", " // Initial guess for W, S_fluid and rhob\n", " CCTK_REAL W0 = sqrt( SQR(hatBbardotS) + SQR(CONSERVS[RHOSTAR]) ) * Psim6;\n", " CCTK_REAL Sf20 = (SQR(W0)*sdots + BbardotS2*(B2bar + 2.0*W0))/SQR(W0+B2bar);\n", " CCTK_REAL rhob0 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 4.d: The main loop \\[Back to [top](#toc)\\]\n", "$$\\label{font_fix_hybrid_eos__main_loop}$$\n", "\n", "We now perform the following iterative process, which is again described in [Appendix A of Zachariah *et al.* (2012)](https://arxiv.org/pdf/1112.0568.pdf). We refer the reader to eqs. (A60), (A61), and (A62).\n", "\n", "1. Store the previously computed values of $W_{n}$, $S_{{\\rm fluid},n}^{2}$, and $\\rho_{n}$\n", "2. Compute $h = 1 + \\epsilon_{\\rm cold} + P_{\\rm cold}/\\rho_{n}$\n", "3. Set\n", "$$\n", "\\boxed{\\rho_{n+1} = \\psi^{-6}\\rho_{\\star}\\left(1 + \\frac{S_{{\\rm fluid},n}^{2}}{\\rho_{\\star} h_{n}}\\right)^{-1/2}}\n", "$$\n", "\n", "4. For a given value of $n$, perform steps 1 (for $\\rho$), 2 and 3 until $\\left|\\rho_{n+1}-\\rho_{n}\\right| < \\rho_{n+1}\\epsilon$, where $\\epsilon$ is a user given tolerance\n", "5. After convergence is obtained, update:\n", "$$\n", "\\boxed{\n", "\\begin{align}\n", "h_{n+1} &= 1 + \\epsilon_{\\rm cold} + P_{\\rm cold}/\\rho_{n+1}\\\\\n", "W_{n+1} &= \\psi^{-6}\\sqrt{\\tilde{S}^{2}_{{\\rm fluid},n} + \\rho_{\\star}^{2} h_{n+1}^{2}}\\\\\n", "S_{{\\rm fluid},n+1}^{2} &= \\frac{W^{2}_{n+1}\\left(\\tilde{S}\\cdot\\tilde{S}\\right) + \\left(\\bar{B}\\cdot\\tilde{S}\\right)^{2}\\left(\\bar{B}^{2} + 2W_{n+1}\\right)}{\\left(W_{n+1} + \\bar{B}^{2}\\right)^{2}}\n", "\\end{align}\n", "}\\ .\n", "$$\n", "6. Repeat steps 1 through 5 until $\\left|W_{n+1}-W_{n}\\right| < W_{n+1}\\epsilon$ *and* $\\left|S^{2}_{{\\rm fluid},n+1}-S^{2}_{{\\rm fluid},n}\\right| < S^{2}_{{\\rm fluid},n+1}\\epsilon$ *or* we reach the maximum number of iterations\n", "7. If font fix fails, increase the tolerance and try again.\n", "\n", "This is done using the function `font_fix__rhob_loop()`, which is documented in the [inlined functions tutorial notebook](Tutorial-IllinoisGRMHD__inlined_functions.ipynb)." ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/font_fix_hybrid_EOS.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__font_fix_hybrid_EOS__C\n", "\n", "\n", " //****************************************************************\n", " // FONT FIX\n", " // Impose Font fix when HARM primitives solver fails to find\n", " // acceptable set of primitives.\n", " //****************************************************************\n", "\n", " /* Set the maximum number of iterations */\n", " int maxits = 500;\n", "\n", " /* Set the allowed tolerance */\n", " CCTK_REAL tol = 1.e-15;\n", "\n", " /* Declare basic variables */\n", " int font_fix_status;\n", "\n", " /**********************\n", " * FONT FIX MAIN LOOP *\n", " **********************\n", " * Perform the font fix routine until convergence\n", " * is obtained and the algorithm returns with no\n", " * error. Every time the Font fix fails, increase\n", " * the tolerance by a factor of 10.\n", " */\n", " int font_fix_attempts = 5;\n", " CCTK_REAL font_fix_tol_factor = 10.0;\n", " for(int n=0; n\n", "\n", "## Step 4.e: Output $\\rho_{b}$ and $u_{i}$ \\[Back to [top](#toc)\\]\n", "$$\\label{font_fix_hybrid_eos__outputs}$$\n", "\n", "Finally, we return $\\rho_{b}$ and $u_{i}$. In the relations below, $N$ indicates the last computed value of $\\rho$ obtained by our iterative process. The quantities evaluated here are\n", "\n", "$$\n", "\\boxed{\n", "\\begin{align}\n", "\\rho_{b} &= \\rho_{N}\\\\\n", "\\gamma_{v} &= \\frac{\\psi^{-6}\\rho_{\\star}}{\\rho_{b}}\\\\\n", "f_{1} &= \\frac{\\psi^{6}\\left(\\bar{B}\\cdot\\tilde{S}\\right)}{\\gamma_{v}\\rho_{\\star} h}\\\\\n", "f_{2} &= \\left(\\rho_{\\star}h + \\psi^{6}\\frac{\\bar{B}^{2}}{\\gamma_{v}}\\right)^{-1}\\\\\n", "u_{i} &= f_{2}\\left(\\tilde{S}_{i} + f_{1}\\bar{B}_{i}\\right)\n", "\\end{align}\n", "}\\ .\n", "$$\n" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/font_fix_hybrid_EOS.C\n" ] } ], "source": [ "%%writefile -a $outfile_path__font_fix_hybrid_EOS__C\n", "\n", "\n", " //**************************************************************************************************************\n", "\n", " /* Font fix works! */\n", " /* First compute P_cold, eps_cold, then h = h_cold */\n", " CCTK_REAL P_cold, eps_cold;\n", " compute_P_cold__eps_cold(eos,rhob, P_cold,eps_cold);\n", " CCTK_REAL h = 1.0 + eps_cold + P_cold/rhob;\n", "\n", " /* Then compute gamma_v using equation (A19) in\n", " * Etienne et al. (2011) [https://arxiv.org/pdf/1112.0568.pdf]\n", " * .-----------------------------------------.\n", " * | gamma_v = psi^{-6} * (rho_star / rho_b) |\n", " * .-----------------------------------------.\n", " */\n", " CCTK_REAL gammav = CONSERVS[RHOSTAR]*Psim6/rhob;\n", "\n", " /* Finally, compute u_{i} */\n", " CCTK_REAL rhosh = CONSERVS[RHOSTAR]*h;\n", " CCTK_REAL fac1 = METRIC_LAP_PSI4[PSI6]*BbardotS/(gammav*rhosh);\n", " CCTK_REAL fac2 = 1.0/(rhosh + METRIC_LAP_PSI4[PSI6]*B2bar/gammav);\n", " u_x = fac2*(CONSERVS[STILDEX] + fac1*Bbar_x);\n", " u_y = fac2*(CONSERVS[STILDEY] + fac1*Bbar_y);\n", " u_z = fac2*(CONSERVS[STILDEZ] + fac1*Bbar_z);\n", "\n", " return 0;\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 5: `harm_primitives_headers.h` \\[Back to [top](#toc)\\]\n", "$$\\label{harm_primitives_headers}$$" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing ../src/harm_primitives_headers.h\n" ] } ], "source": [ "%%writefile $outfile_path__harm_primitives_headers__h\n", "/***********************************************************************************\n", " Copyright 2006 Charles F. Gammie, Jonathan C. McKinney, Scott C. Noble,\n", " Gabor Toth, and Luca Del Zanna\n", "\n", " HARM version 1.0 (released May 1, 2006)\n", "\n", " This file is part of HARM. HARM is a program that solves hyperbolic\n", " partial differential equations in conservative form using high-resolution\n", " shock-capturing techniques. This version of HARM has been configured to\n", " solve the relativistic magnetohydrodynamic equations of motion on a\n", " stationary black hole spacetime in Kerr-Schild coordinates to evolve\n", " an accretion disk model.\n", "\n", " You are morally obligated to cite the following two papers in his/her\n", " scientific literature that results from use of any part of HARM:\n", "\n", " [1] Gammie, C. F., McKinney, J. C., \\& Toth, G.\\ 2003,\n", " Astrophysical Journal, 589, 444.\n", "\n", " [2] Noble, S. C., Gammie, C. F., McKinney, J. C., \\& Del Zanna, L. \\ 2006,\n", " Astrophysical Journal, 641, 626.\n", "\n", "\n", " Further, we strongly encourage you to obtain the latest version of\n", " HARM directly from our distribution website:\n", " http://rainman.astro.uiuc.edu/codelib/\n", "\n", "\n", " HARM is free software; you can redistribute it and/or modify\n", " it under the terms of the GNU General Public License as published by\n", " the Free Software Foundation; either version 2 of the License, or\n", " (at your option) any later version.\n", "\n", " HARM is distributed in the hope that it will be useful,\n", " but WITHOUT ANY WARRANTY; without even the implied warranty of\n", " MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n", " GNU General Public License for more details.\n", "\n", " You should have received a copy of the GNU General Public License\n", " along with HARM; if not, write to the Free Software\n", " Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA\n", "\n", "***********************************************************************************/\n", "#ifndef HARM_PRIMITIVES_HEADERS_H_\n", "#define HARM_PRIMITIVES_HEADERS_H_\n", "\n", "static const int NPR =8;\n", "static const int NDIM=4;\n", "\n", "/* Adiabatic index used for the state equation */\n", "//#define GAMMA (2.0)\n", "\n", "static const CCTK_REAL G_ISOTHERMAL = 1.0;\n", "\n", "/* use K(s)=K(r)=const. (G_ATM = GAMMA) of time or T = T(r) = const. of time (G_ATM = 1.) */\n", "/*\n", " #define USE_ISENTROPIC 1\n", "\n", " #if( USE_ISENTROPIC )\n", " #define G_ATM GAMMA\n", " #else\n", " #define G_ATM G_ISOTHERMAL\n", " #endif\n", "*/\n", "\n", "static const int MAX_NEWT_ITER=30; /* Max. # of Newton-Raphson iterations for find_root_2D(); */\n", "//#define MAX_NEWT_ITER 300 /* Max. # of Newton-Raphson iterations for find_root_2D(); */\n", "static const CCTK_REAL NEWT_TOL =1.0e-10; /* Min. of tolerance allowed for Newton-Raphson iterations */\n", "static const CCTK_REAL MIN_NEWT_TOL=1.0e-10; /* Max. of tolerance allowed for Newton-Raphson iterations */\n", "static const int EXTRA_NEWT_ITER=0; /* ZACH SAYS: Original value = 2. But I don't think this parameter > 0 is warranted. Just slows the code for no reason, since our tolerances are fine. */\n", "\n", "static const CCTK_REAL NEWT_TOL2 =1.0e-15; /* TOL of new 1D^*_{v^2} gnr2 method */\n", "static const CCTK_REAL MIN_NEWT_TOL2=1.0e-10; /* TOL of new 1D^*_{v^2} gnr2 method */\n", "\n", "static const CCTK_REAL W_TOO_BIG =1.e20; /* \\gamma^2 (\\rho_0 + u + p) is assumed\n", " to always be smaller than this. This\n", " is used to detect solver failures */\n", "static const CCTK_REAL UTSQ_TOO_BIG =1.e20; /* \\tilde{u}^2 is assumed to be smaller\n", " than this. Used to detect solver\n", " failures */\n", "\n", "static const CCTK_REAL FAIL_VAL =1.e30; /* Generic value to which we set variables when a problem arises */\n", "\n", "static const CCTK_REAL NUMEPSILON=(2.2204460492503131e-16);\n", "\n", "\n", "/* some mnemonics */\n", "/* for primitive variables */\n", "static const int RHO =0;\n", "static const int UU =1;\n", "static const int UTCON1 =2;\n", "static const int UTCON2 =3;\n", "static const int UTCON3 =4;\n", "static const int BCON1 =5;\n", "static const int BCON2 =6;\n", "static const int BCON3 =7;\n", "\n", "/* for conserved variables */\n", "static const int QCOV0 =1;\n", "static const int QCOV1 =2;\n", "static const int QCOV2 =3;\n", "static const int QCOV3 =4;\n", "\n", "/********************************************************************************************/\n", "// Function prototype declarations:\n", "int Utoprim_2d(eos_struct eos, CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM],\n", " CCTK_REAL gdet, CCTK_REAL prim[NPR], long &n_iter);\n", "\n", "inline int harm_primitives_gammalaw_lowlevel(const int index,const int i,const int j,const int k,CCTK_REAL *X,CCTK_REAL *Y,CCTK_REAL *Z,\n", " CCTK_REAL *METRIC,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,\n", " CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,\n", " CCTK_REAL g4dn[NDIM][NDIM],CCTK_REAL g4up[NDIM][NDIM],\n", " output_stats &stats,eos_struct &eos);\n", "\n", "inline int font_fix__hybrid_EOS(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4, eos_struct eos);\n", "void eigenvalues_3by3_real_sym_matrix(CCTK_REAL & lam1, CCTK_REAL & lam2, CCTK_REAL & lam3,\n", " CCTK_REAL M11, CCTK_REAL M12, CCTK_REAL M13, CCTK_REAL M22, CCTK_REAL M23, CCTK_REAL M33);\n", "\n", "/********************************************************************************************/\n", "\n", "#endif /* HARM_PRIMITIVES_HEADERS_H_ */\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 6: Code validation \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation}$$\n", "\n", "\n", "\n", "## Step 6.a: `driver_conserv_to_prims.C` \\[Back to [top](#toc)\\]\n", "$$\\label{driver_conserv_to_prims_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": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Validation test for driver_conserv_to_prims.C: FAILED!\n", "Diff:\n", "1c1\n", "< /* We evolve forward in time a set of functions called the \n", "---\n", "> /* We evolve forward in time a set of functions called the\n", "3,4c3,4\n", "< * are updated, we must solve for the primitive variables \n", "< * (rho, pressure, velocities) using a Newton-Raphson \n", "---\n", "> * are updated, we must solve for the primitive variables\n", "> * (rho, pressure, velocities) using a Newton-Raphson\n", "6c6\n", "< * of the MHD equations again. \n", "---\n", "> * of the MHD equations again.\n", "9,12c9,12\n", "< * Raphson solver. Truncation errors in conservative \n", "< * variables can lead to no physical solutions in \n", "< * primitive variables. We correct for these errors here \n", "< * through a number of tricks described in the appendices \n", "---\n", "> * Raphson solver. Truncation errors in conservative\n", "> * variables can lead to no physical solutions in\n", "> * primitive variables. We correct for these errors here\n", "> * through a number of tricks described in the appendices\n", "15,16c15,16\n", "< * This is a wrapper for the 2d solver of Noble et al. See \n", "< * harm_utoprim_2d.c for references and copyright notice \n", "---\n", "> * This is a wrapper for the 2d solver of Noble et al. See\n", "> * harm_utoprim_2d.c for references and copyright notice\n", "19,20c19,20\n", "< * \n", "< * For optimal compatibility, this wrapper is licensed under \n", "---\n", "> *\n", "> * For optimal compatibility, this wrapper is licensed under\n", "23c23\n", "< * Note that this code assumes a simple gamma law for the \n", "---\n", "> * Note that this code assumes a simple gamma law for the\n", "27,28c27\n", "< \n", "< #include \"cctk.h\"\n", "---\n", "> // Standard #include's\n", "32d30\n", "< #include \n", "35a34,39\n", "> \n", "> \n", "> #ifdef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", "> #include \"standalone_conserv_to_prims_main_function.h\"\n", "> #else\n", "> #include \"cctk.h\"\n", "41a46\n", "> #include \"harm_u2p_util.c\"\n", "50a56\n", "> #endif\n", "52,54c58,68\n", "< // FIXME: This is definitely one of the places we need to look at later.\n", "< // Will leave this comment to emphasize the \"FIXME\" below.\n", "< // FIXME: only for single gamma-law EOS.\n", "---\n", "> /**********************************\n", "> * Piecewise Polytropic EOS Patch *\n", "> * Setting up the EOS struct *\n", "> **********************************/\n", "> /*\n", "> * The short piece of code below takes care\n", "> * of initializing the EOS parameters.\n", "> * Please refer to the \"inlined_functions.C\"\n", "> * source file for the documentation on the\n", "> * function.\n", "> */\n", "56,63c70,71\n", "< eos.neos=neos;\n", "< eos.K_poly=K_poly;\n", "< eos.rho_tab[0]=rho_tab[0];\n", "< eos.P_tab[0]=P_tab[0];\n", "< eos.gamma_th=gamma_th;\n", "< eos.eps_tab[0]=eps_tab[0];\n", "< eos.k_tab[0]=k_tab[0]; eos.k_tab[1]=k_tab[1];\n", "< eos.gamma_tab[0]=gamma_tab[0]; eos.gamma_tab[1]=gamma_tab[1];\n", "---\n", "> initialize_EOS_struct_from_input(eos);\n", "> \n", "71a80,81\n", "> \n", "> #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", "82a93,94\n", "> #endif\n", "> \n", "93,95d104\n", "< int gamma_equals2 = 1;\n", "< if (fabs(gamma_th-2.0) > 1.e-10) gamma_equals2 = 0;\n", "< \n", "109a119\n", "> \n", "139a150\n", "> \n", "150a162\n", "> \n", "152a165\n", "> \n", "155c168,169\n", "< \n", "---\n", "> \n", "> \n", "169a184\n", "> \n", "176a192\n", "> \n", "190c206\n", "< \n", "---\n", "> \n", "202a219\n", "> \n", "210c227\n", "< \n", "---\n", "> \n", "223a241\n", "> \n", "245,250c263,270\n", "< PRIMS[RHOB] =rho_b_atm;\n", "< if(gamma_equals2) {\n", "< PRIMS[PRESSURE]=K_poly*SQR(rho_b_atm);\n", "< } else {\n", "< PRIMS[PRESSURE]=K_poly*pow(rho_b_atm,gamma_th);\n", "< }\n", "---\n", "> PRIMS[RHOB] = rho_b_atm;\n", "> \n", "> /* Set P = P_cold */\n", "> int polytropic_index = find_polytropic_K_and_Gamma_index(eos, rho_b_atm);\n", "> CCTK_REAL K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];\n", "> CCTK_REAL Gamma_ppoly_tab = eos.K_ppoly_tab[polytropic_index];\n", "> PRIMS[PRESSURE] = K_ppoly_tab*pow(rho_b_atm,Gamma_ppoly_tab);\n", "> \n", "257a278\n", "> \n", "261a283\n", "> \n", "274a297\n", "> \n", "304c327\n", "< error_int_numer += fabs(tau[index] - tau_orig) + fabs(rho_star[index] - rho_star_orig) + \n", "---\n", "> error_int_numer += fabs(tau[index] - tau_orig) + fabs(rho_star[index] - rho_star_orig) +\n", "341c364\n", "< \n", "---\n", "> \n", "346c369\n", "< // Very useful con2prim debugger. If the primitives (con2prim) solver fails, this will output all data needed to \n", "---\n", "> // Very useful con2prim debugger. If the primitives (con2prim) solver fails, this will output all data needed to\n", "348c371\n", "< if(conserv_to_prims_debug==1) {\n", "---\n", "> if(conserv_to_prims_debug==1 && error_int_numer/error_int_denom > 0.05) {\n", "353,354c376,379\n", "< sprintf(filename,\"primitives_debug-%e.dat\",rho_star[CCTK_GFINDEX3D(cctkGH,3,15,6)]);\n", "< //Alternative, for debugging purposes as well: sprintf(filename,\"primitives_debug-%d.dat\",rand());\n", "---\n", "> sprintf(filename,\"primitives_debug-%e.dat\",error_int_numer/error_int_denom);\n", "> //Alternative, for debugging purposes as well:\n", "> //srand(time(NULL));\n", "> //sprintf(filename,\"primitives_debug-%d.dat\",rand());\n", "358a384,386\n", "> myfile.write((char*)&GAMMA_SPEED_LIMIT, 1*sizeof(CCTK_REAL));\n", "> \n", "> myfile.write((char*)&rho_b_max, 1*sizeof(CCTK_REAL));\n", "364,373c392,398\n", "< CCTK_REAL gamma_th=2.0;\n", "< myfile.write((char*)&gamma_th, 1*sizeof(CCTK_REAL));\n", "< myfile.write((char*)&neos, 1*sizeof(int));\n", "< \n", "< myfile.write((char*)gamma_tab, (neos+1)*sizeof(CCTK_REAL));\n", "< myfile.write((char*)k_tab, (neos+1)*sizeof(CCTK_REAL));\n", "< \n", "< myfile.write((char*)eps_tab, neos*sizeof(CCTK_REAL));\n", "< myfile.write((char*)rho_tab, neos*sizeof(CCTK_REAL));\n", "< myfile.write((char*)P_tab, neos*sizeof(CCTK_REAL));\n", "---\n", "> myfile.write((char*)&update_Tmunu, 1*sizeof(int));\n", "> \n", "> myfile.write((char*)&neos, 1*sizeof(int));\n", "> myfile.write((char*)&Gamma_th, 1*sizeof(CCTK_REAL));\n", "> myfile.write((char*)&K_ppoly_tab0, 1*sizeof(CCTK_REAL));\n", "> myfile.write((char*)Gamma_ppoly_tab_in, neos*sizeof(CCTK_REAL));\n", "> myfile.write((char*)rho_ppoly_tab_in, (neos-1)*sizeof(CCTK_REAL));\n", "378a404,424\n", "> \n", "> myfile.write((char *)failure_checker, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)eTtt, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)eTtx, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)eTty, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)eTtz, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)eTxx, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)eTxy, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)eTxz, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)eTyy, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)eTyz, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)eTzz, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)alp, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)gxx, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)gxy, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)gxz, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)gyy, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)gyz, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)gzz, fullsize*sizeof(CCTK_REAL));\n", "> myfile.write((char *)psi_bssn, fullsize*sizeof(CCTK_REAL));\n", "> \n", "399c445\n", "< \n", "---\n", "> \n", "421c467\n", "< CCTK_VInfo(CCTK_THORNSTRING,\"Finished writing...\");\n", "---\n", "> CCTK_VInfo(CCTK_THORNSTRING,\"Finished writing %s\",filename);\n", "422a469,473\n", "> \n", "> #ifdef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", "> return 0; // int main() requires an integer be returned\n", "> #endif\n", "> \n", "425a477\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/driver_conserv_to_prims.C\"\n", "original_IGM_file_name = \"driver_conserv_to_prims-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__driver_conserv_to_prims__C = !diff $original_IGM_file_path $outfile_path__driver_conserv_to_prims__C\n", "\n", "if Validation__driver_conserv_to_prims__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 driver_conserv_to_prims.C: PASSED!\")\n", "else:\n", " # If the validation fails, we keep the original IGM source code file\n", " print(\"Validation test for driver_conserv_to_prims.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__driver_conserv_to_prims__C:\n", " print(diff_line)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 6.b: `harm_primitives_lowlevel.C` \\[Back to [top](#toc)\\]\n", "$$\\label{harm_primitives_lowlevel_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": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Validation test for harm_primitives_lowlevel.C: FAILED!\n", "Diff:\n", "6c6\n", "< \n", "---\n", "> #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", "7a8\n", "> #endif\n", "9,12d9\n", "< CCTK_REAL kpoly = eos.k_tab[0];\n", "< CCTK_REAL gamma = eos.gamma_tab[0];\n", "< int gamma_equals2 = 1;\n", "< if (fabs(gamma-2.0) > 1.e-10) gamma_equals2 = 0;\n", "15c12\n", "< CCTK_REAL U[NPR]; \n", "---\n", "> CCTK_REAL U[NPR];\n", "42a40\n", "> \n", "46c44\n", "< // This is NOT the \\mathcal{B}^i, which differs by \n", "---\n", "> // This is NOT the \\mathcal{B}^i, which differs by\n", "51a50\n", "> \n", "57a57\n", "> \n", "60,61c60,61\n", "< between the two sets of definitions for U and P. The user may \n", "< wish to alter the translation as they see fit. \n", "---\n", "> between the two sets of definitions for U and P. The user may\n", "> wish to alter the translation as they see fit.\n", "64,72c64,72\n", "< // / rho u^t \\ //\n", "< // U = | T^t_t + rho u^t | sqrt(-det(g_{\\mu\\nu})) //\n", "< // | T^t_i | //\n", "< // \\ B^i / //\n", "< // //\n", "< // / rho \\ //\n", "< // P = | uu | //\n", "< // | \\tilde{u}^i | //\n", "< // \\ B^i / //\n", "---\n", "> // / rho u^t \\ //\n", "> // U = | T^t_t + rho u^t | * sqrt(-det(g_{\\mu\\nu})) //\n", "> // | T^t_i | //\n", "> // \\ B^i / //\n", "> // //\n", "> // / rho \\ //\n", "> // P = | uu | //\n", "> // | \\tilde{u}^i | //\n", "> // \\ B^i / //\n", "76c76\n", "< \n", "---\n", "> \n", "81c81\n", "< // prim[NPR] = primitive variables (guess on input, calculated values on \n", "---\n", "> // prim[NPR] = primitive variables (guess on input, calculated values on\n", "84c84\n", "< // U[1] = \n", "---\n", "> // U[1] =\n", "92a93\n", "> \n", "103a105\n", "> CCTK_REAL K_ppoly_tab,Gamma_ppoly_tab;\n", "111,115c113,130\n", "< if (gamma_equals2==1) {\n", "< P_oldL = kpoly*rho_b_oldL*rho_b_oldL;\n", "< } else {\n", "< P_oldL = kpoly*pow(rho_b_oldL,gamma);\n", "< }\n", "---\n", "> \n", "> /**********************************\n", "> * Piecewise Polytropic EOS Patch *\n", "> * Finding Gamma_ppoly_tab and K_ppoly_tab *\n", "> **********************************/\n", "> /* Here we use our newly implemented\n", "> * find_polytropic_K_and_Gamma() function\n", "> * to determine the relevant polytropic\n", "> * Gamma and K parameters to be used\n", "> * within this function.\n", "> */\n", "> int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_oldL);\n", "> K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];\n", "> Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];\n", "> \n", "> // After that, we compute P_cold\n", "> P_oldL = K_ppoly_tab*pow(rho_b_oldL,Gamma_ppoly_tab);\n", "> \n", "125,130c140,157\n", "< // GAMMA=2 ONLY:\n", "< if (gamma_equals2==1) {\n", "< P_oldL = kpoly*rho_b_oldL*rho_b_oldL;\n", "< } else {\n", "< P_oldL = kpoly*pow(rho_b_oldL,gamma);\n", "< }\n", "---\n", "> \n", "> /**********************************\n", "> * Piecewise Polytropic EOS Patch *\n", "> * Finding Gamma_ppoly_tab and K_ppoly_tab *\n", "> **********************************/\n", "> /* Here we use our newly implemented\n", "> * find_polytropic_K_and_Gamma() function\n", "> * to determine the relevant polytropic\n", "> * Gamma and K parameters to be used\n", "> * within this function.\n", "> */\n", "> int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b_oldL);\n", "> K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];\n", "> Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];\n", "> \n", "> // After that, we compute P_cold\n", "> P_oldL = K_ppoly_tab*pow(rho_b_oldL,Gamma_ppoly_tab);\n", "> \n", "136a164\n", "> \n", "139c167\n", "< U[UU] = -CONSERVS[TAUENERGY]*METRIC_LAP_PSI4[LAPSE] - (METRIC_LAP_PSI4[LAPSE]-1.0)*CONSERVS[RHOSTAR] + \n", "---\n", "> U[UU] = -CONSERVS[TAUENERGY]*METRIC_LAP_PSI4[LAPSE] - (METRIC_LAP_PSI4[LAPSE]-1.0)*CONSERVS[RHOSTAR] +\n", "148c176,177\n", "< CCTK_REAL uL = P_oldL/(gamma_th /* <- Should be local polytropic gamma factor */ - 1.0);\n", "---\n", "> \n", "> CCTK_REAL uL = P_oldL/(Gamma_ppoly_tab - 1.0);\n", "161a191\n", "> \n", "164,166c194,196\n", "< check = Utoprim_2d(U, g4dn, g4up, detg, prim,stats.n_iter);\n", "< // Note that we have modified this solver, so that nearly 100% \n", "< // of the time it yields either a good root, or a root with \n", "---\n", "> check = Utoprim_2d(eos, U, g4dn, g4up, detg, prim,stats.n_iter);\n", "> // Note that we have modified this solver, so that nearly 100%\n", "> // of the time it yields either a good root, or a root with\n", "170c200,201\n", "< // Use the new Font fix subroutine \n", "---\n", "> \n", "> // Use the new Font fix subroutine\n", "175,179c206,210\n", "< if (gamma_equals2==1) {\n", "< check = font_fix_gamma_equals2(u_xl,u_yl,u_zl,CONSERVS,PRIMS,METRIC_PHYS,METRIC_LAP_PSI4,eos);\n", "< } else {\n", "< check = font_fix_general_gamma(u_xl,u_yl,u_zl,CONSERVS,PRIMS,METRIC_PHYS,METRIC_LAP_PSI4,eos);\n", "< }\n", "---\n", "> /************************\n", "> * New Font fix routine *\n", "> ************************/\n", "> check = font_fix__hybrid_EOS(u_xl,u_yl,u_zl, CONSERVS,PRIMS,METRIC_PHYS,METRIC_LAP_PSI4, eos);\n", "> \n", "192a224\n", "> \n", "223a256\n", "> \n", "229,234c262,280\n", "< if (gamma_equals2==1) {\n", "< P_cold = kpoly*prim[RHO]*prim[RHO]; // <- FIXME: GAMMA=2 ONLY\n", "< } else {\n", "< P_cold = kpoly*pow(prim[RHO],gamma);\n", "< }\n", "< prim[UU] = P_cold/(gamma_th /* <- Should be local polytropic gamma factor */-1.0);\n", "---\n", "> \n", "> /**********************************\n", "> * Piecewise Polytropic EOS Patch *\n", "> * Finding Gamma_ppoly_tab and K_ppoly_tab *\n", "> **********************************/\n", "> /* Here we use our newly implemented\n", "> * find_polytropic_K_and_Gamma() function\n", "> * to determine the relevant polytropic\n", "> * Gamma and K parameters to be used\n", "> * within this function.\n", "> */\n", "> int polytropic_index = find_polytropic_K_and_Gamma_index(eos,prim[RHO]);\n", "> K_ppoly_tab = eos.K_ppoly_tab[polytropic_index];\n", "> Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];\n", "> \n", "> // After that, we compute P_cold\n", "> P_cold = K_ppoly_tab*pow(prim[RHO],Gamma_ppoly_tab);\n", "> \n", "> prim[UU] = P_cold/(Gamma_ppoly_tab-1.0);\n", "237,239c283,301\n", "< PRIMS[RHOB] = prim[RHO];\n", "< PRIMS[PRESSURE] = (gamma_th /* <- Should be local polytropic gamma factor */-1.0)*prim[UU];\n", "< //Already set u0L.\n", "---\n", "> /* Set rho_b */\n", "> PRIMS[RHOB] = prim[RHO];\n", "> \n", "> /***************\n", "> * PPEOS Patch *\n", "> * Hybrid EOS *\n", "> ***************\n", "> */\n", "> /* We now compute the pressure as a function\n", "> * of rhob, P_cold, eps_cold, and u = rhob*eps,\n", "> * using the function pressure_rho0_u(), which\n", "> * implements the equation:\n", "> * .-------------------------------------------------------------.\n", "> * | p(rho_b,u) = P_cold + (Gamma_th - 1)*(u - rho_b * eps_cold) |\n", "> * .-------------------------------------------------------------.\n", "> */\n", "> PRIMS[PRESSURE] = pressure_rho0_u(eos, prim[RHO],prim[UU]);\n", "> \n", "> /* Already set u0L. */\n", "254c316\n", "< #include \"harm_u2p_util.c\"\n", "---\n", "> //#include \"harm_u2p_util.c\"\n", "257c319,320\n", "< #include \"font_fix_gamma_law.C\"\n", "---\n", "> #include \"font_fix_hybrid_EOS.C\"\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/harm_primitives_lowlevel.C\"\n", "original_IGM_file_name = \"harm_primitives_lowlevel-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__harm_primitives_lowlevel__C = !diff $original_IGM_file_path $outfile_path__harm_primitives_lowlevel__C\n", "\n", "if Validation__harm_primitives_lowlevel__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 harm_primitives_lowlevel.C: PASSED!\")\n", "else:\n", " # If the validation fails, we keep the original IGM source code file\n", " print(\"Validation test for harm_primitives_lowlevel.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__harm_primitives_lowlevel__C:\n", " print(diff_line)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 6.c: `font_fix_gamma_law.C` \\[Back to [top](#toc)\\]\n", "$$\\label{font_fix_gamma_law_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": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Validation test for font_fix_gamma_law.C: FAILED!\n", "Diff:\n", "1,10d0\n", "< /*-----------------------------------------------------------------------------\n", "< *\n", "< * Apply \"Font Fix\", of Font et al., ONLY when the primitives (con2prim) solver\n", "< * fails. Note that this version is OPTIMIZED for gamma=2 polytrope EOS only.\n", "< * Application of this fix GUARANTEES (on an analytic level, at least) \n", "< * that con2prim will succeed.\n", "< *\n", "< * FIXME: There exist better fixes/strategies now. See McKinney et al's work.\n", "< *\n", "< -----------------------------------------------------------------------------*/\n", "12,39d1\n", "< inline int font_fix_gamma_equals2(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,eos_struct &eos) {\n", "< CCTK_REAL rhob;\n", "< CCTK_REAL kpoly2 = 2.0*eos.K_poly;\n", "< CCTK_REAL tol = 1.e-15;\n", "< CCTK_REAL Bxbar = PRIMS[BX_CENTER]*ONE_OVER_SQRT_4PI;\n", "< CCTK_REAL Bybar = PRIMS[BY_CENTER]*ONE_OVER_SQRT_4PI;\n", "< CCTK_REAL Bzbar = PRIMS[BZ_CENTER]*ONE_OVER_SQRT_4PI;\n", "< CCTK_REAL Bbar_x = METRIC_PHYS[GXX]*Bxbar + METRIC_PHYS[GXY]*Bybar + METRIC_PHYS[GXZ]*Bzbar;\n", "< CCTK_REAL Bbar_y = METRIC_PHYS[GXY]*Bxbar + METRIC_PHYS[GYY]*Bybar + METRIC_PHYS[GYZ]*Bzbar;\n", "< CCTK_REAL Bbar_z = METRIC_PHYS[GXZ]*Bxbar + METRIC_PHYS[GYZ]*Bybar + METRIC_PHYS[GZZ]*Bzbar;\n", "< CCTK_REAL B2bar = Bxbar*Bbar_x + Bybar*Bbar_y + Bzbar*Bbar_z;\n", "< CCTK_REAL Bbar = sqrt(B2bar);\n", "< CCTK_REAL check_B_small = fabs(Bxbar)+fabs(Bybar)+fabs(Bzbar);\n", "< if (check_B_small>0 && check_B_small<1.e-150) {\n", "< // need to compute B2bar specially to prevent floating-point underflow\n", "< CCTK_REAL Bmax = fabs(Bxbar);\n", "< if (Bmax < fabs(Bybar)) Bmax=fabs(Bybar);\n", "< if (Bmax < fabs(Bzbar)) Bmax=fabs(Bzbar);\n", "< CCTK_REAL Bxtmp=Bxbar/Bmax, Bytemp=Bybar/Bmax, Bztemp=Bzbar/Bmax;\n", "< CCTK_REAL B_xtemp=Bbar_x/Bmax, B_ytemp=Bbar_y/Bmax, B_ztemp=Bbar_z/Bmax;\n", "< Bbar = sqrt(Bxtmp*B_xtemp + Bytemp*B_ytemp + Bztemp*B_ztemp)*Bmax;\n", "< }\n", "< \n", "< CCTK_REAL BbardotS = Bxbar*CONSERVS[STILDEX] + Bybar*CONSERVS[STILDEY] + Bzbar*CONSERVS[STILDEZ];\n", "< CCTK_REAL BbardotS2 = BbardotS*BbardotS;\n", "< CCTK_REAL hatBbardotS = BbardotS/Bbar;\n", "< if (Bbar<1.e-300) hatBbardotS = 0.0;\n", "< CCTK_REAL Psim6 = 1.0/METRIC_LAP_PSI4[PSI6];\n", "41,165c3,7\n", "< // Limit hatBbardotS\n", "< //CCTK_REAL max_gammav = 100.0;\n", "< //CCTK_REAL rhob_max = CONSERVS[RHOSTAR]*Psim6;\n", "< //CCTK_REAL hmax = 1.0 + kpoly2*rhob_max;\n", "< //CCTK_REAL abs_hatBbardotS_max = sqrt(SQR(max_gammav)-1.0)*CONSERVS[RHOSTAR]*hmax;\n", "< //if (fabs(hatBbardotS) > abs_hatBbardotS_max) {\n", "< // CCTK_REAL fac_reduce = abs_hatBbardotS_max/fabs(hatBbardotS);\n", "< // CCTK_REAL hatBbardotS_max = hatBbardotS*fac_reduce;\n", "< // CCTK_REAL Bbar_inv = 1.0/Bbar;\n", "< // CCTK_REAL hat_Bbar_x = Bbar_x*Bbar_inv;\n", "< // CCTK_REAL hat_Bbar_y = Bbar_y*Bbar_inv;\n", "< // CCTK_REAL hat_Bbar_z = Bbar_z*Bbar_inv;\n", "< // CCTK_REAL sub_fact = hatBbardotS_max - hatBbardotS;\n", "< // CONSERVS[STILDEX] += sub_fact*hat_Bbar_x;\n", "< // CONSERVS[STILDEY] += sub_fact*hat_Bbar_y;\n", "< // CONSERVS[STILDEZ] += sub_fact*hat_Bbar_z;\n", "< // hatBbardotS = hatBbardotS_max;\n", "< // BbardotS *= fac_reduce;\n", "< // BbardotS2 = BbardotS*BbardotS;\n", "< //}\n", "< \n", "< CCTK_REAL sdots = METRIC_PHYS[GUPXX]*SQR(CONSERVS[STILDEX]) + METRIC_PHYS[GUPYY]*SQR(CONSERVS[STILDEY]) + METRIC_PHYS[GUPZZ]*SQR(CONSERVS[STILDEZ])\n", "< + 2.0*( METRIC_PHYS[GUPXY]*CONSERVS[STILDEX]*CONSERVS[STILDEY] + METRIC_PHYS[GUPXZ]*CONSERVS[STILDEX]*CONSERVS[STILDEZ]\n", "< + METRIC_PHYS[GUPYZ]*CONSERVS[STILDEY]*CONSERVS[STILDEZ]);\n", "< \n", "< if (sdots<1.e-300) {\n", "< rhob = CONSERVS[RHOSTAR]*Psim6;\n", "< u_x=0.0; u_y=0.0; u_z=0.0;\n", "< return 0;\n", "< }\n", "< \n", "< /*\n", "< if (fabs(BbardotS2 - sdots*B2bar) > 1e-8) {\n", "< CCTK_VInfo(CCTK_THORNSTRING,\"(Bbar dot S)^2, Bbar^2 * sdotS, %e %e\",SQR(BbardotS),sdots*B2bar);\n", "< CCTK_VInfo(CCTK_THORNSTRING,\"Cauchy-Schwartz inequality is violated!\");\n", "< }\n", "< */\n", "< \n", "< // Initial guess for W, S_fluid and rhob\n", "< CCTK_REAL W0 = sqrt( SQR(hatBbardotS) + SQR(CONSERVS[RHOSTAR]) ) * Psim6;\n", "< CCTK_REAL Sf20 = (SQR(W0)*sdots + BbardotS2*(B2bar + 2.0*W0))/SQR(W0+B2bar);\n", "< CCTK_REAL rhob0 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]));\n", "< CCTK_REAL W=W0,Sf2=Sf20,rhob1=rhob0;\n", "< \n", "< //****************************************************************\n", "< // FONT FIX\n", "< // Impose Font fix when HARM primitives solver fails to find\n", "< // acceptable set of primitives.\n", "< //****************************************************************\n", "< bool fontcheck=true;\n", "< \n", "< int itcount = 0, maxits=300;\n", "< while(fontcheck && itcount < maxits) {\n", "< itcount++;\n", "< W0 = W;\n", "< Sf20 = Sf2;\n", "< rhob0 = rhob1;\n", "< \n", "< // first, find rhob for the given S_fluid^2\n", "< CCTK_REAL h = 1.0 + kpoly2*rhob0;\n", "< rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h));\n", "< while( fabs(rhob1-rhob0) > rhob1*tol) {\n", "< rhob0 = rhob1;\n", "< h = 1.0 + kpoly2*rhob0;\n", "< rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h));\n", "< }\n", "< \n", "< h = 1.0 + kpoly2*rhob1;\n", "< W = sqrt( Sf20 + SQR(CONSERVS[RHOSTAR]*h))*Psim6;\n", "< Sf2 = (SQR(W)*sdots + BbardotS2*(B2bar + 2.0*W))/SQR(W+B2bar);\n", "< if ( fabs(W-W0) < W*tol && fabs(Sf20-Sf2) < Sf2*tol) fontcheck=false;\n", "< }\n", "< \n", "< if (itcount>=maxits) {\n", "< // Increase tol and try again \n", "< maxits*=100;\n", "< tol *=10.0;\n", "< itcount = 0;\n", "< fontcheck=true;\n", "< while(fontcheck && itcount < maxits) {\n", "< itcount++;\n", "< W0 = W;\n", "< Sf20 = Sf2;\n", "< rhob0 = rhob1;\n", "< \n", "< // first, find rhob for the given S_fluid^2\n", "< CCTK_REAL h = 1.0 + kpoly2*rhob0; \n", "< rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h));\n", "< while( fabs(rhob1-rhob0) > rhob1*tol) {\n", "< rhob0 = rhob1;\n", "< h = 1.0 + kpoly2*rhob0;\n", "< rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h));\n", "< }\n", "< \n", "< h = 1.0 + kpoly2*rhob1;\n", "< W = sqrt( Sf20 + SQR(CONSERVS[RHOSTAR]*h))*Psim6;\n", "< Sf2 = (SQR(W)*sdots + BbardotS2*(B2bar + 2.0*W))/SQR(W+B2bar);\n", "< if ( fabs(W-W0) < W*tol && fabs(Sf20-Sf2) < Sf2*tol) fontcheck=false;\n", "< }\n", "< }\n", "< //************************************************************************************************************** \n", "< \n", "< if(fontcheck==true) {\n", "< return 1;\n", "< }\n", "< \n", "< // Font fix works, now compute u_i\n", "< rhob = rhob1;\n", "< CCTK_REAL h = 1.0+kpoly2*rhob;\n", "< CCTK_REAL gammav = CONSERVS[RHOSTAR]*Psim6/rhob;\n", "< CCTK_REAL rhosh = CONSERVS[RHOSTAR]*h;\n", "< CCTK_REAL fac1 = METRIC_LAP_PSI4[PSI6]*BbardotS/(gammav*rhosh);\n", "< CCTK_REAL fac2 = 1.0/(rhosh + METRIC_LAP_PSI4[PSI6]*B2bar/gammav);\n", "< u_x = fac2*(CONSERVS[STILDEX] + fac1*Bbar_x);\n", "< u_y = fac2*(CONSERVS[STILDEY] + fac1*Bbar_y);\n", "< u_z = fac2*(CONSERVS[STILDEZ] + fac1*Bbar_z);\n", "< return 0;\n", "< }\n", "< \n", "< inline int font_fix_general_gamma(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,eos_struct &eos) {\n", "< CCTK_REAL gamma=eos.gamma_tab[0];\n", "< CCTK_REAL rhob;\n", "< CCTK_REAL tol = 1.e-15;\n", "< CCTK_REAL gam1 = gamma-1.0;\n", "< CCTK_REAL gam_gamm1_kpoly = gamma/gam1*eos.K_poly;\n", "---\n", "> /**********************************\n", "> * Piecewise Polytropic EOS Patch *\n", "> * Font fix: function call *\n", "> **********************************/\n", "> inline int font_fix__hybrid_EOS(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4, eos_struct eos) {\n", "216a59,60\n", "> \n", "> CCTK_REAL rhob;\n", "227a72,73\n", "> \n", "> \n", "229,230c75,76\n", "< CCTK_REAL W0 = sqrt( SQR(hatBbardotS) + SQR(CONSERVS[RHOSTAR]) ) * Psim6;\n", "< CCTK_REAL Sf20 = (SQR(W0)*sdots + BbardotS2*(B2bar + 2.0*W0))/SQR(W0+B2bar);\n", "---\n", "> CCTK_REAL W0 = sqrt( SQR(hatBbardotS) + SQR(CONSERVS[RHOSTAR]) ) * Psim6;\n", "> CCTK_REAL Sf20 = (SQR(W0)*sdots + BbardotS2*(B2bar + 2.0*W0))/SQR(W0+B2bar);\n", "232c78\n", "< CCTK_REAL W=W0,Sf2=Sf20,rhob1=rhob0;\n", "---\n", "> \n", "240c86,87\n", "< bool fontcheck=true;\n", "---\n", "> /* Set the maximum number of iterations */\n", "> int maxits = 500;\n", "242,262c89,90\n", "< int itcount = 0, maxits=500;\n", "< while(fontcheck && itcount < maxits) {\n", "< itcount++;\n", "< W0 = W;\n", "< Sf20 = Sf2;\n", "< rhob0 = rhob1;\n", "< \n", "< // first, find rhob for the given S_fluid^2\n", "< CCTK_REAL h = 1.0 + gam_gamm1_kpoly*pow(rhob0,gam1); \n", "< rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h));\n", "< while( fabs(rhob1-rhob0) > rhob1*tol) {\n", "< rhob0 = rhob1;\n", "< h = 1.0 + gam_gamm1_kpoly*pow(rhob0,gam1);\n", "< rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h));\n", "< }\n", "< \n", "< h = 1.0 + gam_gamm1_kpoly*pow(rhob1,gam1);\n", "< W = sqrt( Sf20 + SQR(CONSERVS[RHOSTAR]*h))*Psim6;\n", "< Sf2 = (SQR(W)*sdots + BbardotS2*(B2bar + 2.0*W))/SQR(W+B2bar);\n", "< if ( fabs(W-W0) < W*tol && fabs(Sf20-Sf2) < Sf2*tol) fontcheck=false;\n", "< }\n", "---\n", "> /* Set the allowed tolerance */\n", "> CCTK_REAL tol = 1.e-15;\n", "264,290c92,110\n", "< if (itcount>=maxits) {\n", "< // Increase tol and try again\n", "< fontcheck=true;\n", "< tol *=100.0;\n", "< itcount = 0;\n", "< while(fontcheck && itcount < maxits) {\n", "< itcount++;\n", "< W0 = W;\n", "< Sf20 = Sf2;\n", "< rhob0 = rhob1;\n", "< \n", "< // first, find rhob for the given S_fluid^2\n", "< CCTK_REAL h = 1.0 + gam_gamm1_kpoly*pow(rhob0,gam1);\n", "< rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h));\n", "< while( fabs(rhob1-rhob0) > rhob1*tol) {\n", "< rhob0 = rhob1;\n", "< h = 1.0 + gam_gamm1_kpoly*pow(rhob0,gam1);\n", "< rhob1 = CONSERVS[RHOSTAR]*Psim6/sqrt(1.0+Sf20/SQR(CONSERVS[RHOSTAR]*h));\n", "< }\n", "< \n", "< h = 1.0 + gam_gamm1_kpoly*pow(rhob1,gam1);\n", "< W = sqrt( Sf20 + SQR(CONSERVS[RHOSTAR]*h))*Psim6;\n", "< Sf2 = (SQR(W)*sdots + BbardotS2*(B2bar + 2.0*W))/SQR(W+B2bar);\n", "< if ( fabs(W-W0) < W*tol && fabs(Sf20-Sf2) < Sf2*tol) fontcheck=false;\n", "< }\n", "< }\n", "< //************************************************************************************************************** \n", "---\n", "> /* Declare basic variables */\n", "> int font_fix_status;\n", "> \n", "> /**********************\n", "> * FONT FIX MAIN LOOP *\n", "> **********************\n", "> * Perform the font fix routine until convergence\n", "> * is obtained and the algorithm returns with no\n", "> * error. Every time the Font fix fails, increase\n", "> * the tolerance by a factor of 10.\n", "> */\n", "> int font_fix_attempts = 5;\n", "> CCTK_REAL font_fix_tol_factor = 10.0;\n", "> for(int n=0; n \n", "> tol *= pow(font_fix_tol_factor,n);\n", "> font_fix_status = font_fix__rhob_loop(maxits,tol, W0,Sf20,Psim6,sdots,BbardotS2,B2bar, CONSERVS,eos, rhob0,rhob);\n", "> rhob0 = rhob;\n", "> if(font_fix_status==0) break;\n", "292,293d111\n", "< if(fontcheck==true) {\n", "< return 1;\n", "296,298c114,128\n", "< // Font fix works, now compute u_i\n", "< rhob = rhob1;\n", "< CCTK_REAL h = 1.0+gam_gamm1_kpoly*pow(rhob,gam1);\n", "---\n", "> \n", "> //**************************************************************************************************************\n", "> \n", "> /* Font fix works! */\n", "> /* First compute P_cold, eps_cold, then h = h_cold */\n", "> CCTK_REAL P_cold, eps_cold;\n", "> compute_P_cold__eps_cold(eos,rhob, P_cold,eps_cold);\n", "> CCTK_REAL h = 1.0 + eps_cold + P_cold/rhob;\n", "> \n", "> /* Then compute gamma_v using equation (A19) in\n", "> * Etienne et al. (2011) [https://arxiv.org/pdf/1112.0568.pdf]\n", "> * .-----------------------------------------.\n", "> * | gamma_v = psi^{-6} * (rho_star / rho_b) |\n", "> * .-----------------------------------------.\n", "> */\n", "299a130,131\n", "> \n", "> /* Finally, compute u_{i} */\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/font_fix_gamma_law.C\"\n", "original_IGM_file_name = \"font_fix_gamma_law-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__font_fix_gamma_law__C = !diff $original_IGM_file_path $outfile_path__font_fix_hybrid_EOS__C\n", "\n", "if Validation__font_fix_gamma_law__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 font_fix_gamma_law.C: PASSED!\")\n", "else:\n", " # If the validation fails, we keep the original IGM source code file\n", " print(\"Validation test for font_fix_gamma_law.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__font_fix_gamma_law__C:\n", " print(diff_line)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 6.d: `harm_primitives_headers.h` \\[Back to [top](#toc)\\]\n", "$$\\label{harm_primitives_headers_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": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Validation test for harm_primitives_headers.h: FAILED!\n", "Diff:\n", "2c2\n", "< Copyright 2006 Charles F. Gammie, Jonathan C. McKinney, Scott C. Noble, \n", "---\n", "> Copyright 2006 Charles F. Gammie, Jonathan C. McKinney, Scott C. Noble,\n", "7c7\n", "< This file is part of HARM. HARM is a program that solves hyperbolic \n", "---\n", "> This file is part of HARM. HARM is a program that solves hyperbolic\n", "9,10c9,10\n", "< shock-capturing techniques. This version of HARM has been configured to \n", "< solve the relativistic magnetohydrodynamic equations of motion on a \n", "---\n", "> shock-capturing techniques. This version of HARM has been configured to\n", "> solve the relativistic magnetohydrodynamic equations of motion on a\n", "12c12\n", "< an accretion disk model. \n", "---\n", "> an accretion disk model.\n", "14c14\n", "< You are morally obligated to cite the following two papers in his/her \n", "---\n", "> You are morally obligated to cite the following two papers in his/her\n", "17c17\n", "< [1] Gammie, C. F., McKinney, J. C., \\& Toth, G.\\ 2003, \n", "---\n", "> [1] Gammie, C. F., McKinney, J. C., \\& Toth, G.\\ 2003,\n", "20c20\n", "< [2] Noble, S. C., Gammie, C. F., McKinney, J. C., \\& Del Zanna, L. \\ 2006, \n", "---\n", "> [2] Noble, S. C., Gammie, C. F., McKinney, J. C., \\& Del Zanna, L. \\ 2006,\n", "23,24c23,24\n", "< \n", "< Further, we strongly encourage you to obtain the latest version of \n", "---\n", "> \n", "> Further, we strongly encourage you to obtain the latest version of\n", "59c59\n", "< #if( USE_ISENTROPIC ) \n", "---\n", "> #if( USE_ISENTROPIC )\n", "89c89\n", "< static const int RHO =0; \n", "---\n", "> static const int RHO =0;\n", "106c106\n", "< int Utoprim_2d(CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM], \n", "---\n", "> int Utoprim_2d(eos_struct eos, CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM],\n", "115,117c115\n", "< inline int font_fix_general_gamma(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,eos_struct &eos);\n", "< inline int font_fix_gamma_equals2(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4,eos_struct &eos);\n", "< \n", "---\n", "> inline int font_fix__hybrid_EOS(CCTK_REAL &u_x, CCTK_REAL &u_y, CCTK_REAL &u_z,CCTK_REAL *CONSERVS,CCTK_REAL *PRIMS,CCTK_REAL *METRIC_PHYS,CCTK_REAL *METRIC_LAP_PSI4, eos_struct eos);\n", "123a122\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/harm_primitives_headers.h\"\n", "original_IGM_file_name = \"harm_primitives_headers-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__harm_primitives_headers__h = !diff $original_IGM_file_path $outfile_path__harm_primitives_headers__h\n", "\n", "if Validation__harm_primitives_headers__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 harm_primitives_headers.h: PASSED!\")\n", "else:\n", " # If the validation fails, we keep the original IGM source code file\n", " print(\"Validation test for harm_primitives_headers.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__harm_primitives_headers__h:\n", " print(diff_line)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 7: 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__driver_conserv_to_prims.pdf](Tutorial-IllinoisGRMHD__driver_conserv_to_prims.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": 41, "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__the_conservative_to_primitive_algorithm.ipynb\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__the_conservative_to_primitive_algorithm.tex\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__the_conservative_to_primitive_algorithm.tex\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__the_conservative_to_primitive_algorithm.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.7.3" } }, "nbformat": 4, "nbformat_minor": 4 }