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