{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Tutorial-IllinoisGRMHD: harm_utoprim_2d.c\n", "\n", "## Authors: Leo Werneck & Zach Etienne\n", "\n", "**This module is currently under development**\n", "\n", "## This tutorial provides a comprehensive explaination of the conservative-to-primitive algorithm used by `HARM`. \n", "\n", "### Note: This module will likely be absorbed by another one once we finish documenting the code.\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](#harm_utoprim_2d__c__eos_indep): **EOS independent routines**\n", " 1. [Step 2.a](#utoprim_2d): *The `Utoprim_2d()` function*\n", " 1. [Step 2.a.i](#utoprim_2d__bi_and_alpha): Setting $B^{i}_{\\rm HARM}$ and $\\alpha$\n", " 1. [Step 2.a.ii](#utoprim_2d__converting): Preparing the variables to be used by the `Utoprim_new_body()` function\n", " 1. [Step 2.b](#utoprim_new_body): *The `Utoprim_new_body()` function*\n", " 1. [Step 2.b.i](#utoprim_new_body__basic_quantities): Computing basic quantities\n", " 1. [Step 2.b.ii](#utoprim_new_body__wlast): Determining $W$ from the previous iteration, $W_{\\rm last}$\n", " 1. [Step 2.b.iii](#utoprim_new_body__vsqlast_and_recompute_w_and_vsq): Compute $v^{2}_{\\rm last}$, then update $v^{2}$ and $W$\n", " 1. [Step 2.b.iv](#utoprim_new_body__compute_prims): Computing the primitive variables\n", " 1. [Step 2.c](#vsq_calc): *The `vsq_calc()` function*\n", " 1. [Step 2.d](#x1_of_x0): *The `x1_of_x0()` function*\n", " 1. [Step 2.e](#validate_x): *The `validate_x()` function*\n", " 1. [Step 2.f](#general_newton_raphson): *The `general_newton_raphson()` function*\n", " 1. [Step 2.g](#func_vsq): *The `func_vsq()` function*\n", "1. [Step 3](#harm_utoprim_2d__c__eos_dep): **EOS dependent routines**\n", " 1. [Step 3.a](#pressure_w_vsq): *The `pressure_W_vsq()` function*\n", " 1. [Step 3.b](#dpdw_calc_vsq): *The `dpdW_calc_vsq()` function*\n", " 1. [Step 3.c](#dpdvsq_calc): *The `dpdvsq_calc()` function*\n", " 1. [Step 3.c.i](#dpdvsq_calc__basic_quantities): Setting basic quantities and computing $P_{\\rm cold}$ and $\\epsilon_{\\rm cold}$\n", " 1. [Step 3.c.ii](#dpdvsq_calc__dpcolddvsq): Computing $\\frac{\\partial P_{\\rm cold}}{\\partial\\left(v^{2}\\right)}$\n", " 1. [Step 3.c.iii](#dpdvsq_calc__depscolddvsq): Computing $\\frac{\\partial \\epsilon_{\\rm cold}}{\\partial\\left(v^{2}\\right)}$\n", " 1. [Step 3.c.iv](#dpdvsq_calc__dpdvsq): Computing $\\frac{\\partial p_{\\rm hybrid}}{\\partial\\left(v^{2}\\right)}$\n", "1. [Step 4](#code_validation): **Code validation**\n", "1. [Step 5](#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__harm_utoprim_2d__c = os.path.join(IGM_src_dir_path,\"harm_utoprim_2d.c\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 1: Introduction \\[Back to [top](#toc)\\]\n", "$$\\label{introduction}$$\n", "\n", "Comment on license: `HARM` uses GPL, while `IllinoisGRMHD` uses BSD." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: EOS independent routines \\[Back to [top](#toc)\\]\n", "$$\\label{harm_utoprim_2d__c__eos_indep}$$\n", "\n", "Let us now start documenting the `harm_utoprim_2d.c`, which is a part of the `Harm` code. Our main reference throughout this discussion will be the required citation [Noble *et al.* (2006)](https://arxiv.org/abs/astro-ph/0512420).\n", "\n", "We will start with the code's required preamble." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Writing ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile $outfile_path__harm_utoprim_2d__c\n", "#ifndef __HARM_UTOPRIM_2D__C__\n", "#define __HARM_UTOPRIM_2D__C__\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", "\n", "/*************************************************************************************/\n", "/*************************************************************************************/\n", "/*************************************************************************************\n", "\n", "utoprim_2d.c:\n", "---------------\n", "\n", " Uses the 2D method:\n", " -- solves for two independent variables (W,v^2) via a 2D\n", " Newton-Raphson method\n", " -- can be used (in principle) with a general equation of state.\n", "\n", " -- Currently returns with an error state (>0) if a negative rest-mass\n", " density or internal energy density is calculated. You may want\n", " to change this aspect of the code so that it still calculates the\n", " velocity and so that you can floor the densities. If you want to\n", " change this aspect of the code please comment out the \"return(retval)\"\n", " statement after \"retval = 5;\" statement in Utoprim_new_body();\n", "\n", "******************************************************************************/\n", "\n", "static const int NEWT_DIM=2;\n", "\n", "// Declarations:\n", "static CCTK_REAL vsq_calc(CCTK_REAL W,CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D);\n", "static int Utoprim_new_body(eos_struct eos, CCTK_REAL U[], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL gdet, CCTK_REAL prim[],long &n_iter);\n", "static int general_newton_raphson( eos_struct eos, CCTK_REAL x[], int n, long &n_iter, void (*funcd) (eos_struct, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *, CCTK_REAL *, int,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &),CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D);\n", "static void func_vsq( eos_struct eos, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *f, CCTK_REAL *df, int n,CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D);\n", "static CCTK_REAL x1_of_x0(CCTK_REAL x0, CCTK_REAL &Bsq, CCTK_REAL &QdotBsq, CCTK_REAL &Qtsq, CCTK_REAL &Qdotn, CCTK_REAL &D ) ;\n", "static CCTK_REAL pressure_W_vsq(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) ;\n", "static CCTK_REAL dpdW_calc_vsq(CCTK_REAL W, CCTK_REAL vsq);\n", "static CCTK_REAL dpdvsq_calc(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D);\n", "\n", "/**********************************************************************/\n", "/******************************************************************\n", "\n", " Utoprim_2d():\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. Note that Greek\n", " indices run 0,1,2,3 and Latin indices run 1,2,3 (spatial only).\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", "\n", " Arguments:\n", " U[NPR] = conserved variables (current values on input/output);\n", " gcov[NDIM][NDIM] = covariant form of the metric ;\n", " gcon[NDIM][NDIM] = contravariant form of the metric ;\n", " gdet = sqrt( - determinant of the metric) ;\n", " prim[NPR] = primitive variables (guess on input, calculated values on\n", " output if there are no problems);\n", "\n", " -- NOTE: for those using this routine for special relativistic MHD and are\n", " unfamiliar with metrics, merely set\n", " gcov = gcon = diag(-1,1,1,1) and gdet = 1. ;\n", "\n", "******************************************************************/" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.a: The `Utoprim_2d()` function \\[Back to [top](#toc)\\]\n", "$$\\label{utoprim_2d}$$\n", "\n", "The `Utoprim_2d()` function is the driver function of the `HARM` conservative-to-primitive algorithm. We remind you from the definitions of primitive and conservative variables used in the code:\n", "\n", "$$\n", "\\begin{align}\n", "\\boldsymbol{P}_{\\rm HARM} &= \\left\\{\\rho_{b},u,\\tilde{u}^{i},B^{i}_{\\rm HARM}\\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\\}\\ .\n", "\\end{align}\n", "$$\n", "\n", "\n", "\n", "### Step 2.a.i: Setting $B^{i}_{\\rm HARM}$ and $\\alpha$ \\[Back to [top](#toc)\\]\n", "$$\\label{utoprim_2d__bi_and_alpha}$$\n", "\n", "Let\n", "\n", "$$\n", "\\tilde{B}^{i}_{\\rm HARM} \\equiv \\sqrt{-g}B^{i}_{\\rm HARM}\\ .\n", "$$\n", "\n", "The code starts by relating\n", "\n", "$$\n", "\\boxed{B^{i}_{\\rm HARM} = \\frac{\\tilde{B}^{i}_{\\rm HARM}}{\\sqrt{-g}}}\\ ,\n", "$$\n", "\n", "and setting\n", "\n", "$$\n", "\\boxed{\\alpha = \\frac{1}{\\sqrt{-g^{00}}}} \\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\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", "\n", " CCTK_REAL U_tmp[NPR], prim_tmp[NPR];\n", " int i, ret;\n", " CCTK_REAL alpha;\n", "\n", " if( U[0] <= 0. ) {\n", " return(-100);\n", " }\n", "\n", " /* First update the primitive B-fields */\n", " for(i = BCON1; i <= BCON3; i++) prim[i] = U[i] / gdet ;\n", "\n", " /* Set the geometry variables: */\n", " alpha = 1.0/sqrt(-gcon[0][0]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 2.a.ii: Preparing the variables to be used by the `Utoprim_new_body()` function \\[Back to [top](#toc)\\]\n", "$$\\label{utoprim_2d__converting}$$\n", "\n", "The conservative-to-primitive algorithm uses the `Utoprim_new_body()` function. However, this function assumes a *different* set of primitive/conservative variables. Thus, we must perform the proper conversion. First, let us ease on the notation a bit by defining:\n", "\n", "$$\n", "\\boldsymbol{C} \\equiv \\left\\{\\rho_{\\star},u_{\\star},\\tilde{S}_{i},\\tilde{B}^{i}_{\\rm HARM}\\right\\} \\equiv \\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\\}\\ .\n", "$$\n", "\n", "\n", "\n", "Below we list the main differences in the conservative variables:\n", "\n", "| `Utoprim_2d()` | `Utoprim_new_body()` |\n", "|------------------------------------------|---------------------------------------------------------------------------|\n", "| $\\color{blue}{\\textbf{Conservatives}}$ | $\\color{red}{\\textbf{Conservatives}}$ |\n", "| $\\color{blue}{\\rho_{\\star}}$ | $\\color{red}{\\frac{\\alpha}{\\sqrt{-g}}\\rho_{\\star}}$ |\n", "| $\\color{blue}{u_{\\star}}$ | $\\color{red}{\\frac{\\alpha}{\\sqrt{-g}}\\left(u_{\\star}-\\rho_{\\star}\\right)}$|\n", "| $\\color{blue}{\\tilde{S}_{i}}$ | $\\color{red}{\\frac{\\alpha}{\\sqrt{-g}}\\tilde{S}_{i}}$ |\n", "| $\\color{blue}{\\tilde{B}^{i}_{\\rm HARM}}$ | $\\color{red}{\\frac{\\alpha}{\\sqrt{-g}}\\tilde{B}^{i}_{\\rm HARM}}$ |\n", "\n", "These are necessary conversions because while `Utoprim_2d()` assumes the set of conservatives above, `Utoprim_new_body()` assumes\n", "\n", "$$\n", "\\left\\{\\gamma\\rho_{b},\\alpha T^{0}_{\\ \\ 0}, \\alpha T^{0}_{\\ \\ i}, \\alpha B^{i}_{\\rm HARM}\\right\\}\\ .\n", "$$\n", "\n", "Let us first pause to understand the table above. From definition (15) in [Noble *et al.* (2006)](https://arxiv.org/abs/astro-ph/0512420) and the discussion just below it, we know that $\\gamma = \\alpha u^{0}$. Thus\n", "\n", "$$\n", "\\rho_{\\star} = \\sqrt{-g}\\rho_{b}u^{0} = \\sqrt{-g}\\left(\\frac{\\gamma}{\\alpha}\\rho_{b}\\right)\\implies\\boxed{\\gamma \\rho_{b} = \\frac{\\alpha}{\\sqrt{-g}}\\rho_{\\star}}\\ .\n", "$$\n", "\n", "Then we have\n", "\n", "$$\n", "u_{\\star} = \\sqrt{-g}\\left(T^{0}_{\\ \\ 0} + \\rho_{b}u^{0}\\right)= \\sqrt{-g}\\left(T^{0}_{\\ \\ 0} + \\frac{\\rho_{\\star}}{\\sqrt{-g}}\\right) = \\sqrt{-g}T^{0}_{\\ \\ 0} + \\rho_{\\star} \\implies \\boxed{\\alpha T^{0}_{\\ \\ 0} = \\frac{\\alpha}{\\sqrt{-g}}\\left(u_{\\star}-\\rho_{\\star}\\right)}\\ .\n", "$$\n", "\n", "The other two relations are more straightforward. We have\n", "\n", "$$\n", "\\tilde{S}_{i} = \\sqrt{-g}T^{0}_{\\ \\ i} \\implies \\boxed{\\alpha T^{0}_{\\ \\ i} = \\frac{\\alpha}{\\sqrt{-g}}\\tilde{S}_{i}}\\ ,\n", "$$\n", "\n", "and\n", "\n", "$$\n", "\\tilde{B}^{i}_{\\rm HARM} = \\sqrt{-g}B^{i}_{\\rm HARM}\\implies \\boxed{\\alpha B^{i}_{\\rm HARM} = \\frac{\\alpha}{\\sqrt{-g}}\\tilde{B}^{i}_{\\rm HARM}}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", " /* Transform the CONSERVED variables into the new system */\n", " U_tmp[RHO] = alpha * U[RHO] / gdet;\n", " U_tmp[UU] = alpha * (U[UU] - U[RHO]) / gdet ;\n", " for( i = UTCON1; i <= UTCON3; i++ ) {\n", " U_tmp[i] = alpha * U[i] / gdet ;\n", " }\n", " for( i = BCON1; i <= BCON3; i++ ) {\n", " U_tmp[i] = alpha * U[i] / gdet ;\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Below we list the necessary transformations on the primitive variables:\n", "\n", "| `Utoprim_2d()` | `Utoprim_new_body()` |\n", "|-------------------------------------|----------------------------------------|\n", "| $\\color{blue}{\\textbf{Primitives}}$ | $\\color{red}{\\textbf{Primitives}}$ |\n", "| $\\color{blue}{\\rho_{b}}$ | $\\color{red}{\\rho_{b}}$ |\n", "| $\\color{blue}{u}$ | $\\color{red}{u}$ |\n", "| $\\color{blue}{\\tilde{u}^{i}}$ | $\\color{red}{\\tilde{u}^{i}}$ |\n", "| $\\color{blue}{B^{i}_{\\rm HARM}}$ | $\\color{red}{\\alpha B^{i}_{\\rm HARM}}$ |\n", "\n", "After this slight modification, we call the `Utoprim_new_body()` function. If it returns without errors, then the variables ${\\rm prim\\_tmp}$ will now contain the values of the primitives. We then update the ${\\rm prim}$ variables with these newly computed values." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", " /* Transform the PRIMITIVE variables into the new system */\n", " for( i = 0; i < BCON1; i++ ) {\n", " prim_tmp[i] = prim[i];\n", " }\n", " for( i = BCON1; i <= BCON3; i++ ) {\n", " prim_tmp[i] = alpha*prim[i];\n", " }\n", "\n", " ret = Utoprim_new_body(eos, U_tmp, gcov, gcon, gdet, prim_tmp,n_iter);\n", "\n", " /* Transform new primitive variables back if there was no problem : */\n", " if( ret == 0 || ret == 5 || ret==101 ) {\n", " for( i = 0; i < BCON1; i++ ) {\n", " prim[i] = prim_tmp[i];\n", " }\n", " }\n", "\n", " return( ret ) ;\n", "\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.b: The `Utoprim_new_body()` function \\[Back to [top](#toc)\\]\n", "$$\\label{utoprim_new_body}$$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", "\n", "/**********************************************************************/\n", "/**********************************************************************************\n", "\n", " Utoprim_new_body():\n", "\n", " -- Attempt an inversion from U to prim using the initial guess prim.\n", "\n", " -- This is the main routine that calculates auxiliary quantities for the\n", " Newton-Raphson routine.\n", "\n", " -- assumes that\n", " / rho gamma \\\n", " U = | alpha T^t_\\mu |\n", " \\ alpha B^i /\n", "\n", "\n", "\n", " / rho \\\n", " prim = | uu |\n", " | \\tilde{u}^i |\n", " \\ alpha B^i /\n", "\n", "\n", "return: (i*100 + j) where\n", " i = 0 -> Newton-Raphson solver either was not called (yet or not used)\n", " or returned successfully;\n", " 1 -> Newton-Raphson solver did not converge to a solution with the\n", " given tolerances;\n", " 2 -> Newton-Raphson procedure encountered a numerical divergence\n", " (occurrence of \"nan\" or \"+/-inf\" ;\n", "\n", " j = 0 -> success\n", " 1 -> failure: some sort of failure in Newton-Raphson;\n", " 2 -> failure: utsq<0 w/ initial p[] guess;\n", " 3 -> failure: W<0 or W>W_TOO_BIG\n", " 4 -> failure: v^2 > 1\n", " 5 -> failure: rho,uu <= 0 ;\n", "\n", "**********************************************************************************/\n", "\n", "static int Utoprim_new_body(eos_struct eos, CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM],\n", " CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL gdet, CCTK_REAL prim[NPR], long &n_iter)\n", "{\n", "\n", " CCTK_REAL x_2d[NEWT_DIM];\n", " CCTK_REAL QdotB,Bcon[NDIM],Bcov[NDIM],Qcov[NDIM],Qcon[NDIM],ncov[NDIM],ncon[NDIM],Qsq,Qtcon[NDIM];\n", " CCTK_REAL rho0,u,p,w,gammasq,gamma,gtmp,W_last,W,utsq,vsq;\n", " int i,j, n, retval, i_increase;\n", "\n", " n = NEWT_DIM ;\n", "\n", " // Assume ok initially:\n", " retval = 0;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.b.i: Computing basic quantities \\[Back to [top](#toc)\\]\n", "$$\\label{utoprim_new_body__basic_quantities}$$\n", "\n", "We start by computing basic quantities from the input variables. Notice that this conservative-to-primitive algorithm does not need to update the magnetic field, thus\n", "\n", "$$\n", "\\boxed{B_{\\rm prim}^{i} = B_{\\rm conserv}^{i}}\\ .\n", "$$\n", "\n", "Since they are both equal, we will not distinguish between prim and conserve in what follows. We also set $B^{0} = 0$. Then we define\n", "\n", "$$\n", "\\boxed{Q_{\\mu} \\equiv \\alpha T^{0}_{\\ \\ \\mu}}\\ .\n", "$$\n", "\n", "From these, the following quantities are then computed:\n", "\n", "$$\n", "\\boxed{\n", "\\begin{align}\n", "B_{i} &= g_{i\\mu}B^{\\mu}\\\\\n", "Q^{\\mu} &= g^{\\mu\\nu}Q_{\\nu}\\\\\n", "B^{2} &= B_{i}B^{i}\\\\\n", "Q\\cdot B &= Q_{\\mu}B^{\\mu}\\\\\n", "\\left(Q\\cdot B\\right)^{2} &= \\left(Q\\cdot B\\right)\\left(Q\\cdot B\\right)\\\\\n", "n_{\\mu} &= \\left(-\\alpha,0,0,0\\right)\\\\\n", "n^{\\mu} &= g^{\\mu\\nu}n_{\\nu}\\\\\n", "\\left(Q\\cdot n\\right) &= Q^{\\mu}n_{\\mu}\\\\\n", "Q^{2} &= Q_{\\mu}Q^{\\mu}\\\\\n", "\\tilde{Q}^{2} &= Q^{2} + \\left(Q\\cdot n\\right)\\left(Q\\cdot n\\right)\\\\\n", "D &\\equiv \\gamma \\rho_{b}\n", "\\end{align}\n", "}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", " for(i = BCON1; i <= BCON3; i++) prim[i] = U[i] ;\n", "\n", " // Calculate various scalars (Q.B, Q^2, etc) from the conserved variables:\n", " Bcon[0] = 0. ;\n", " for(i=1;i<4;i++) Bcon[i] = U[BCON1+i-1] ;\n", "\n", " lower_g(Bcon,gcov,Bcov) ;\n", "\n", " for(i=0;i<4;i++) Qcov[i] = U[QCOV0+i] ;\n", " raise_g(Qcov,gcon,Qcon) ;\n", "\n", "\n", " CCTK_REAL Bsq = 0. ;\n", " for(i=1;i<4;i++) Bsq += Bcon[i]*Bcov[i] ;\n", "\n", " QdotB = 0. ;\n", " for(i=0;i<4;i++) QdotB += Qcov[i]*Bcon[i] ;\n", " CCTK_REAL QdotBsq = QdotB*QdotB ;\n", "\n", " ncov_calc(gcon,ncov) ;\n", " // FIXME: The exact form of n^{\\mu} can be found\n", " // in eq. (2.116) and implementing it\n", " // directly is a lot more efficient than\n", " // performing n^{\\mu} = g^{\\mu\\nu}n_{nu}\n", " raise_g(ncov,gcon,ncon);\n", "\n", " CCTK_REAL Qdotn = Qcon[0]*ncov[0] ;\n", "\n", " Qsq = 0. ;\n", " for(i=0;i<4;i++) Qsq += Qcov[i]*Qcon[i] ;\n", "\n", " CCTK_REAL Qtsq = Qsq + Qdotn*Qdotn ;\n", "\n", " CCTK_REAL D = U[RHO] ;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.b.ii: Determining $W$ from the previous iteration, $W_{\\rm last}$ \\[Back to [top](#toc)\\]\n", "$$\\label{utoprim_new_body__wlast}$$\n", "\n", "The quantity $W$ is defined as\n", "\n", "$$\n", "W \\equiv w\\gamma^{2}\\ ,\n", "$$\n", "\n", "where\n", "\n", "$$\n", "\\begin{align}\n", "w &= \\rho_{b} + u + p\\ ,\\\\\n", "\\gamma^{2} &= 1 + g_{ij}\\tilde{u}^{i}\\tilde{u}^{j}\\ .\n", "\\end{align}\n", "$$\n", "\n", "Thus the quantities $g_{ij}\\tilde{u}^{i}\\tilde{u}^{j}$ and then $\\gamma^{2}$ and $\\gamma$. Thus, by computing $\\rho_{b}$ and $p$ from the input variables, i.e. $D$, one can determine $w$ and then compute the value of $W$ from the input values (previous iteration), which we denote by $W_{\\rm last}$.\n", "\n", "**Dependecy note:** Note that this function depends on the `pressure_rho0_u()` function, which is *not* EOS independent." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", " /* calculate W from last timestep and use for guess */\n", " utsq = 0. ;\n", " for(i=1;i<4;i++)\n", " for(j=1;j<4;j++) utsq += gcov[i][j]*prim[UTCON1+i-1]*prim[UTCON1+j-1] ;\n", "\n", "\n", " if( (utsq < 0.) && (fabs(utsq) < 1.0e-13) ) {\n", " utsq = fabs(utsq);\n", " }\n", " if(utsq < 0. || utsq > UTSQ_TOO_BIG) {\n", " retval = 2;\n", " return(retval) ;\n", " }\n", "\n", " gammasq = 1. + utsq ;\n", " gamma = sqrt(gammasq);\n", "\n", " // Always calculate rho from D and gamma so that using D in EOS remains consistent\n", " // i.e. you don't get positive values for dP/d(vsq) .\n", " rho0 = D / gamma ;\n", " u = prim[UU] ;\n", " p = pressure_rho0_u(eos, rho0,u) ;\n", " w = rho0 + u + p ;\n", "\n", " W_last = w*gammasq ;\n", "\n", "\n", " // Make sure that W is large enough so that v^2 < 1 :\n", " i_increase = 0;\n", " while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq )\n", " - QdotBsq*(2.*W_last + Bsq) ) <= W_last*W_last*(Qtsq-Bsq*Bsq))\n", " && (i_increase < 10) ) {\n", " W_last *= 10.;\n", " i_increase++;\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.b.iii: Compute $v^{2}_{\\rm last}$, then update $v^{2}$ and $W$ \\[Back to [top](#toc)\\]\n", "$$\\label{utoprim_new_body__vsqlast_and_recompute_w_and_vsq}$$\n", "\n", "Then we use equation (28) in [Noble *et al.* (2006)](https://arxiv.org/abs/astro-ph/0512420) to determine $v^{2}$:\n", "\n", "$$\n", "\\boxed{v^{2} = \\frac{\\tilde{Q}^{2}W^{2} + \\left(Q\\cdot B\\right)^{2}\\left(B^{2}+2W\\right)}{\\left(B^{2}+W\\right)^{2}W^{2}}}\\ .\n", "$$\n", "\n", "This is done by calling the `x1_of_x0()` function, where $x_{0} = W$ and $x_{1} = v^{2}$, which itself calls the `vsq_calc()` function which implements the boxed equation above.\n", "\n", "After we have $\\left\\{W_{\\rm last},v^{2}_{\\rm last}\\right\\}$ we use them as the initial guess for the `general_newton_raphson()`, which returns the updated values $\\left\\{W,v^{2}\\right\\}$.\n", "\n", "All functions mentioned above are documented in this tutorial notebook, so look at the [Table of Contents](#toc) for more information." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", " // Calculate W and vsq:\n", " x_2d[0] = fabs( W_last );\n", " x_2d[1] = x1_of_x0( W_last , Bsq,QdotBsq,Qtsq,Qdotn,D) ;\n", " retval = general_newton_raphson( eos, x_2d, n, n_iter, func_vsq, Bsq,QdotBsq,Qtsq,Qdotn,D) ;\n", "\n", " W = x_2d[0];\n", " vsq = x_2d[1];\n", "\n", " /* Problem with solver, so return denoting error before doing anything further */\n", " if( (retval != 0) || (W == FAIL_VAL) ) {\n", " retval = retval*100+1;\n", " return(retval);\n", " }\n", " else{\n", " if(W <= 0. || W > W_TOO_BIG) {\n", " retval = 3;\n", " return(retval) ;\n", " }\n", " }\n", "\n", " // Calculate v^2:\n", " if( vsq >= 1. ) {\n", " vsq = 1.-2.e-16;\n", " //retval = 4;\n", " //return(retval) ;\n", " }" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.b.iv: Computing the primitive variables \\[Back to [top](#toc)\\]\n", "$$\\label{utoprim_new_body__compute_prims}$$\n", "\n", "Now that we have $\\left\\{W,v^{2}\\right\\}$, we recompute the primitive variables. We start with\n", "\n", "$$\n", "\\left\\{\n", "\\begin{align}\n", "\\tilde{g} &\\equiv \\sqrt{1-v^{2}}\\\\\n", "\\gamma &= \\frac{1}{\\tilde{g}}\n", "\\end{align}\n", "\\right.\n", "\\implies\n", "\\boxed{\\rho_{b} = D\\tilde{g}}\\ .\n", "$$\n", "\n", "Then, we determine the pressure $p$ using the `pressure_rho0_w()` function and\n", "\n", "$$\n", "w = W\\left(1-v^{2}\\right)\n", "\\implies\n", "\\boxed{u = w - \\left(\\rho_{b} + p\\right)}\\ .\n", "$$\n", "\n", "**Dependecy note:** Note that this function depends on the `pressure_rho0_w()` function, which is *not* EOS independent.\n", "\n", "Finally, we can obtain $\\tilde{u}^{i}$ using eq. 31 in [Noble *et al.* (2006)](https://arxiv.org/abs/astro-ph/0512420)\n", "\n", "$$\n", "\\boxed{\n", "\\tilde{u}^{i} = \\frac{\\gamma}{\\left(W+B^{2}\\right)}\\left[\\tilde{Q}^{i} + \\frac{\\left(Q\\cdot B\\right)}{W}B^{i}\\right]\n", "}\\ ,\n", "$$\n", "\n", "where\n", "\n", "$$\n", "\\tilde{Q}^{i} = Q^{i} + \\left(Q\\cdot n\\right)n^{i}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", " // Recover the primitive variables from the scalars and conserved variables:\n", " gtmp = sqrt(1. - vsq);\n", " gamma = 1./gtmp ;\n", " rho0 = D * gtmp;\n", "\n", " w = W * (1. - vsq) ;\n", " p = pressure_rho0_w(eos, rho0,w) ;\n", " u = w - (rho0 + p) ; // u = rho0 eps, w = rho0 h\n", "\n", " if( (rho0 <= 0.) || (u <= 0.) ) {\n", " // User may want to handle this case differently, e.g. do NOT return upon\n", " // a negative rho/u, calculate v^i so that rho/u can be floored by other routine:\n", "\n", " retval = 5;\n", " //return(retval) ;\n", " }\n", "\n", " /*\n", " if(retval==5 && fabs(u)<1e-16) {\n", " u = fabs(u);\n", " CCTK_VInfo(CCTK_THORNSTRING,\"%e\\t%e\\t%e\",1.0-w/(rho0 + p),rho0,p);\n", " retval=0;\n", " }\n", " */\n", "\n", " prim[RHO] = rho0 ;\n", " prim[UU] = u ;\n", "\n", " for(i=1;i<4;i++) Qtcon[i] = Qcon[i] + ncon[i] * Qdotn;\n", " for(i=1;i<4;i++) prim[UTCON1+i-1] = gamma/(W+Bsq) * ( Qtcon[i] + QdotB*Bcon[i]/W ) ;\n", "\n", " /* set field components */\n", " for(i = BCON1; i <= BCON3; i++) prim[i] = U[i] ;\n", "\n", "\n", " /* done! */\n", " return(retval) ;\n", "\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.c: The `vsq_calc()` function \\[Back to [top](#toc)\\]\n", "$$\\label{vsq_calc}$$\n", "\n", "This function implements eq. (28) in [Noble *et al.* (2006)](https://arxiv.org/abs/astro-ph/0512420) to determine $v^{2}$:\n", "\n", "$$\n", "\\boxed{v^{2} = \\frac{\\tilde{Q}^{2}W^{2} + \\left(Q\\cdot B\\right)^{2}\\left(B^{2}+2W\\right)}{\\left(B^{2}+W\\right)^{2}W^{2}}}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", "\n", "/**********************************************************************/\n", "/****************************************************************************\n", " vsq_calc():\n", "\n", " -- evaluate v^2 (spatial, normalized velocity) from\n", " W = \\gamma^2 w\n", "\n", "****************************************************************************/\n", "static CCTK_REAL vsq_calc(CCTK_REAL W,CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D)\n", "{\n", " CCTK_REAL Wsq,Xsq;\n", "\n", " Wsq = W*W ;\n", " Xsq = (Bsq + W) * (Bsq + W);\n", "\n", " return( ( Wsq * Qtsq + QdotBsq * (Bsq + 2.*W)) / (Wsq*Xsq) );\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.d: The `x1_of_x0()` function \\[Back to [top](#toc)\\]\n", "$$\\label{x1_of_x0}$$\n", "\n", "This function computes $v^{2}$, as described [above](#vsq_calc), then performs physical checks on $v^{2}$ (i.e. whether or not it is superluminal). This function assumes $W$ is physical." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", "\n", "/********************************************************************\n", "\n", " x1_of_x0():\n", "\n", " -- calculates v^2 from W with some physical bounds checking;\n", " -- asumes x0 is already physical\n", " -- makes v^2 physical if not;\n", "\n", "*********************************************************************/\n", "\n", "static CCTK_REAL x1_of_x0(CCTK_REAL x0, CCTK_REAL &Bsq, CCTK_REAL &QdotBsq, CCTK_REAL &Qtsq, CCTK_REAL &Qdotn, CCTK_REAL &D )\n", "{\n", " CCTK_REAL vsq;\n", " CCTK_REAL dv = 1.e-15;\n", "\n", " vsq = fabs(vsq_calc(x0,Bsq,QdotBsq,Qtsq,Qdotn,D)) ; // guaranteed to be positive\n", "\n", "\n", " return( ( vsq > 1. ) ? (1.0 - dv) : vsq );\n", "\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.e: The `validate_x()` function \\[Back to [top](#toc)\\]\n", "$$\\label{validate_x}$$\n", "\n", "This function performs physical tests on $\\left\\{W,v^{2}\\right\\}$ based on their definitions." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", "/********************************************************************\n", "\n", " validate_x():\n", "\n", " -- makes sure that x[0,1] have physical values, based upon\n", " their definitions:\n", "\n", "*********************************************************************/\n", "\n", "static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] )\n", "{\n", "\n", " CCTK_REAL dv = 1.e-15;\n", "\n", " /* Always take the absolute value of x[0] and check to see if it's too big: */\n", " x[0] = fabs(x[0]);\n", " x[0] = (x[0] > W_TOO_BIG) ? x0[0] : x[0];\n", "\n", "\n", " x[1] = (x[1] < 0.) ? 0. : x[1]; /* if it's too small */\n", " x[1] = (x[1] > 1.) ? (1. - dv) : x[1]; /* if it's too big */\n", "\n", " return;\n", "\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.f: The `general_newton_raphson()` function \\[Back to [top](#toc)\\]\n", "$$\\label{general_newton_raphson}$$\n", "\n", "This function implements a [multidimensional Newton-Raphson method](https://en.wikipedia.org/wiki/Newton%27s_method#k_variables,_k_functions). We will not make the effort of explaining the algorithm exhaustively since it is pretty standard, so we will settle for a summary of the method.\n", "\n", "Given a system of $N$ non-linear of equations and $N$ variables, $\\left\\{\\vec{F}\\!\\left(\\vec{x}\\right),\\vec{x}\\right\\}$, the Newton-Raphson method attempts to determine the root vector, $\\vec{x}_{\\star}$, iteratively through\n", "\n", "$$\n", "\\begin{align}\n", "\\vec{x}_{n+1} = \\vec{x}_{n} - J^{-1}_{F}\\!\\left(\\vec{x}_{n}\\right)\\vec{F}\\!\\left(\\vec{x}\\right)\\ ,\n", "\\end{align}\n", "$$\n", "\n", "where $J^{-1}_{F}$ is the Jacobian matrix\n", "\n", "$$\n", "\\left(J_{F}\\right)^{i}_{\\ \\ j} = \\frac{\\partial F^{i}}{\\partial x^{j}}\\ .\n", "$$\n", "\n", "The index $n$ above is an *iteration* index and $\\vec{x}_{n+1}$ represents an improved approximation to $\\vec{x}_{\\star}$ when compared to $\\vec{x}_{n}$." ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", "/************************************************************\n", "\n", " general_newton_raphson():\n", "\n", " -- performs Newton-Rapshon method on an arbitrary system.\n", "\n", " -- inspired in part by Num. Rec.'s routine newt();\n", "\n", "*****************************************************************/\n", "static int general_newton_raphson( eos_struct eos, CCTK_REAL x[], int n, long &n_iter,\n", " void (*funcd) (eos_struct, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [],\n", " CCTK_REAL [][NEWT_DIM], CCTK_REAL *,\n", " CCTK_REAL *, int,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &),CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D)\n", "{\n", " CCTK_REAL f, df, dx[NEWT_DIM], x_old[NEWT_DIM];\n", " CCTK_REAL resid[NEWT_DIM], jac[NEWT_DIM][NEWT_DIM];\n", " CCTK_REAL errx, x_orig[NEWT_DIM];\n", " int id, i_extra, doing_extra;\n", "\n", " int keep_iterating;\n", "\n", "\n", " // Initialize various parameters and variables:\n", " errx = 1. ;\n", " df = f = 1.;\n", " i_extra = doing_extra = 0;\n", " for( id = 0; id < n ; id++) x_old[id] = x_orig[id] = x[id] ;\n", "\n", " n_iter = 0;\n", "\n", " /* Start the Newton-Raphson iterations : */\n", " keep_iterating = 1;\n", " while( keep_iterating ) {\n", "\n", " (*funcd) (eos, x, dx, resid, jac, &f, &df, n, Bsq,QdotBsq,Qtsq,Qdotn,D); /* returns with new dx, f, df */\n", "\n", "\n", " /* Save old values before calculating the new: */\n", " errx = 0.;\n", " for( id = 0; id < n ; id++) {\n", " x_old[id] = x[id] ;\n", " }\n", "\n", " /* Make the newton step: */\n", " for( id = 0; id < n ; id++) {\n", " x[id] += dx[id] ;\n", " }\n", "\n", " /****************************************/\n", " /* Calculate the convergence criterion */\n", " /****************************************/\n", " errx = (x[0]==0.) ? fabs(dx[0]) : fabs(dx[0]/x[0]);\n", "\n", "\n", " /****************************************/\n", " /* Make sure that the new x[] is physical : */\n", " /****************************************/\n", " validate_x( x, x_old ) ;\n", "\n", "\n", " /*****************************************************************************/\n", " /* If we've reached the tolerance level, then just do a few extra iterations */\n", " /* before stopping */\n", " /*****************************************************************************/\n", "\n", " if( (fabs(errx) <= NEWT_TOL) && (doing_extra == 0) && (EXTRA_NEWT_ITER > 0) ) {\n", " doing_extra = 1;\n", " }\n", "\n", " if( doing_extra == 1 ) i_extra++ ;\n", "\n", " if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0))\n", " || (i_extra > EXTRA_NEWT_ITER) || (n_iter >= (MAX_NEWT_ITER-1)) ) {\n", " keep_iterating = 0;\n", " }\n", "\n", " n_iter++;\n", "\n", " } // END of while(keep_iterating)\n", "\n", " /* Check for bad untrapped divergences : */\n", " if( (finite(f)==0) || (finite(df)==0) ) {\n", " return(2);\n", " }\n", "\n", "\n", " if( fabs(errx) > MIN_NEWT_TOL){\n", " //CCTK_VInfo(CCTK_THORNSTRING,\"%d %e %e %e %e\",n_iter,f,df,errx,MIN_NEWT_TOL);\n", " return(1);\n", " }\n", " if( (fabs(errx) <= MIN_NEWT_TOL) && (fabs(errx) > NEWT_TOL) ){\n", " return(0);\n", " }\n", " if( fabs(errx) <= NEWT_TOL ){\n", " return(0);\n", " }\n", "\n", " return(0);\n", "\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 2.g: The `func_vsq()` function \\[Back to [top](#toc)\\]\n", "$$\\label{func_vsq}$$\n", "\n", "This function is used by the `general_newton_raphson()` function to compute the residuals and stepping. We will again not describe it in great detail since the method itself is relatively straightforward." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", "\n", "\n", "/**********************************************************************/\n", "/*********************************************************************************\n", " func_vsq():\n", "\n", " -- calculates the residuals, and Newton step for general_newton_raphson();\n", " -- for this method, x=W,vsq here;\n", "\n", " Arguments:\n", " x = current value of independent var's (on input & output);\n", " dx = Newton-Raphson step (on output);\n", " resid = residuals based on x (on output);\n", " jac = Jacobian matrix based on x (on output);\n", " f = resid.resid/2 (on output)\n", " df = -2*f; (on output)\n", " n = dimension of x[];\n", "*********************************************************************************/\n", "\n", "static void func_vsq(eos_struct eos, CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],\n", " CCTK_REAL jac[][NEWT_DIM], CCTK_REAL *f, CCTK_REAL *df, int n,\n", " CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D)\n", "{\n", "\n", "\n", " CCTK_REAL W, vsq, Wsq, p_tmp, dPdvsq, dPdW;\n", " CCTK_REAL t11;\n", " CCTK_REAL t16;\n", " CCTK_REAL t18;\n", " CCTK_REAL t2;\n", " CCTK_REAL t21;\n", " CCTK_REAL t23;\n", " CCTK_REAL t24;\n", " CCTK_REAL t25;\n", " CCTK_REAL t3;\n", " CCTK_REAL t35;\n", " CCTK_REAL t36;\n", " CCTK_REAL t4;\n", " CCTK_REAL t40;\n", " CCTK_REAL t9;\n", "\n", " // vv TESTING vv\n", " // CCTK_REAL D,gtmp,gamma,rho0,w,p,u;\n", " // ^^ TESTING ^^\n", "\n", " W = x[0];\n", " vsq = x[1];\n", "\n", " Wsq = W*W;\n", "\n", " // vv TESTING vv\n", " /*\n", " D = U[RHO] ;\n", " gtmp = sqrt(1. - vsq);\n", " gamma = 1./gtmp ;\n", " rho0 = D * gtmp;\n", "\n", " w = W * (1. - vsq) ;\n", " p = pressure_rho0_w(eos, rho0,w) ;\n", " u = w - (rho0 + p) ;\n", "\n", " if(u<=0 && 1==1) {\n", " vsq = 0.9999999 * (1.0-(rho0+p)/W);\n", "\n", " w = W * (1. - vsq) ;\n", " p = pressure_rho0_w(eos, rho0,w) ;\n", " u = w - (rho0 + p) ;\n", "\n", " //CCTK_VInfo(CCTK_THORNSTRING,\"%e check\",u);\n", " }\n", " */\n", " // ^^ TESTING ^^\n", "\n", "\n", " p_tmp = pressure_W_vsq( eos, W, vsq , D);\n", " dPdW = dpdW_calc_vsq( W, vsq );\n", " dPdvsq = dpdvsq_calc( eos, W, vsq, D );\n", "\n", " // These expressions were calculated using Mathematica, but made into efficient\n", " // code using Maple. Since we know the analytic form of the equations, we can\n", " // explicitly calculate the Newton-Raphson step:\n", "\n", " t2 = -0.5*Bsq+dPdvsq;\n", " t3 = Bsq+W;\n", " t4 = t3*t3;\n", " t9 = 1/Wsq;\n", " t11 = Qtsq-vsq*t4+QdotBsq*(Bsq+2.0*W)*t9;\n", " t16 = QdotBsq*t9;\n", " t18 = -Qdotn-0.5*Bsq*(1.0+vsq)+0.5*t16-W+p_tmp;\n", " t21 = 1/t3;\n", " t23 = 1/W;\n", " t24 = t16*t23;\n", " t25 = -1.0+dPdW-t24;\n", " t35 = t25*t3+(Bsq-2.0*dPdvsq)*(QdotBsq+vsq*Wsq*W)*t9*t23;\n", " t36 = 1/t35;\n", " dx[0] = -(t2*t11+t4*t18)*t21*t36;\n", " t40 = (vsq+t24)*t3;\n", " dx[1] = -(-t25*t11-2.0*t40*t18)*t21*t36;\n", " //detJ = t3*t35; // <- set but not used...\n", " jac[0][0] = -2.0*t40;\n", " jac[0][1] = -t4;\n", " jac[1][0] = t25;\n", " jac[1][1] = t2;\n", " resid[0] = t11;\n", " resid[1] = t18;\n", "\n", "\n", "\n", " *df = -resid[0]*resid[0] - resid[1]*resid[1];\n", "\n", " *f = -0.5 * ( *df );\n", "\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 3: EOS dependent routines \\[Back to [top](#toc)\\]\n", "$$\\label{harm_utoprim_2d__c__eos_dep}$$" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", "\n", "/**********************************************************************\n", " **********************************************************************\n", "\n", " The following routines specify the equation of state. All routines\n", " above here should be indpendent of EOS. If the user wishes\n", " to use another equation of state, the below functions must be replaced\n", " by equivalent routines based upon the new EOS.\n", "\n", " **********************************************************************\n", " **********************************************************************/" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.a: The `pressure_W_vsq()` function \\[Back to [top](#toc)\\]\n", "$$\\label{pressure_w_vsq}$$\n", "\n", "This function computes $p\\left(W,v^{2}\\right)$. For a $\\Gamma$-law equation of state,\n", "\n", "$$\n", "p_{\\Gamma} = \\left(\\Gamma-1\\right)u\\ ,\n", "$$\n", "\n", "and with the definitions\n", "\n", "$$\n", "\\begin{align}\n", "\\gamma^{2} &= \\frac{1}{1-v^{2}}\\ ,\\\\\n", "W &= \\gamma^{2}w\\ ,\\\\\n", "D &= \\gamma\\rho_{b}\\ ,\\\\\n", "w &= \\rho_{b} + u + p\\ ,\n", "\\end{align}\n", "$$\n", "\n", "we have\n", "\n", "$$\n", "\\begin{align}\n", "p_{\\Gamma} &= \\left(\\Gamma-1\\right)u\\\\\n", " &= \\left(\\Gamma-1\\right)\\left(w - \\rho_{b} - p_{\\Gamma}\\right)\\\\\n", " &= \\left(\\Gamma-1\\right)\\left(\\frac{W}{\\gamma^{2}} - \\frac{D}{\\gamma}\\right) - \\left(\\Gamma-1\\right)p_{\\Gamma}\\\\\n", "\\implies\n", "&\\boxed{\n", "p_{\\Gamma} = \\frac{\\left(\\Gamma-1\\right)}{\\Gamma}\\left(\\frac{W}{\\gamma^{2}} - \\frac{D}{\\gamma}\\right)\n", "}\\ .\n", "\\end{align}\n", "$$\n", "\n", "Thus, the pre-PPEOS Patch version of this function was\n", "\n", "```c\n", "/**********************************************************************/\n", "/********************************************************************** \n", " pressure_W_vsq(): \n", " \n", " -- Gamma-law equation of state;\n", " -- pressure as a function of W, vsq, and D:\n", "**********************************************************************/\n", "static CCTK_REAL pressure_W_vsq(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) \n", "{\n", "\n", " CCTK_REAL gtmp;\n", " gtmp = 1. - vsq;\n", " \n", " return( (GAMMA - 1.) * ( W * gtmp - D * sqrt(gtmp) ) / GAMMA );\n", "\n", "}\n", "```\n", "\n", "We are now, however, interested in the hybrid EOS of the form\n", "\n", "$$\n", "p_{\\rm hybrid} = P_{\\rm cold} + P_{\\rm th}\\ ,\n", "$$\n", "\n", "where $P_{\\rm cold}$ is given by a single or piecewise polytrope EOS,\n", "\n", "$$\n", "P_{\\rm cold} = K_{i}\\rho_{b}^{\\Gamma_{i}}\\ ,\n", "$$\n", "\n", "$P_{\\rm th}$ accounts for thermal effects and is given by\n", "\n", "$$\n", "P_{\\rm th} = \\left(\\Gamma_{\\rm th} - 1\\right)\\epsilon_{\\rm th}\\ ,\n", "$$\n", "\n", "and\n", "\n", "$$\n", "\\begin{align}\n", "\\epsilon \\equiv \\frac{u}{\\rho_{b}} &= \\epsilon_{\\rm th}+\\epsilon_{\\rm cold}\\ ,\\\\\n", "\\epsilon_{\\rm cold} &= \\int d\\rho \\frac{P_{\\rm cold}(\\rho)}{\\rho^{2}}\\ .\n", "\\end{align}\n", "$$\n", "\n", "We then have\n", "\n", "$$\n", "\\begin{align}\n", "p_{\\rm hybrid} &= P_{\\rm cold} + P_{\\rm th}\\\\\n", " &= P_{\\rm cold} + \\left(\\Gamma_{\\rm th}-1\\right)\\rho_{b}\\epsilon_{\\rm th}\\\\\n", " &= P_{\\rm cold} + \\left(\\Gamma_{\\rm th}-1\\right)\\rho_{b}\\left(\\epsilon - \\epsilon_{\\rm cold}\\right)\\\\\n", " &= P_{\\rm cold} + \\left(\\Gamma_{\\rm th}-1\\right)\\left(u - \\frac{D}{\\gamma}\\epsilon_{\\rm cold}\\right)\\\\\n", " &= P_{\\rm cold} + \\left(\\Gamma_{\\rm th}-1\\right)\\left(w - \\rho_{b} - p_{\\rm hybrid} - \\frac{D}{\\gamma}\\epsilon_{\\rm cold}\\right)\\\\\n", " &= P_{\\rm cold} + \\left(\\Gamma_{\\rm th}-1\\right)\\left(\\frac{W}{\\gamma^{2}} - \\frac{D}{\\gamma} - \\frac{D}{\\gamma}\\epsilon_{\\rm cold}\\right)-\\left(\\Gamma_{\\rm th}-1\\right)p_{\\rm hybrid}\\\\\n", " &= P_{\\rm cold} + \\left(\\Gamma_{\\rm th}-1\\right)\\left[\\frac{W}{\\gamma^{2}} - \\frac{D}{\\gamma}\\left(1+\\epsilon_{\\rm cold}\\right)\\right]-\\left(\\Gamma_{\\rm th}-1\\right)p_{\\rm hybrid}\\\\\n", "\\implies\n", "&\\boxed{ p_{\\rm hybrid} = \\frac{P_{\\rm cold}}{\\Gamma_{\\rm th}} + \\frac{\\left(\\Gamma_{\\rm th}-1\\right)}{\\Gamma_{\\rm th}}\\left[\\frac{W}{\\gamma^{2}} - \\frac{D}{\\gamma}\\left(1+\\epsilon_{\\rm cold}\\right)\\right] }\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", "/**********************************************************************/\n", "/**********************************************************************\n", " pressure_W_vsq():\n", "\n", " -- Hybrid single and piecewise polytropic equation of state;\n", " -- pressure as a function of P_cold, eps_cold, W, vsq, and D:\n", "**********************************************************************/\n", "static CCTK_REAL pressure_W_vsq(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D)\n", "{\n", "\n", "#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", " DECLARE_CCTK_PARAMETERS;\n", "#endif\n", "\n", " // Compute gamma^{-2} = 1 - v^{2} and gamma^{-1}\n", " CCTK_REAL inv_gammasq = 1.0 - vsq;\n", " CCTK_REAL inv_gamma = sqrt(inv_gammasq);\n", "\n", " // Compute rho_b = D / gamma\n", " CCTK_REAL rho_b = D*inv_gamma;\n", "\n", " // Compute P_cold and eps_cold\n", " CCTK_REAL P_cold, eps_cold;\n", " compute_P_cold__eps_cold(eos,rho_b, P_cold,eps_cold);\n", "\n", " // Compute p = P_{cold} + P_{th}\n", " return( ( P_cold + (Gamma_th - 1.0)*( W*inv_gammasq - D*inv_gamma*( 1.0 + eps_cold ) ) )/Gamma_th );\n", "\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.b: The `dpdW_calc_vsq()` function \\[Back to [top](#toc)\\]\n", "$$\\label{dpdw_calc_vsq}$$\n", "\n", "This function computes $\\frac{\\partial p\\left(W,v^{2}\\right)}{\\partial W}$. For a $\\Gamma$-law equation of state, remember that\n", "\n", "$$\n", "p_{\\Gamma} = \\frac{\\left(\\Gamma-1\\right)}{\\Gamma}\\left(\\frac{W}{\\gamma^{2}} - \\frac{D}{\\gamma}\\right)\\ ,\n", "$$\n", "\n", "which then implies\n", "\n", "$$\n", "\\boxed{\\frac{\\partial p_{\\Gamma}}{\\partial W} = \\frac{\\Gamma-1}{\\Gamma \\gamma^{2}} = \\frac{\\left(\\Gamma-1\\right)\\left(1-v^{2}\\right)}{\\Gamma}}\\ .\n", "$$\n", "\n", "Thus, the pre-PPEOS Patch version of this function was\n", "\n", "```c\n", "/**********************************************************************/\n", "/********************************************************************** \n", " dpdW_calc_vsq(): \n", " \n", " -- partial derivative of pressure with respect to W;\n", "**********************************************************************/\n", "static CCTK_REAL dpdW_calc_vsq(CCTK_REAL W, CCTK_REAL vsq)\n", "{\n", "\n", " return( (GAMMA - 1.) * (1. - vsq) / GAMMA ) ;\n", "\n", "}\n", "```\n", "\n", "For the case of a hybrid, single or piecewise polytropic EOS, we have\n", "\n", "$$\n", "p_{\\rm hybrid} = \\frac{P_{\\rm cold}}{\\Gamma_{\\rm th}} + \\frac{\\left(\\Gamma_{\\rm th}-1\\right)}{\\Gamma_{\\rm th}}\\left[\\frac{W}{\\gamma^{2}} - \\frac{D}{\\gamma}\\left(1+\\epsilon_{\\rm cold}\\right)\\right]\\ .\n", "$$\n", "\n", "It is important to notice that the cold components of $p_{\\rm hybrid}$ are *not* functions of $W$, but instead functions of $D$: $P_{\\rm cold} = P_{\\rm cold}(\\rho_{b}) = P_{\\rm cold}(D)$ and $\\epsilon_{\\rm cold} = \\epsilon_{\\rm cold}(\\rho_{b}) = \\epsilon_{\\rm cold}(D)$. Thus\n", "\n", "$$\n", "\\boxed{\\frac{\\partial p_{\\rm hybrid}}{\\partial W} = \\frac{\\Gamma_{\\rm th}-1}{\\Gamma_{\\rm th} \\gamma^{2}} = \\frac{\\left(\\Gamma_{\\rm th}-1\\right)\\left(1-v^{2}\\right)}{\\Gamma_{\\rm th}}}\\ .\n", "$$" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", "\n", "/**********************************************************************/\n", "/**********************************************************************\n", " dpdW_calc_vsq():\n", "\n", " -- partial derivative of pressure with respect to W;\n", "**********************************************************************/\n", "static CCTK_REAL dpdW_calc_vsq(CCTK_REAL W, CCTK_REAL vsq)\n", "{\n", "\n", "#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", " DECLARE_CCTK_PARAMETERS;\n", "#endif\n", "\n", " return( (Gamma_th - 1.0) * (1.0 - vsq) / Gamma_th ) ;\n", "\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Step 3.c: The `dpdvsq_calc()` function \\[Back to [top](#toc)\\]\n", "$$\\label{dpdvsq_calc}$$\n", "\n", "This function computes $\\frac{\\partial p\\left(W,v^{2}\\right)}{\\partial W}$. For a $\\Gamma$-law equation of state, remember that\n", "\n", "$$\n", "p_{\\Gamma} = \\frac{\\left(\\Gamma-1\\right)}{\\Gamma}\\left(\\frac{W}{\\gamma^{2}} - \\frac{D}{\\gamma}\\right) = \\frac{\\left(\\Gamma-1\\right)}{\\Gamma}\\left[W\\left(1-v^{2}\\right) - D\\sqrt{1-v^{2}}\\right]\\ ,\n", "$$\n", "\n", "which then implies\n", "\n", "$$\n", "\\boxed{\\frac{\\partial p_{\\Gamma}}{\\partial\\left(v^{2}\\right)} = \\frac{\\Gamma-1}{\\Gamma}\\left(\\frac{D}{2\\sqrt{1-v^{2}}}-W\\right)} \\ .\n", "$$\n", "\n", "Thus, the pre-PPEOS Patch version of this function was\n", "\n", "```c\n", "/**********************************************************************/\n", "/********************************************************************** \n", " dpdvsq_calc(): \n", " \n", " -- partial derivative of pressure with respect to vsq\n", "**********************************************************************/\n", "static CCTK_REAL dpdvsq_calc(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D)\n", "{\n", " return( (GAMMA - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / GAMMA ) ;\n", "}\n", "```\n", "\n", "\n", "\n", "### Step 3.c.i: Setting basic quantities and computing $P_{\\rm cold}$ and $\\epsilon_{\\rm cold}$ \\[Back to [top](#toc)\\]\n", "$$\\label{dpdvsq_calc__basic_quantities}$$\n", "\n", "For the case of a hybrid, single or piecewise polytropic EOS, we have\n", "\n", "$$\n", "p_{\\rm hybrid} = \\frac{P_{\\rm cold}}{\\Gamma_{\\rm th}} + \\frac{\\left(\\Gamma_{\\rm th}-1\\right)}{\\Gamma_{\\rm th}}\\left[\\frac{W}{\\gamma^{2}} - \\frac{D}{\\gamma}\\left(1+\\epsilon_{\\rm cold}\\right)\\right]\\ .\n", "$$\n", "\n", "Let us thus begin by setting the necessary parameters from the hybrid EOS." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", "/**********************************************************************/\n", "/**********************************************************************\n", " dpdvsq_calc():\n", "\n", " -- partial derivative of pressure with respect to vsq\n", "**********************************************************************/\n", "static CCTK_REAL dpdvsq_calc(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D)\n", "{\n", "\n", " // This sets Gamma_th\n", "#ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", " DECLARE_CCTK_PARAMETERS;\n", "#endif\n", "\n", "\n", " // Set gamma and rho\n", " CCTK_REAL gamma = 1.0/sqrt(1.0 - vsq);\n", " CCTK_REAL rho_b = D/gamma;\n", "\n", " // Compute P_cold and eps_cold\n", " CCTK_REAL P_cold, eps_cold;\n", " compute_P_cold__eps_cold(eos,rho_b, P_cold,eps_cold);\n", "\n", " // Set basic polytropic quantities\n", " int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b);\n", " CCTK_REAL Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.c.ii: Computing $\\frac{\\partial P_{\\rm cold}}{\\partial\\left(v^{2}\\right)}$ \\[Back to [top](#toc)\\]\n", "$$\\label{dpdvsq_calc__dpcolddvsq}$$\n", "\n", "Next, remember that $P_{\\rm cold} = P_{\\rm cold}(\\rho_{b}) = P_{\\rm cold}(D,v^{2})$ and also $\\epsilon_{\\rm cold} = \\epsilon_{\\rm cold}(D,v^{2})$. Therefore, we must start by finding the derivatives of $P_{\\rm cold}$ and $\\epsilon_{\\rm cold}$ with respect to $v^{2}$.\n", "\n", "Let us first notice that\n", "\n", "$$\n", "\\frac{\\partial\\gamma}{\\partial\\left(v^{2}\\right)} = \\frac{\\partial}{\\partial\\left(v^{2}\\right)}\\left[\\frac{1}{\\sqrt{1-v^{2}}}\\right] = \\frac{1}{2}\\left(1-v^{2}\\right)^{-3/2} = \\frac{\\gamma^{3}}{2}\\ .\n", "$$\n", "\n", "Thus, for a general power\n", "\n", "$$\n", "\\frac{\\partial\\gamma^{a}}{\\partial\\left(v^{2}\\right)} = a\\gamma^{a-1}\\frac{\\partial\\gamma}{\\partial\\left(v^{2}\\right)} = a\\gamma^{a-1}\\left(\\frac{\\gamma^{3}}{2}\\right) = \\frac{a}{2}\\gamma^{a+2}\n", "$$\n", "\n", "Thus we have\n", "\n", "$$\n", "\\begin{align}\n", "\\frac{\\partial P_{\\rm cold}}{\\partial \\left(v^{2}\\right)}\n", "&= \\frac{\\partial}{\\partial\\left(v^{2}\\right)}\\left(K_{\\rm poly}\\rho_{b}^{\\Gamma_{\\rm poly}}\\right)\\\\\n", "&= \\frac{\\partial}{\\partial\\left(v^{2}\\right)}\\left[K_{\\rm poly}\\left(\\frac{D}{\\gamma}\\right)^{\\Gamma_{\\rm poly}}\\right]\\\\\n", "&= K_{\\rm poly}D^{\\Gamma_{\\rm poly}}\\frac{\\partial}{\\partial\\left(v^{2}\\right)}\\left[\\gamma^{-\\Gamma_{\\rm poly}/2}\\right]\\\\\n", "&=K_{\\rm poly}D^{\\Gamma_{\\rm poly}}\\left[\\frac{-\\Gamma_{\\rm poly}/2}{2}\\gamma^{-\\Gamma_{\\rm poly}/2 + 2}\\right]\\\\\n", "&=K_{\\rm poly}\\left(\\frac{D}{\\gamma}\\right)^{\\Gamma_{\\rm poly}}\\gamma^{-\\frac{\\Gamma_{\\rm poly}}{2} + 2 + \\Gamma_{\\rm poly}}\\\\\n", "\\implies &\\boxed{ \\frac{\\partial P_{\\rm cold}}{\\partial \\left(v^{2}\\right)} = \\gamma^{2+\\frac{\\Gamma_{\\rm poly}}{2}}P_{\\rm cold}}\\ .\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", " /* Now we implement the derivative of P_cold with respect\n", " * to v^{2}, given by\n", " * ----------------------------------------------------\n", " * | dP_cold/dvsq = gamma^{2 + Gamma_{poly}/2} P_{cold} |\n", " * ----------------------------------------------------\n", " */\n", " CCTK_REAL dPcold_dvsq = P_cold * pow(gamma,2.0 + 0.5*Gamma_ppoly_tab);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.c.iii: Computing $\\frac{\\partial \\epsilon_{\\rm cold}}{\\partial\\left(v^{2}\\right)}$ \\[Back to [top](#toc)\\]\n", "$$\\label{dpdvsq_calc__depscolddvsq}$$\n", "\n", "Now, obtaining $\\epsilon_{\\rm cold}$ from $P_{\\rm cold}$ requires an integration and, therefore, generates an integration constant. Since we are interested in a *derivative* of $\\epsilon_{\\rm cold}$, however, we will simply drop the constant altogether. Remember that:\n", "\n", "$$\n", "\\epsilon_{\\rm cold} = K_{\\rm poly}\\int d\\rho_{b} \\rho_{b}^{\\Gamma_{\\rm poly}-2} = \\frac{K_{\\rm poly}\\rho_{b}^{\\Gamma_{\\rm poly}-1}}{\\Gamma_{\\rm poly}-1} = \\frac{P_{\\rm cold}}{\\rho_{b}\\left(\\Gamma_{\\rm poly}-1\\right)} = \\frac{\\gamma P_{\\rm cold}}{D\\left(\\Gamma_{\\rm poly}-1\\right)}\\ .\n", "$$\n", "\n", "Thus\n", "\n", "$$\n", "\\begin{align}\n", "\\frac{\\partial \\epsilon_{\\rm cold}}{\\partial \\left(v^{2}\\right)}\n", "&= \\frac{1}{D\\left(\\Gamma_{\\rm poly}-1\\right)}\\left[\\gamma\\frac{\\partial P_{\\rm cold}}{\\partial \\left(v^{2}\\right)} + P_{\\rm cold}\\frac{\\partial\\gamma}{\\partial \\left(v^{2}\\right)}\\right]\\\\\n", "&=\\frac{1}{D\\left(\\Gamma_{\\rm poly}-1\\right)}\\left[\\gamma\\frac{\\partial P_{\\rm cold}}{\\partial \\left(v^{2}\\right)} + P_{\\rm cold}\\left(\\frac{\\gamma^{3}}{2}\\right)\\right]\\\\\n", "\\implies &\\boxed{\n", "\\frac{\\partial \\epsilon_{\\rm cold}}{\\partial \\left(v^{2}\\right)} = \\frac{\\gamma}{D\\left(\\Gamma_{\\rm poly}-1\\right)}\\left[\\frac{\\partial P_{\\rm cold}}{\\partial \\left(v^{2}\\right)} + \\frac{\\gamma^{2} P_{\\rm cold}}{2}\\right]\\ .\n", "}\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", "\n", " /* Now we implement the derivative of eps_cold with respect\n", " * to v^{2}, given by\n", " * -----------------------------------------------------------------------------------\n", " * | deps_cold/dvsq = gamma/(D*(Gamma_ppoly_tab-1)) * (dP_cold/dvsq + gamma^{2} P_cold / 2) |\n", " * -----------------------------------------------------------------------------------\n", " */\n", " CCTK_REAL depscold_dvsq = ( gamma/(D*(Gamma_ppoly_tab-1.0)) ) * ( dPcold_dvsq + 0.5*gamma*gamma*P_cold );" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "### Step 3.c.iv: Computing $\\frac{\\partial p_{\\rm hybrid}}{\\partial\\left(v^{2}\\right)}$ \\[Back to [top](#toc)\\]\n", "$$\\label{dpdvsq_calc__dpdvsq}$$\n", "\n", "Finally, remembering that\n", "\n", "$$\n", "\\begin{align}\n", "p_{\\rm hybrid} &= \\frac{P_{\\rm cold}}{\\Gamma_{\\rm th}} + \\frac{\\left(\\Gamma_{\\rm th}-1\\right)}{\\Gamma_{\\rm th}}\\left[\\frac{W}{\\gamma^{2}} - \\frac{D}{\\gamma}\\left(1+\\epsilon_{\\rm cold}\\right)\\right]\\ ,\\\\\n", "\\frac{\\partial\\gamma^{a}}{\\partial\\left(v^{2}\\right)} &= \\frac{a}{2}\\gamma^{a+2}\\ ,\n", "\\end{align}\n", "$$\n", "\n", "we have\n", "\n", "$$\n", "\\boxed{\n", "\\frac{\\partial p_{\\rm hybrid}}{\\partial\\left(v^{2}\\right)}\n", "= \\frac{1}{\\Gamma_{\\rm th}}\\left\\{\\frac{\\partial P_{\\rm cold}}{\\partial\\left(v^{2}\\right)} + \\left(\\Gamma_{\\rm th}-1\\right)\\left[-W + \\frac{D\\gamma}{2}\\left(1+\\epsilon_{\\rm cold}\\right) - \\frac{D}{\\gamma}\\frac{\\partial \\epsilon_{\\rm cold}}{\\partial\\left(v^{2}\\right)}\\right]\\right\\}\\ .\n", "}\n", "$$" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Appending to ../src/harm_utoprim_2d.c\n" ] } ], "source": [ "%%writefile -a $outfile_path__harm_utoprim_2d__c\n", "\n", " /* Now we implement the derivative of p_hybrid with respect\n", " * to v^{2}, given by\n", " * -----------------------------------------------------------------------------\n", " * | dp/dvsq = Gamma_th^{-1}( dP_cold/dvsq |\n", " * | + (Gamma_{th}-1)*(-W |\n", " * | + D gamma (1 + eps_cold)/2 |\n", " * | - (D/gamma) * deps_cold/dvsq) ) |\n", " * -----------------------------------------------------------------------------\n", " */\n", " return( ( dPcold_dvsq + (Gamma_th-1.0)*( -W + D*gamma*(1+eps_cold)/2.0 - D*depscold_dvsq/gamma ) )/Gamma_th );\n", "}\n", "\n", "\n", "/******************************************************************************\n", " END OF UTOPRIM_2D.C\n", "******************************************************************************/\n", "#endif\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 4: Code validation \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation}$$\n", "\n", "First we download the original `IllinoisGRMHD` source code and then compare it to the source code generated by this tutorial notebook." ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Validation test for harm_utoprim_2d.c: FAILED!\n", "Diff:\n", "0a1,2\n", "> #ifndef __HARM_UTOPRIM_2D__C__\n", "> #define __HARM_UTOPRIM_2D__C__\n", "2c4\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", "7c9\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,10c11,12\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", "12c14\n", "< an accretion disk model. \n", "---\n", "> an accretion disk model.\n", "14c16\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", "17c19\n", "< [1] Gammie, C. F., McKinney, J. C., \\& Toth, G.\\ 2003, \n", "---\n", "> [1] Gammie, C. F., McKinney, J. C., \\& Toth, G.\\ 2003,\n", "20c22\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,24c25,26\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", "49c51\n", "< utoprim_2d.c: \n", "---\n", "> utoprim_2d.c:\n", "52c54\n", "< Uses the 2D method: \n", "---\n", "> Uses the 2D method:\n", "54,55c56,57\n", "< Newton-Raphson method \n", "< -- can be used (in principle) with a general equation of state. \n", "---\n", "> Newton-Raphson method\n", "> -- can be used (in principle) with a general equation of state.\n", "58,60c60,62\n", "< density or internal energy density is calculated. You may want \n", "< to change this aspect of the code so that it still calculates the \n", "< velocity and so that you can floor the densities. If you want to \n", "---\n", "> density or internal energy density is calculated. You may want\n", "> to change this aspect of the code so that it still calculates the\n", "> velocity and so that you can floor the densities. If you want to\n", "68c70\n", "< // Declarations: \n", "---\n", "> // Declarations:\n", "70,72c72,74\n", "< static int Utoprim_new_body(CCTK_REAL U[], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL gdet, CCTK_REAL prim[],long &n_iter);\n", "< static int general_newton_raphson( CCTK_REAL x[], int n, long &n_iter, void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *, CCTK_REAL *, int,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &),CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D);\n", "< static void func_vsq( CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *f, CCTK_REAL *df, int n,CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D);\n", "---\n", "> static int Utoprim_new_body(eos_struct eos, CCTK_REAL U[], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM], CCTK_REAL gdet, CCTK_REAL prim[],long &n_iter);\n", "> static int general_newton_raphson( eos_struct eos, CCTK_REAL x[], int n, long &n_iter, void (*funcd) (eos_struct, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *, CCTK_REAL *, int,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &,CCTK_REAL &),CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D);\n", "> static void func_vsq( eos_struct eos, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], CCTK_REAL [][NEWT_DIM], CCTK_REAL *f, CCTK_REAL *df, int n,CCTK_REAL &Bsq,CCTK_REAL &QdotBsq,CCTK_REAL &Qtsq,CCTK_REAL &Qdotn,CCTK_REAL &D);\n", "74c76\n", "< static CCTK_REAL pressure_W_vsq(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) ;\n", "---\n", "> static CCTK_REAL pressure_W_vsq(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) ;\n", "76c78\n", "< static CCTK_REAL dpdvsq_calc(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D);\n", "---\n", "> static CCTK_REAL dpdvsq_calc(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D);\n", "82c84\n", "< \n", "---\n", "> \n", "84c86\n", "< between the two sets of definitions for U and P. The user may \n", "---\n", "> between the two sets of definitions for U and P. The user may\n", "89,97c91,99\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", "101c103\n", "< U[NPR] = conserved variables (current values on input/output);\n", "---\n", "> U[NPR] = conserved variables (current values on input/output);\n", "105,107c107,109\n", "< prim[NPR] = primitive variables (guess on input, calculated values on \n", "< output if there are no problems);\n", "< \n", "---\n", "> prim[NPR] = primitive variables (guess on input, calculated values on\n", "> output if there are no problems);\n", "> \n", "109c111\n", "< unfamiliar with metrics, merely set \n", "---\n", "> unfamiliar with metrics, merely set\n", "114c116,117\n", "< int Utoprim_2d(CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM], \n", "---\n", "> \n", "> int Utoprim_2d(eos_struct eos, CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], CCTK_REAL gcon[NDIM][NDIM],\n", "119c122\n", "< int i, ret; \n", "---\n", "> int i, ret;\n", "122c125\n", "< if( U[0] <= 0. ) { \n", "---\n", "> if( U[0] <= 0. ) {\n", "131c134,135\n", "< \n", "---\n", "> \n", "> \n", "141a146\n", "> \n", "150c155\n", "< ret = Utoprim_new_body(U_tmp, gcov, gcon, gdet, prim_tmp,n_iter);\n", "---\n", "> ret = Utoprim_new_body(eos, U_tmp, gcov, gcon, gdet, prim_tmp,n_iter);\n", "152c157\n", "< /* Transform new primitive variables back if there was no problem : */ \n", "---\n", "> /* Transform new primitive variables back if there was no problem : */\n", "163a169\n", "> \n", "171,172c177,178\n", "< -- This is the main routine that calculates auxiliary quantities for the \n", "< Newton-Raphson routine. \n", "---\n", "> -- This is the main routine that calculates auxiliary quantities for the\n", "> Newton-Raphson routine.\n", "174,177c180,183\n", "< -- assumes that \n", "< / rho gamma \\\n", "< U = | alpha T^t_\\mu |\n", "< \\ alpha B^i /\n", "---\n", "> -- assumes that\n", "> / rho gamma \\\n", "> U = | alpha T^t_\\mu |\n", "> \\ alpha B^i /\n", "181,184c187,190\n", "< / rho \\\n", "< prim = | uu |\n", "< | \\tilde{u}^i |\n", "< \\ alpha B^i /\n", "---\n", "> / rho \\\n", "> prim = | uu |\n", "> | \\tilde{u}^i |\n", "> \\ alpha B^i /\n", "187,188c193,194\n", "< return: (i*100 + j) where \n", "< i = 0 -> Newton-Raphson solver either was not called (yet or not used) \n", "---\n", "> return: (i*100 + j) where\n", "> i = 0 -> Newton-Raphson solver either was not called (yet or not used)\n", "190c196\n", "< 1 -> Newton-Raphson solver did not converge to a solution with the \n", "---\n", "> 1 -> Newton-Raphson solver did not converge to a solution with the\n", "192c198\n", "< 2 -> Newton-Raphson procedure encountered a numerical divergence \n", "---\n", "> 2 -> Newton-Raphson procedure encountered a numerical divergence\n", "194,196c200,202\n", "< \n", "< j = 0 -> success \n", "< 1 -> failure: some sort of failure in Newton-Raphson; \n", "---\n", "> \n", "> j = 0 -> success\n", "> 1 -> failure: some sort of failure in Newton-Raphson;\n", "198,199c204,205\n", "< 3 -> failure: W<0 or W>W_TOO_BIG\n", "< 4 -> failure: v^2 > 1 \n", "---\n", "> 3 -> failure: W<0 or W>W_TOO_BIG\n", "> 4 -> failure: v^2 > 1\n", "204c210\n", "< static int Utoprim_new_body(CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM], \n", "---\n", "> static int Utoprim_new_body(eos_struct eos, CCTK_REAL U[NPR], CCTK_REAL gcov[NDIM][NDIM],\n", "217a224\n", "> \n", "237a245,248\n", "> // FIXME: The exact form of n^{\\mu} can be found\n", "> // in eq. (2.116) and implementing it\n", "> // directly is a lot more efficient than\n", "> // performing n^{\\mu} = g^{\\mu\\nu}n_{nu}\n", "248a260\n", "> \n", "255c267\n", "< if( (utsq < 0.) && (fabs(utsq) < 1.0e-13) ) { \n", "---\n", "> if( (utsq < 0.) && (fabs(utsq) < 1.0e-13) ) {\n", "265c277\n", "< \n", "---\n", "> \n", "267c279\n", "< // i.e. you don't get positive values for dP/d(vsq) . \n", "---\n", "> // i.e. you don't get positive values for dP/d(vsq) .\n", "270c282\n", "< p = pressure_rho0_u(rho0,u) ;\n", "---\n", "> p = pressure_rho0_u(eos, rho0,u) ;\n", "276c288\n", "< // Make sure that W is large enough so that v^2 < 1 : \n", "---\n", "> // Make sure that W is large enough so that v^2 < 1 :\n", "278c290\n", "< while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq ) \n", "---\n", "> while( (( W_last*W_last*W_last * ( W_last + 2.*Bsq )\n", "284,285c296,298\n", "< \n", "< // Calculate W and vsq: \n", "---\n", "> \n", "> \n", "> // Calculate W and vsq:\n", "288c301\n", "< retval = general_newton_raphson( x_2d, n, n_iter, func_vsq, Bsq,QdotBsq,Qtsq,Qdotn,D) ; \n", "---\n", "> retval = general_newton_raphson( eos, x_2d, n, n_iter, func_vsq, Bsq,QdotBsq,Qtsq,Qdotn,D) ;\n", "292c305\n", "< \n", "---\n", "> \n", "311a325\n", "> \n", "318c332\n", "< p = pressure_rho0_w(rho0,w) ;\n", "---\n", "> p = pressure_rho0_w(eos, rho0,w) ;\n", "322c336\n", "< // User may want to handle this case differently, e.g. do NOT return upon \n", "---\n", "> // User may want to handle this case differently, e.g. do NOT return upon\n", "342c356\n", "< \n", "---\n", "> \n", "353c367,368\n", "< /**********************************************************************/ \n", "---\n", "> \n", "> /**********************************************************************/\n", "355,358c370,373\n", "< vsq_calc(): \n", "< \n", "< -- evaluate v^2 (spatial, normalized velocity) from \n", "< W = \\gamma^2 w \n", "---\n", "> vsq_calc():\n", "> \n", "> -- evaluate v^2 (spatial, normalized velocity) from\n", "> W = \\gamma^2 w\n", "364c379\n", "< \n", "---\n", "> \n", "371a387\n", "> \n", "374,375c390,391\n", "< x1_of_x0(): \n", "< \n", "---\n", "> x1_of_x0():\n", "> \n", "382c398\n", "< static CCTK_REAL x1_of_x0(CCTK_REAL x0, CCTK_REAL &Bsq, CCTK_REAL &QdotBsq, CCTK_REAL &Qtsq, CCTK_REAL &Qdotn, CCTK_REAL &D ) \n", "---\n", "> static CCTK_REAL x1_of_x0(CCTK_REAL x0, CCTK_REAL &Bsq, CCTK_REAL &QdotBsq, CCTK_REAL &Qtsq, CCTK_REAL &Qdotn, CCTK_REAL &D )\n", "386,387d401\n", "< \n", "< vsq = fabs(vsq_calc(x0,Bsq,QdotBsq,Qtsq,Qdotn,D)) ; // guaranteed to be positive \n", "388a403\n", "> vsq = fabs(vsq_calc(x0,Bsq,QdotBsq,Qtsq,Qdotn,D)) ; // guaranteed to be positive\n", "390c405,406\n", "< return( ( vsq > 1. ) ? (1.0 - dv) : vsq ); \n", "---\n", "> \n", "> return( ( vsq > 1. ) ? (1.0 - dv) : vsq );\n", "393a410\n", "> \n", "396,398c413,415\n", "< validate_x(): \n", "< \n", "< -- makes sure that x[0,1] have physical values, based upon \n", "---\n", "> validate_x():\n", "> \n", "> -- makes sure that x[0,1] have physical values, based upon\n", "400c417\n", "< \n", "---\n", "> \n", "403c420\n", "< static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] ) \n", "---\n", "> static void validate_x(CCTK_REAL x[2], CCTK_REAL x0[2] )\n", "405c422\n", "< \n", "---\n", "> \n", "408c425\n", "< /* Always take the absolute value of x[0] and check to see if it's too big: */ \n", "---\n", "> /* Always take the absolute value of x[0] and check to see if it's too big: */\n", "411c428\n", "< \n", "---\n", "> \n", "419a437\n", "> \n", "422c440\n", "< general_newton_raphson(): \n", "---\n", "> general_newton_raphson():\n", "429,431c447,449\n", "< static int general_newton_raphson( CCTK_REAL x[], int n, long &n_iter,\n", "< void (*funcd) (CCTK_REAL [], CCTK_REAL [], CCTK_REAL [], \n", "< CCTK_REAL [][NEWT_DIM], CCTK_REAL *, \n", "---\n", "> static int general_newton_raphson( eos_struct eos, CCTK_REAL x[], int n, long &n_iter,\n", "> void (*funcd) (eos_struct, CCTK_REAL [], CCTK_REAL [], CCTK_REAL [],\n", "> CCTK_REAL [][NEWT_DIM], CCTK_REAL *,\n", "443c461\n", "< errx = 1. ; \n", "---\n", "> errx = 1. ;\n", "452c470,472\n", "< while( keep_iterating ) { \n", "---\n", "> while( keep_iterating ) {\n", "> \n", "> (*funcd) (eos, x, dx, resid, jac, &f, &df, n, Bsq,QdotBsq,Qtsq,Qdotn,D); /* returns with new dx, f, df */\n", "454,455d473\n", "< (*funcd) (x, dx, resid, jac, &f, &df, n, Bsq,QdotBsq,Qtsq,Qdotn,D); /* returns with new dx, f, df */\n", "< \n", "484c502\n", "< \n", "---\n", "> \n", "491c509\n", "< if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0)) \n", "---\n", "> if( ((fabs(errx) <= NEWT_TOL)&&(doing_extra == 0))\n", "509c527\n", "< } \n", "---\n", "> }\n", "522a541\n", "> \n", "525c544\n", "< func_vsq(): \n", "---\n", "> func_vsq():\n", "540c559\n", "< static void func_vsq(CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[], \n", "---\n", "> static void func_vsq(eos_struct eos, CCTK_REAL x[], CCTK_REAL dx[], CCTK_REAL resid[],\n", "545c564\n", "< \n", "---\n", "> \n", "568c587\n", "< \n", "---\n", "> \n", "579c598\n", "< p = pressure_rho0_w(rho0,w) ;\n", "---\n", "> p = pressure_rho0_w(eos, rho0,w) ;\n", "586c605\n", "< p = pressure_rho0_w(rho0,w) ;\n", "---\n", "> p = pressure_rho0_w(eos, rho0,w) ;\n", "594,595c613,614\n", "< \n", "< p_tmp = pressure_W_vsq( W, vsq , D);\n", "---\n", "> \n", "> p_tmp = pressure_W_vsq( eos, W, vsq , D);\n", "597c616\n", "< dPdvsq = dpdvsq_calc( W, vsq, D );\n", "---\n", "> dPdvsq = dpdvsq_calc( eos, W, vsq, D );\n", "599,601c618,620\n", "< // These expressions were calculated using Mathematica, but made into efficient \n", "< // code using Maple. Since we know the analytic form of the equations, we can \n", "< // explicitly calculate the Newton-Raphson step: \n", "---\n", "> // These expressions were calculated using Mathematica, but made into efficient\n", "> // code using Maple. Since we know the analytic form of the equations, we can\n", "> // explicitly calculate the Newton-Raphson step:\n", "627c646\n", "< \n", "---\n", "> \n", "636,642c655,662\n", "< /********************************************************************** \n", "< ********************************************************************** \n", "< \n", "< The following routines specify the equation of state. All routines \n", "< above here should be indpendent of EOS. If the user wishes \n", "< to use another equation of state, the below functions must be replaced \n", "< by equivalent routines based upon the new EOS. \n", "---\n", "> \n", "> /**********************************************************************\n", "> **********************************************************************\n", "> \n", "> The following routines specify the equation of state. All routines\n", "> above here should be indpendent of EOS. If the user wishes\n", "> to use another equation of state, the below functions must be replaced\n", "> by equivalent routines based upon the new EOS.\n", "646a667\n", "> \n", "648,652c669,673\n", "< /********************************************************************** \n", "< pressure_W_vsq(): \n", "< \n", "< -- Gamma-law equation of state;\n", "< -- pressure as a function of W, vsq, and D:\n", "---\n", "> /**********************************************************************\n", "> pressure_W_vsq():\n", "> \n", "> -- Hybrid single and piecewise polytropic equation of state;\n", "> -- pressure as a function of P_cold, eps_cold, W, vsq, and D:\n", "654c675\n", "< static CCTK_REAL pressure_W_vsq(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D) \n", "---\n", "> static CCTK_REAL pressure_W_vsq(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D)\n", "655a677,678\n", "> \n", "> #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", "656a680,691\n", "> #endif\n", "> \n", "> // Compute gamma^{-2} = 1 - v^{2} and gamma^{-1}\n", "> CCTK_REAL inv_gammasq = 1.0 - vsq;\n", "> CCTK_REAL inv_gamma = sqrt(inv_gammasq);\n", "> \n", "> // Compute rho_b = D / gamma\n", "> CCTK_REAL rho_b = D*inv_gamma;\n", "> \n", "> // Compute P_cold and eps_cold\n", "> CCTK_REAL P_cold, eps_cold;\n", "> compute_P_cold__eps_cold(eos,rho_b, P_cold,eps_cold);\n", "658,661c693,694\n", "< CCTK_REAL gtmp;\n", "< gtmp = 1. - vsq;\n", "< \n", "< return( (gamma_th /* <- Should be local polytropic Gamma factor */ - 1.) * ( W * gtmp - D * sqrt(gtmp) ) / gamma_th /* <- Should be local polytropic Gamma factor */ );\n", "---\n", "> // Compute p = P_{cold} + P_{th}\n", "> return( ( P_cold + (Gamma_th - 1.0)*( W*inv_gammasq - D*inv_gamma*( 1.0 + eps_cold ) ) )/Gamma_th );\n", "665a699\n", "> \n", "667,669c701,703\n", "< /********************************************************************** \n", "< dpdW_calc_vsq(): \n", "< \n", "---\n", "> /**********************************************************************\n", "> dpdW_calc_vsq():\n", "> \n", "673a708,709\n", "> \n", "> #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", "675c711,713\n", "< return( (gamma_th /* <- Should be local polytropic Gamma factor */ - 1.) * (1. - vsq) / gamma_th /* <- Should be local polytropic Gamma factor */ ) ;\n", "---\n", "> #endif\n", "> \n", "> return( (Gamma_th - 1.0) * (1.0 - vsq) / Gamma_th ) ;\n", "678a717\n", "> \n", "680,682c719,721\n", "< /********************************************************************** \n", "< dpdvsq_calc(): \n", "< \n", "---\n", "> /**********************************************************************\n", "> dpdvsq_calc():\n", "> \n", "685c724\n", "< static CCTK_REAL dpdvsq_calc(CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D)\n", "---\n", "> static CCTK_REAL dpdvsq_calc(eos_struct eos, CCTK_REAL W, CCTK_REAL vsq, CCTK_REAL &D)\n", "686a726,728\n", "> \n", "> // This sets Gamma_th\n", "> #ifndef ENABLE_STANDALONE_IGM_C2P_SOLVER\n", "688c730,772\n", "< return( (gamma_th /* <- Should be local polytropic Gamma factor */ - 1.) * ( 0.5 * D / sqrt(1.-vsq) - W ) / gamma_th /* <- Should be local polytropic Gamma factor */ ) ;\n", "---\n", "> #endif\n", "> \n", "> \n", "> // Set gamma and rho\n", "> CCTK_REAL gamma = 1.0/sqrt(1.0 - vsq);\n", "> CCTK_REAL rho_b = D/gamma;\n", "> \n", "> // Compute P_cold and eps_cold\n", "> CCTK_REAL P_cold, eps_cold;\n", "> compute_P_cold__eps_cold(eos,rho_b, P_cold,eps_cold);\n", "> \n", "> // Set basic polytropic quantities\n", "> int polytropic_index = find_polytropic_K_and_Gamma_index(eos,rho_b);\n", "> CCTK_REAL Gamma_ppoly_tab = eos.Gamma_ppoly_tab[polytropic_index];\n", "> \n", "> \n", "> /* Now we implement the derivative of P_cold with respect\n", "> * to v^{2}, given by\n", "> * ----------------------------------------------------\n", "> * | dP_cold/dvsq = gamma^{2 + Gamma_{poly}/2} P_{cold} |\n", "> * ----------------------------------------------------\n", "> */\n", "> CCTK_REAL dPcold_dvsq = P_cold * pow(gamma,2.0 + 0.5*Gamma_ppoly_tab);\n", "> \n", "> \n", "> /* Now we implement the derivative of eps_cold with respect\n", "> * to v^{2}, given by\n", "> * -----------------------------------------------------------------------------------\n", "> * | deps_cold/dvsq = gamma/(D*(Gamma_ppoly_tab-1)) * (dP_cold/dvsq + gamma^{2} P_cold / 2) |\n", "> * -----------------------------------------------------------------------------------\n", "> */\n", "> CCTK_REAL depscold_dvsq = ( gamma/(D*(Gamma_ppoly_tab-1.0)) ) * ( dPcold_dvsq + 0.5*gamma*gamma*P_cold );\n", "> \n", "> /* Now we implement the derivative of p_hybrid with respect\n", "> * to v^{2}, given by\n", "> * -----------------------------------------------------------------------------\n", "> * | dp/dvsq = Gamma_th^{-1}( dP_cold/dvsq |\n", "> * | + (Gamma_{th}-1)*(-W |\n", "> * | + D gamma (1 + eps_cold)/2 |\n", "> * | - (D/gamma) * deps_cold/dvsq) ) |\n", "> * -----------------------------------------------------------------------------\n", "> */\n", "> return( ( dPcold_dvsq + (Gamma_th-1.0)*( -W + D*gamma*(1+eps_cold)/2.0 - D*depscold_dvsq/gamma ) )/Gamma_th );\n", "692c776\n", "< /****************************************************************************** \n", "---\n", "> /******************************************************************************\n", "694a779\n", "> #endif\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_utoprim_2d.c\"\n", "original_IGM_file_name = \"harm_utoprim_2d-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_utoprim_2d__c = !diff $original_IGM_file_path $outfile_path__harm_utoprim_2d__c\n", "\n", "if Validation__harm_utoprim_2d__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_utoprim_2d.c: PASSED!\")\n", "else:\n", " # If the validation fails, we keep the original IGM source code file\n", " print(\"Validation test for harm_utoprim_2d.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_utoprim_2d__c:\n", " print(diff_line)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 5: 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__harm_utoprim_2d.pdf](Tutorial-IllinoisGRMHD__harm_utoprim_2d.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": 24, "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__harm_utoprim_2d.ipynb\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_utoprim_2d.tex\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_utoprim_2d.tex\n", "#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__harm_utoprim_2d.tex\n", "!rm -f Tut*.out Tut*.aux Tut*.log" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.1" } }, "nbformat": 4, "nbformat_minor": 4 }