{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
"<script>\n",
" window.dataLayer = window.dataLayer || [];\n",
" function gtag(){dataLayer.push(arguments);}\n",
" gtag('js', new Date());\n",
"\n",
" gtag('config', 'UA-59152712-8');\n",
"</script>\n",
"\n",
"# Tutorial-IllinoisGRMHD: mhdflux.C\n",
"\n",
"## Authors: Leo Werneck & Zach Etienne\n",
"\n",
"<font color='red'>**This module is currently under development**</font>\n",
"\n",
"## In this tutorial module we explain how the MHD fluxes are evaluated within IllinoisGRMHD.\n",
"\n",
"### Required and recommended citations:\n",
"\n",
"* **(Required)** Etienne, Z. B., Paschalidis, V., Haas R., Mösta P., and Shapiro, S. L. IllinoisGRMHD: an open-source, user-friendly GRMHD code for dynamical spacetimes. Class. Quantum Grav. 32 (2015) 175009. ([arxiv:1501.07276](http://arxiv.org/abs/1501.07276)).\n",
"* **(Required)** Noble, S. C., Gammie, C. F., McKinney, J. C., Del Zanna, L. Primitive Variable Solvers for Conservative General Relativistic Magnetohydrodynamics. Astrophysical Journal, 641, 626 (2006) ([astro-ph/0512420](https://arxiv.org/abs/astro-ph/0512420)).\n",
"* **(Recommended)** Del Zanna, L., Bucciantini N., Londrillo, P. An efficient shock-capturing central-type scheme for multidimensional relativistic flows - II. Magnetohydrodynamics. A&A 400 (2) 397-413 (2003). DOI: 10.1051/0004-6361:20021641 ([astro-ph/0210618](https://arxiv.org/abs/astro-ph/0210618))."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='toc'></a>\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](#basic_quantities): **Basic quantities**\n",
"1. [Step 3](#eos_quantities): **Equation of State (EOS) quantities**\n",
"1. [Step 4](#impose_speed_limit_output_u0): **Impose a speed limit and determine $u^{0}$**\n",
"1. [Step 5](#magnetic_field_comoving_frame): **The magnetic field measured in the comoving fluid frame, $b^{\\mu}$**\n",
"1. [Step 6](#min_max_characteristic_speeds): **The minimum and maximum characteristic speeds**\n",
"1. [Step 7](#rho_star_flux): **MHD flux: $\\rho_{\\star}$**\n",
"1. [Step 8](#energy_flux): **MHD flux: $\\tilde{\\tau}$**\n",
"1. [Step 9](#momentum_flux): **MHD flux: $\\tilde{S}_{i}$**\n",
"1. [Step 10](#code_validation): **Code validation**\n",
"1. [Step 11](#latex_pdf_output): **Output this notebook to $\\LaTeX$-formatted PDF file**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='src_dir'></a>\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__mhdflux__C = os.path.join(IGM_src_dir_path,\"mhdflux.C\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='introduction'></a>\n",
"\n",
"# Step 1: Introduction \\[Back to [top](#toc)\\]\n",
"$$\\label{introduction}$$\n",
"\n",
"In this tutorial module we will compute the flux terms for the conservative variables $U=\\left\\{\\rho_{\\star},\\tilde{\\tau},\\tilde{S}_{i}\\right\\}$, which are defined in terms of the primitive variables as\n",
"\n",
"$$\n",
"\\left(\n",
"\\begin{matrix}\n",
"\\rho_{\\star}\\\\\n",
"\\tilde{\\tau}\\\\\n",
"\\tilde{S}_{i}\n",
"\\end{matrix}\n",
"\\right)\n",
"=\n",
"\\left(\n",
"\\begin{matrix}\n",
"\\alpha\\sqrt{\\gamma}\\rho_{b}u^{0}\\\\\n",
"\\alpha^{2}\\sqrt{\\gamma}T^{00}-\\rho_{\\star}\\\\\n",
"\\left(\\rho_{\\star}h + \\alpha u^{0}\\sqrt{\\gamma}b^{2}\\right)u_{i} - \\alpha\\sqrt{\\gamma}b^{0}b_{i}\n",
"\\end{matrix}\n",
"\\right)\\ .\n",
"$$\n",
"\n",
"The flux terms for these conservative variables are\n",
"\n",
"$$\n",
"\\boldsymbol{F} \n",
"= \n",
"\\left(\n",
"\\begin{matrix}\n",
"F^{i}_{\\rho_{\\star}}\\\\\n",
"F^{i}_{\\tilde{\\tau}}\\\\\n",
"\\left(F_{\\tilde{S}}\\right)^{j}_{\\ i}\n",
"\\end{matrix}\n",
"\\right)\n",
"=\n",
"\\left(\n",
"\\begin{matrix}\n",
"\\rho_{\\star}v^{i}\\\\\n",
"\\alpha^{2}\\sqrt{\\gamma}T^{0j} - \\rho_{\\star}v^{j}\\\\\n",
"\\alpha\\sqrt{\\gamma}T^{j}_{\\ i}\n",
"\\end{matrix}\n",
"\\right)\\ .\n",
"$$\n",
"\n",
"The MHD flux algorithm computes, for each of the fluxes above, the standard Harten-Lax-van Leer (HLL) flux\n",
"\n",
"$$\n",
"\\boxed{F^{\\rm HLL} = \\frac{c^{-}F_{r} + c^{+}F_{l} - c^{+}c^{-}\\left(U_{r} - U_{l}\\right)}{c^{+} + c^{-}}}\\ .\n",
"$$\n",
"\n",
"As we go through the implementation, we will notice that will need a lot of the functions defined within the `inlined_functions.C` code file of `IllinoisGRMHD`. We will therefore explain the algorithm of each function in quite some detail, so that the reader can go through this tutorial without having to stop and read the [`inlined_functions.C` tutorial module](Tutorial-IllinoisGRMHD_inlined_functions.ipynb). We do encourage the reader to go through the module as well, though, since we will be presenting the math behind each function, but not the code itself."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='basic_quantities'></a>\n",
"\n",
"# Step 2: Basic quantities \\[Back to [top](#toc)\\]\n",
"$$\\label{basic_quantities}$$\n",
"\n",
"We start by defining useful quantities, such as $\\psi^{\\pm4}$, $\\psi^{6}$, $\\alpha\\sqrt{\\gamma}$, $\\alpha^{-1}$, and $\\alpha^{-2}$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Writing ../src/mhdflux.C\n"
]
}
],
"source": [
"%%writefile $outfile_path__mhdflux__C\n",
"//-----------------------------------------------------------------------------\n",
"// Compute the flux for advecting rho_star, tau (Font's energy variable),\n",
"// and S_i .\n",
"//-----------------------------------------------------------------------------\n",
"static inline void mhdflux(int i,int j,int k,const int flux_dirn,CCTK_REAL *Ul,CCTK_REAL *Ur, CCTK_REAL *FACEVAL,CCTK_REAL *FACEVAL_LAPSE_PSI4,eos_struct &eos, CCTK_REAL Gamma_th,\n",
" CCTK_REAL &cmax,CCTK_REAL &cmin,\n",
" CCTK_REAL &rho_star_flux,CCTK_REAL &tau_flux,CCTK_REAL &st_x_flux,CCTK_REAL &st_y_flux,CCTK_REAL &st_z_flux) {\n",
"\n",
" CCTK_REAL psi4 = FACEVAL_LAPSE_PSI4[PSI4];\n",
" CCTK_REAL psi6 = FACEVAL_LAPSE_PSI4[PSI4]*FACEVAL_LAPSE_PSI4[PSI2];\n",
" CCTK_REAL psim4 = 1.0/(psi4);\n",
"\n",
" CCTK_REAL alpha_sqrt_gamma = FACEVAL_LAPSE_PSI4[LAPSE]*psi6;\n",
" CCTK_REAL ONE_OVER_LAPSE = 1.0/FACEVAL_LAPSE_PSI4[LAPSE];\n",
" CCTK_REAL ONE_OVER_LAPSE_SQUARED=SQR(ONE_OVER_LAPSE);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='eos_quantities'></a>\n",
"\n",
"# Step 3: Equation of State (EOS) quantities \\[Back to [top](#toc)\\]\n",
"$$\\label{eos_quantities}$$\n",
"\n",
"Next we compute the quantities $P_{\\rm cold}$, $\\epsilon_{\\rm cold}$, $\\frac{dP_{\\rm cold}}{d\\rho}$, $\\epsilon_{\\rm th}$, $h$, and $\\Gamma_{\\rm cold}$. We have written $\\rho_{b}\\equiv\\rho$ so that the discussion constains a notation that is slightly less cluttered with indices.\n",
"\n",
"The function `compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold()` is defined within the `inlined_functions.C` code file of `IllinoisGRMHD`. It currently implementes 3 different cases:\n",
"\n",
"1. **Case 1: $\\boldsymbol{\\rho_{b} = 0}$**\n",
"\n",
" In this case, we set $\\Gamma_{\\rm cold} = \\Gamma_{\\rm tab}$, i.e. it's tabulated input value, and all other quantities to zero: $P_{\\rm cold}=\\epsilon_{\\rm cold}=\\frac{dP_{\\rm cold}}{d\\rho}=h=\\epsilon_{\\rm th}=0$\n",
" \n",
"2. **Case 2: Single Polytrope EOS**\n",
"\n",
" In this case we have have the relations:\n",
" \n",
" $$\n",
" \\boxed{\n",
" \\begin{align}\n",
" P_{\\rm cold} &= \\kappa \\rho^{\\Gamma_{\\rm cold}}\\ ,\\\\\n",
" \\epsilon_{\\rm cold} &= \\left.\\left(\\frac{P_{\\rm cold}}{\\rho}\\right)\\middle/\\left(\\Gamma_{\\rm cold}-1\\right)\\right.\\ ,\\\\\n",
" \\frac{dP_{\\rm cold}}{d\\rho} &= \\frac{\\Gamma_{\\rm cold}P_{\\rm cold}}{\\rho}\\ ,\\\\\n",
" \\epsilon_{\\rm th} &= \\left.\\left(\\frac{P-P_{\\rm cold}}{\\rho}\\right)\\middle/\\left(\\Gamma_{\\rm cold}-1\\right)\\right.\\ ,\\\\\n",
" h &= 1+\\epsilon_{\\rm cold}+\\epsilon_{\\rm th} + \\frac{P}{\\rho}\\ ,\\\\\n",
" \\Gamma_{\\rm cold} &= \\Gamma_{\\rm tab}\\ .\n",
" \\end{align}}\n",
" $$\n",
" \n",
"2. **Case 3: Piecewise Polytrope EOS**\n",
"\n",
" We now consider a set of dividing densities $\\rho_{\\min} < \\rho_{1} < \\cdots < \\rho_{n-1} < \\rho_{\\max}$ such that the pressure and energy density are everywhere continuous. In each interval, the EOS is a single polytrope, with its own $K_{i}$ and $\\left(\\Gamma_{cold}\\right)_{i}\\equiv\\Gamma_{i}$, so that\n",
" \n",
" $$\n",
" \\boxed{\n",
" P_{\\rm cold} =\n",
" \\left\\{\n",
" \\begin{matrix}\n",
" K_{0}\\rho^{\\Gamma_{0}} & , & \\rho \\leq \\rho_{0}\\\\\n",
" K_{1}\\rho^{\\Gamma_{1}} & , & \\rho_{0} \\leq \\rho \\leq \\rho_{1}\\\\\n",
" \\vdots & & \\vdots\\\\\n",
" K_{j}\\rho^{\\Gamma_{j}} & , & \\rho_{j-1} \\leq \\rho \\leq \\rho_{j}\\\\\n",
" \\vdots & & \\vdots\\\\\n",
" K_{N-1}\\rho^{\\Gamma_{N-1}} & , & \\rho_{N-2} \\leq \\rho \\leq \\rho_{N-1}\\\\\n",
" K_{N}\\rho^{\\Gamma_{N}} & , & \\rho \\geq \\rho_{N-1}\n",
" \\end{matrix}\n",
" \\right.\n",
" }\\ .\n",
" $$\n",
" \n",
" Then, for each set $\\left\\{P_{i},K_{i},\\Gamma_{i}\\right\\}$ the desired quantities are computed using the same equations used in Case 2."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/mhdflux.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__mhdflux__C\n",
"\n",
"\n",
" // First compute P_{cold}, \\epsilon_{cold}, dP_{cold}/drho, \\epsilon_{th}, h, and \\Gamma_{cold},\n",
" // for right and left faces:\n",
" CCTK_REAL P_coldr,eps_coldr,dPcold_drhor=0,eps_thr=0,h_r=0,Gamma_coldr;\n",
" compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(Ur,eos,Gamma_th,P_coldr,eps_coldr,dPcold_drhor,eps_thr,h_r,Gamma_coldr);\n",
" CCTK_REAL P_coldl,eps_coldl,dPcold_drhol=0,eps_thl=0,h_l=0,Gamma_coldl;\n",
" compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(Ul,eos,Gamma_th,P_coldl,eps_coldl,dPcold_drhol,eps_thl,h_l,Gamma_coldl);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='impose_speed_limit_output_u0'></a>\n",
"\n",
"# Step 4: Impose a speed limit and determine $u^{0}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{impose_speed_limit_output_u0}$$\n",
"\n",
"We now call upon the `impose_speed_limit_output_u0()` function inside the `inlined_functions.C` code file of `IllinoisGRMHD`. The basic algorithm performed by this function is summarized here. We start by evaluating the quantity\n",
"\n",
"$$\n",
"\\begin{align}\n",
"{\\rm one\\_minus\\_one\\_over\\_alpha\\_u0\\_squared} \\equiv A \n",
"&= \\gamma_{ij}\\left(\\frac{v^{i}+\\beta^{i}}{\\alpha}\\right)\\left(\\frac{v^{j}+\\beta^{j}}{\\alpha}\\right)\\\\\n",
"&= \\frac{\\gamma_{ij}}{\\alpha^{2}}\\left[\\frac{\\gamma^{ik}u_{k}}{u^{0}} - \\beta^{i} + \\beta^{i}\\right]\\left[\\frac{\\gamma^{j\\ell}u_{\\ell}}{u^{0}} - \\beta^{j} + \\beta^{j}\\right]\\\\\n",
"&=\\frac{\\gamma_{ij}u^{i}u^{j}}{\\left(\\alpha u^{0}\\right)^{2}}\\\\\n",
"&=\\frac{\\left(\\alpha u^{0}\\right)^{2}-1}{\\left(\\alpha u^{0}\\right)^{2}}\\\\\n",
"&=1 - \\frac{1}{\\left(\\alpha u^{0}\\right)^{2}}\\ \\\\\n",
"\\implies \\boxed{A = 1 - \\frac{1}{\\left(\\alpha u^{0}\\right)^{2}}}\\ ,\n",
"\\end{align}\n",
"$$\n",
"\n",
"where when going from line 1 to 2 and from line 3 to 4 we have used eqs. (53) and (56) from [Duez *et al.*](https://arxiv.org/pdf/astro-ph/0503420.pdf), respectively. Then we construct the \"speed limit quantity\"\n",
"\n",
"$$\n",
"{\\rm ONE\\_MINUS\\_ONE\\_OVER\\_GAMMA\\_SPEED\\_LIMIT\\_SQUARED} \\equiv B = 1-\\frac{1}{\\gamma^{2}_{\\rm speed\\ limit}}\\ .\n",
"$$\n",
"\n",
"If $A > B$, then we construct the correction factor $C\\equiv A / B$, and adjust the velocities using\n",
"\n",
"$$\n",
"\\boxed{v^{i} \\to \\left(v^{i}+\\beta^{i}\\right)C - \\beta^{i}}\\ .\n",
"$$\n",
"\n",
"Finally, since $A$ is evaluated using the first line above, namely\n",
"\n",
"$$\n",
"A = \\gamma_{ij}\\left(\\frac{v^{i}+\\beta^{i}}{\\alpha}\\right)\\left(\\frac{v^{j}+\\beta^{j}}{\\alpha}\\right)\\ ,\n",
"$$\n",
"\n",
"but we know, from the last line how $A$ and $u^{0}$ are related, we compute\n",
"\n",
"$$\n",
"\\boxed{u^{0} = \\frac{1}{\\alpha\\sqrt{1-A}}}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/mhdflux.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__mhdflux__C\n",
"\n",
"\n",
" //Compute face velocities\n",
" // Begin by computing u0\n",
" output_stats stats; stats.failure_checker=0;\n",
" CCTK_REAL u0_r,u0_l;\n",
" impose_speed_limit_output_u0(FACEVAL,Ur,psi4,ONE_OVER_LAPSE,stats,u0_r);\n",
" impose_speed_limit_output_u0(FACEVAL,Ul,psi4,ONE_OVER_LAPSE,stats,u0_l);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='magnetic_field_comoving_frame'></a>\n",
"\n",
"# Step 5: The magnetic field measured in the comoving fluid frame, $b^{\\mu}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{magnetic_field_comoving_frame}$$\n",
"\n",
"Now we use th function `compute_smallba_b2_and_u_i_over_u0_psi4()` from the `inlined_functions.C` code file of `IllinoisGRMHD`.\n",
"\n",
"We will need the following identities\n",
"\n",
"$$\n",
"\\begin{align}\n",
"v^{i} &= \\frac{u^{i}}{u^{0}}\\ ,\\\\\n",
"B^{0}_{(u)} &= \\frac{u_{i}B^{i}}{\\alpha}\\ ,\\\\\n",
"B^{i}_{(u)} &= \\frac{1}{u^{0}}\\left(\\frac{B^{i}}{\\alpha} + u^{i}B^{0}_{(u)}\\right)\\ ,\\\\\n",
"b^{\\mu} &= \\frac{B^{\\mu}_{(u)}}{\\sqrt{4\\pi}}\\ .\n",
"\\end{align}\n",
"$$\n",
"\n",
"We start by setting the relation\n",
"\n",
"$$\n",
"b^{0} = \\frac{u_{i}B^{i}}{\\alpha\\sqrt{4\\pi}} \\implies \\boxed{\\alpha\\sqrt{4\\pi}b^{0} = u_{i}B^{i}}\\ .\n",
"$$\n",
"\n",
"Then we compute\n",
"\n",
"$$\n",
"\\begin{align}\n",
"b^{i} &= \\frac{B^{i}_{(u)}}{\\sqrt{4\\pi}}\\\\\n",
" &= \\frac{1}{u^{0}\\sqrt{4\\pi}}\\left(\\frac{B^{i}}{\\alpha} + B^{0}_{(u)}u^{i}\\right)\\\\\n",
" &= \\frac{1}{u^{0}\\sqrt{4\\pi}}\\left(\\frac{B^{i}}{\\alpha} + \\sqrt{4\\pi}b^{0}u^{i}\\right)\\\\\n",
" &= \\frac{1}{\\alpha\\sqrt{4\\pi}}\\left(\\frac{B^{i}}{u^{0}} + \\alpha\\sqrt{4\\pi}b^{0}\\frac{u^{i}}{u^{0}}\\right)\\\\\n",
"\\implies &\\boxed{b^{i} = \\frac{1}{\\alpha\\sqrt{4\\pi}}\\left(\\frac{B^{i}}{u^{0}} + \\alpha\\sqrt{4\\pi}b^{0}v^{i}\\right)}\\ .\n",
"\\end{align}\n",
"$$\n",
"\n",
"Finally, we compute\n",
"\n",
"$$\n",
"\\begin{align}\n",
"b^{2} &= g_{\\mu\\nu}b^{\\mu}b^{\\nu}\\\\\n",
" &= g_{00}\\left(b^{0}\\right)^{2} + g_{ij}b^{i}b^{j} + 2g_{0i}b^{0}b^{i}\\\\\n",
" &= \\left(-\\alpha^{2} + \\gamma_{ij}\\beta^{i}\\beta^{j}\\right)\\left(b^{0}\\right)^{2} + \\gamma_{ij}b^{i}b^{j} + 2b^{0}\\gamma_{ij}\\beta^{j}b^{i}\\\\\n",
" &= -\\left(\\alpha b^{0}\\right)^{2} + \\gamma_{ij}\\left[b^{i}b^{j} + 2b^{0}b^{i}\\beta^{j} + \\left(b^{0}\\right)^{2}\\beta^{i}\\beta^{j}\\right]\\\\\n",
"\\implies &\\boxed{b^{2} = -\\left(\\alpha b^{0}\\right)^{2} + \\gamma_{ij}\\left(b^{i} + b^{0}\\beta^{i}\\right)\\left(b^{j} + b^{0}\\beta^{j}\\right)}\n",
"\\end{align}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/mhdflux.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__mhdflux__C\n",
"\n",
"\n",
" //Next compute b^{\\mu}, the magnetic field measured in the comoving fluid frame:\n",
" CCTK_REAL ONE_OVER_LAPSE_SQRT_4PI = ONE_OVER_LAPSE*ONE_OVER_SQRT_4PI;\n",
" /***********************************************************/\n",
" /********** RIGHT FACE ************/\n",
" // Note that smallbr[4] = b^a defined in Gammie's paper, on the right face.\n",
" CCTK_REAL u_x_over_u0_psi4r,u_y_over_u0_psi4r,u_z_over_u0_psi4r;\n",
" CCTK_REAL smallbr[NUMVARS_SMALLB];\n",
" // Compute b^{a}, b^2, and u_i over u^0\n",
" compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ur,u0_r,ONE_OVER_LAPSE_SQRT_4PI,\n",
" u_x_over_u0_psi4r,u_y_over_u0_psi4r,u_z_over_u0_psi4r,smallbr);\n",
" // Then compute u_xr,u_yr, and u_zr. We need to set the zeroth component so we can specify U_LOWER{r,l}[{UX,UY,UZ}] (UX=1,UY=2,UZ=3).\n",
" CCTK_REAL U_LOWERr[4] = { 0.0, u_x_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4],\n",
" u_z_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4] };\n",
" /********** LEFT FACE ************/\n",
" // Note that smallbl[4] = b^a defined in Gammie's paper, on the left face.\n",
" CCTK_REAL u_x_over_u0_psi4l,u_y_over_u0_psi4l,u_z_over_u0_psi4l;\n",
" CCTK_REAL smallbl[NUMVARS_SMALLB];\n",
" // Compute b^{a}, b^2, and u_i over u^0\n",
" compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ul,u0_l,ONE_OVER_LAPSE_SQRT_4PI,\n",
" u_x_over_u0_psi4l,u_y_over_u0_psi4l,u_z_over_u0_psi4l,smallbl);\n",
" // Then compute u_xr,u_yr, and u_zr. We need to set the zeroth component so we can specify U_LOWER{r,l}[{UX,UY,UZ}]\n",
" CCTK_REAL U_LOWERl[4] = { 0.0, u_x_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4],\n",
" u_z_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4] };\n",
" /***********************************************************/"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='min_max_characteristic_speeds'></a>\n",
"\n",
"# Step 6: The minimum and maximum characteristic speeds \\[Back to [top](#toc)\\]\n",
"$$\\label{min_max_characteristic_speeds}$$\n",
"\n",
"We will now explain two different functions contained in the `inlined_functions.C` code file of `IllinoisGRMHD`. These functions are `compute_v02()` and `find_cp_cm`. These functions are called with the objective of computing the minimum ($-$) and maximum ($+$) characteristic speeds at each cell interface, $c_{\\pm}^{r,l}$.\n",
"\n",
"We approximate the general GRMHD dispersion relation (eq. 27 of [Gammie & McKinney (2003)](https://arxiv.org/pdf/astro-ph/0301509.pdf)) by the simpler expression\n",
"\n",
"$$\n",
"\\omega_{\\rm cm}^{2} = \\left[v_{\\rm A}^{2} + c_{\\rm s}^{2}\\left(1-v_{\\rm A}^{2}\\right)\\right]k_{\\rm cm}^{2}\\ ,\n",
"$$\n",
"\n",
"where $\\omega_{\\rm cm}=-k_{\\mu}u^{\\mu}$ is the frequency and $k_{\\rm cm}^{2} = K_{\\mu}K^{\\mu}$ the wavenumber of an MHD wave mode in the frame comoving with the fluid, where $K_{\\mu}$ is defined as the projection of the wave vector $k^{\\nu}$ onto the direction normal to $u^{\\nu}$: $K_{\\mu} = \\left(g_{\\mu\\nu}+u_{\\mu}u_{\\nu}\\right)k^{\\nu}$. $c_{\\rm s}$ is the sound speed, and $v_{\\rm A}$ is the Alfvén speed, given by\n",
"\n",
"$$\n",
"v_{\\rm A} = \\sqrt{\\frac{b^{2}}{\\rho_{b}h + b^{2}}}\\ .\n",
"$$\n",
"\n",
"With these definitions, we may then solve the approximate dispersion relation above along direction $i$, noting that in the comoving frame $k_{\\mu} = \\left(-\\omega,k_{j}\\delta^{j}_{\\ i}\\right)$ and the wave (phase) velocity is $c_{\\pm} = \\left.\\omega\\middle/\\left(k_{j}\\delta^{j}_{\\ i}\\right)\\right.$. The dispersion can then be written as a quadratic equation for $c_{\\pm}$:\n",
"\n",
"$$\n",
"ac_{\\pm}^{2} + bc_{\\pm} + c = 0\\ ,\n",
"$$\n",
"\n",
"with\n",
"\n",
"$$\n",
"\\boxed{\n",
"\\begin{align}\n",
"a &= \\left(1-v_{0}^{2}\\right)\\left(u^{0}\\right)^{2} - v_{0}^{2}g^{00}\\ ,\\\\\n",
"b &= 2v_{0}^{2}g^{i0} - 2u^{i}u^{0}\\left(1-v^{2}_{0}\\right)\\ ,\\\\\n",
"c &= \\left(1-v_{0}^{2}\\right)\\left(u^{i}\\right)^{2} - v_{0}^{2}g^{ii}\\ ,\\\\\n",
"v_{0}^{2} &= v_{\\rm A}^{2} + c_{\\rm s}^{2}\\left(1-v_{\\rm A}^{2}\\right)\\ ,\\\\\n",
"c_{\\rm s} &= \\left.\\left[\\frac{dP_{\\rm cold}}{d\\rho_{b}} + \\Gamma_{\\rm th}\\left(\\Gamma_{\\rm th}-1\\right)\\epsilon_{\\rm th}\\right]\\middle/h\\right.\\ ,\\\\\n",
"c_{+} &= \\max\\left(\\frac{-b \\pm \\sqrt{b^{2}-4ac}}{2a}\\right)\\ ,\\\\\n",
"c_{-} &= \\min\\left(\\frac{-b \\pm \\sqrt{b^{2}-4ac}}{2a}\\right)\\ .\n",
"\\end{align}\n",
"}\n",
"$$\n",
"\n",
"Finally, after the maximum and minimum speeds $c_{\\pm}$ have been computed at left and right facs, the standard Harten-Lax-van Leer (HLL), approximate Riemann solve is applied to compute fluxes for the three conservative variables $U = \\left\\{\\rho_{\\star},\\tilde{\\tau},\\tilde{S}_{i}\\right\\}$:\n",
"\n",
"$$\n",
"F^{\\rm HLL} = \\frac{c^{-}F_{r} + c^{+}F_{l} - c^{+}c^{-}\\left(U_{r} - U_{l}\\right)}{c^{+} + c^{-}}\\ ,\n",
"$$\n",
"\n",
"where\n",
"\n",
"$$\n",
"\\boxed{c^{\\pm} = \\pm\\max\\left(0,c^{r}_{\\pm},c^{l}_{\\pm}\\right)}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/mhdflux.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__mhdflux__C\n",
"\n",
"\n",
" // Compute v02 = v_A^2 + c_s^2*(1.0-v_A^2), where c_s = sound speed, and v_A = Alfven velocity\n",
" CCTK_REAL v02r,v02l;\n",
" // First right face\n",
" compute_v02(dPcold_drhor,Gamma_th,eps_thr,h_r,smallbr,Ur,v02r);\n",
" // Then left face.\n",
" compute_v02(dPcold_drhol,Gamma_th,eps_thl,h_l,smallbl,Ul,v02l);\n",
"\n",
" int offset=flux_dirn-1;\n",
"\n",
" // Now that we have computed v02 = v_A^2 + c_s^2*(1.0-v_A^2), we can\n",
" // next compute c_+ and c_- using a simplified dispersion relation.\n",
" // Note that, in solving the simplified disp. relation, we overestimate\n",
" // c_+ and c_- by around a factor of 2, making the MHD evolution more\n",
" // diffusive (and potentially more *stable*) than it could be.\n",
" CCTK_REAL cplusr,cminusr,cplusl,cminusl;\n",
" find_cp_cm(cplusr,cminusr,v02r,u0_r,\n",
" Ur[VX+offset],ONE_OVER_LAPSE_SQUARED,FACEVAL[SHIFTX+offset],psim4,FACEVAL[GUPXX+offset]);\n",
" find_cp_cm(cplusl,cminusl,v02l,u0_l,\n",
" Ul[VX+offset],ONE_OVER_LAPSE_SQUARED,FACEVAL[SHIFTX+offset],psim4,FACEVAL[GUPXX+offset]);\n",
"\n",
" // Then compute cmax, cmin. This is required for the HLL flux.\n",
" CCTK_REAL cmaxL = MAX(0.0,MAX(cplusl,cplusr));\n",
" CCTK_REAL cminL = -MIN(0.0,MIN(cminusl,cminusr));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='rho_star_flux'></a>\n",
"\n",
"# Step 7: MHD flux: $\\rho_{\\star}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{rho_star_flux}$$\n",
"\n",
"We now evaluate the flux for $U=\\rho_{\\star}$. First, remember that\n",
"\n",
"$$\n",
"\\rho_{\\star} = \\alpha\\sqrt{\\gamma}\\rho_{b}u^{0}\\ .\n",
"$$\n",
"\n",
"$$\n",
"F^{i}_{\\rho_{\\star}} = \\rho_{\\star} v^{i}\\ ,\n",
"$$\n",
"\n",
"where $i$ is the current flux direction. We can then evaluate the HLL flux in the $i$-direction,\n",
"\n",
"$$\n",
"\\boxed{F^{\\rm HLL} = \\frac{c^{-}F_{r} + c^{+}F_{l} - c^{+}c^{-}\\left(U_{r} - U_{l}\\right)}{c^{+} + c^{-}}}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/mhdflux.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__mhdflux__C\n",
"\n",
"\n",
" //*********************************************************************\n",
" // density flux = \\rho_* v^m, where m is the current flux direction (the m index)\n",
" //*********************************************************************\n",
" CCTK_REAL rho_star_r = alpha_sqrt_gamma*Ur[RHOB]*u0_r;\n",
" CCTK_REAL rho_star_l = alpha_sqrt_gamma*Ul[RHOB]*u0_l;\n",
" CCTK_REAL Fr = rho_star_r*Ur[VX+offset]; // flux_dirn = 2, so offset = 1, implies Ur[VX] -> Ur[VY]\n",
" CCTK_REAL Fl = rho_star_l*Ul[VX+offset]; // flux_dirn = 2, so offset = 1, implies Ul[VX] -> Ul[VY]\n",
"\n",
" // HLL step for rho_star:\n",
" rho_star_flux = (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(rho_star_r-rho_star_l) )/(cmaxL + cminL);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='energy_flux'></a>\n",
"\n",
"# Step 8: MHD flux: $\\tilde{\\tau}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{energy_flux}$$\n",
"\n",
"We then evaluate the flux for $U=\\tilde{\\tau}$. First let us remember that\n",
"\n",
"$$\n",
"\\tilde{\\tau} = \\alpha^{2}\\sqrt{\\gamma}T^{00} - \\rho_{\\star}\\ .\n",
"$$\n",
"\n",
"We then have the flux term\n",
"\n",
"$$\n",
"F^{i}_{\\tilde{\\tau}} = \\alpha^{2}\\sqrt{\\gamma}T^{0j} - \\rho_{\\star}v^{j}\\ ,\n",
"$$\n",
"\n",
"where $i$ is the current flux direction. We can then evaluate the HLL flux in the $i$-direction,\n",
"\n",
"$$\n",
"\\boxed{F^{\\rm HLL} = \\frac{c^{-}F_{r} + c^{+}F_{l} - c^{+}c^{-}\\left(U_{r} - U_{l}\\right)}{c^{+} + c^{-}}}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/mhdflux.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__mhdflux__C\n",
"\n",
"\n",
" //*********************************************************************\n",
" // energy flux = \\alpha^2 \\sqrt{\\gamma} T^{0m} - \\rho_* v^m, where m is the current flux direction (the m index)\n",
" //*********************************************************************\n",
" // First compute some useful metric quantities:\n",
" CCTK_REAL alpha_squared_sqrt_gamma = FACEVAL_LAPSE_PSI4[LAPSE]*alpha_sqrt_gamma;\n",
" CCTK_REAL g4uptm = ONE_OVER_LAPSE_SQUARED*FACEVAL[SHIFTX+offset];\n",
" CCTK_REAL g4uptt = -ONE_OVER_LAPSE_SQUARED;\n",
" /********** RIGHT FACE ************/\n",
" // Compute a couple useful hydro quantities:\n",
" CCTK_REAL rho0_h_plus_b2_r = (Ur[RHOB]*h_r + smallbr[SMALLB2]);\n",
" CCTK_REAL P_plus_half_b2_r = (Ur[PRESSURE]+0.5*smallbr[SMALLB2]);\n",
" // Then compute T^{0m} and the flux:\n",
" CCTK_REAL TUP0m_r = rho0_h_plus_b2_r*SQR(u0_r)*Ur[VX+offset] + P_plus_half_b2_r*g4uptm - smallbr[SMALLBT]*smallbr[SMALLBX+offset];\n",
" Fr = alpha_squared_sqrt_gamma * TUP0m_r - rho_star_r * Ur[VX+offset];\n",
" // Finally compute tau\n",
" CCTK_REAL TUP00_r = rho0_h_plus_b2_r*u0_r*u0_r + P_plus_half_b2_r*g4uptt - smallbr[SMALLBT]*smallbr[SMALLBT];\n",
" CCTK_REAL tau_r = alpha_squared_sqrt_gamma * TUP00_r - rho_star_r;\n",
" /********** LEFT FACE *************/\n",
" // Compute a couple useful hydro quantities:\n",
" CCTK_REAL rho0_h_plus_b2_l = (Ul[RHOB]*h_l + smallbl[SMALLB2]);\n",
" CCTK_REAL P_plus_half_b2_l = (Ul[PRESSURE]+0.5*smallbl[SMALLB2]);\n",
" // Then compute T^{0m} and the flux:\n",
" CCTK_REAL TUP0m_l = rho0_h_plus_b2_l*SQR(u0_l)*Ul[VX+offset] + P_plus_half_b2_l*g4uptm - smallbl[SMALLBT]*smallbl[SMALLBX+offset];\n",
" Fl = alpha_squared_sqrt_gamma * TUP0m_l - rho_star_l * Ul[VX+offset];\n",
" // Finally compute tau\n",
" CCTK_REAL TUP00_l = rho0_h_plus_b2_l*u0_l*u0_l + P_plus_half_b2_l*g4uptt - smallbl[SMALLBT]*smallbl[SMALLBT];\n",
" CCTK_REAL tau_l = alpha_squared_sqrt_gamma * TUP00_l - rho_star_l;\n",
"\n",
" // HLL step for tau:\n",
" tau_flux = (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(tau_r-tau_l) )/(cmaxL + cminL);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='momentum_flux'></a>\n",
"\n",
"# Step 9: MHD flux: $\\tilde{S}_{i}$ \\[Back to [top](#toc)\\]\n",
"$$\\label{momentum_flux}$$\n",
"\n",
"We then evaluate the flux for $U=\\tilde{S}_{i}$. First let us remember that\n",
"\n",
"$$\n",
"\\tilde{S}_{i} = \\left(\\rho_{\\star} h + \\alpha u^{0} \\sqrt{\\gamma}b^{2}\\right)u_{i} - \\alpha\\sqrt{\\gamma}b^{0}b_{i}\\ .\n",
"$$\n",
"\n",
"We then have the flux term\n",
"\n",
"$$\n",
"\\left(F_{\\tilde{S}}\\right)^{j}_{\\ i} = \\alpha\\sqrt{\\gamma}T^{j}_{\\ i}\\ ,\n",
"$$\n",
"\n",
"where $j$ is the current flux direction and $i$ correspond to the component of $\\tilde{S}_{i}$ we are interested in. We can then evaluate the HLL flux in the $i$-direction,\n",
"\n",
"$$\n",
"\\boxed{F^{\\rm HLL} = \\frac{c^{-}F_{r} + c^{+}F_{l} - c^{+}c^{-}\\left(U_{r} - U_{l}\\right)}{c^{+} + c^{-}}}\\ .\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Appending to ../src/mhdflux.C\n"
]
}
],
"source": [
"%%writefile -a $outfile_path__mhdflux__C\n",
"\n",
"\n",
" //*********************************************************************\n",
" // momentum flux = \\alpha \\sqrt{\\gamma} T^m_j, where m is the current flux direction (the m index)\n",
" //*********************************************************************\n",
" // b_j = g_{ij} (b^i + b^t shift^i), g_{ij} = physical metric\n",
" //CCTK_REAL sbtr=0,sbtl=0;\n",
" CCTK_REAL smallb_lowerr[NUMVARS_SMALLB],smallb_lowerl[NUMVARS_SMALLB];\n",
" lower_4vector_output_spatial_part(psi4,FACEVAL,smallbr,smallb_lowerr);\n",
" lower_4vector_output_spatial_part(psi4,FACEVAL,smallbl,smallb_lowerl);\n",
"\n",
" /********** Flux for S_x **********/\n",
" // [S_x flux] = \\alpha \\sqrt{\\gamma} T^m_x, where m is the current flux direction (the m index)\n",
" // Again, offset = 0 for reconstruction in x direction, 1 for y, and 2 for z\n",
" // Note that kronecker_delta[flux_dirn][0] = { 1 if flux_dirn==1, 0 otherwise }.\n",
" Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UX]\n",
" + P_plus_half_b2_r*kronecker_delta[flux_dirn][0] - smallbr[SMALLBX+offset]*smallb_lowerr[SMALLBX] );\n",
" Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UX]\n",
" + P_plus_half_b2_l*kronecker_delta[flux_dirn][0] - smallbl[SMALLBX+offset]*smallb_lowerl[SMALLBX] );\n",
"\n",
" // S_x =\\alpha\\sqrt{\\gamma}( T^0_x )\n",
" CCTK_REAL st_x_r = alpha_sqrt_gamma*( rho0_h_plus_b2_r*u0_r*U_LOWERr[UX] - smallbr[SMALLBT]*smallb_lowerr[SMALLBX] );\n",
" CCTK_REAL st_x_l = alpha_sqrt_gamma*( rho0_h_plus_b2_l*u0_l*U_LOWERl[UX] - smallbl[SMALLBT]*smallb_lowerl[SMALLBX] );\n",
"\n",
" // HLL step for Sx:\n",
" st_x_flux = (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(st_x_r-st_x_l) )/(cmaxL + cminL);\n",
"\n",
" /********** Flux for S_y **********/\n",
" // [S_y flux] = \\alpha \\sqrt{\\gamma} T^m_y, where m is the current flux direction (the m index)\n",
" // Again, offset = 1 for reconstruction in x direction, 2 for y, and 3 for z\n",
" // Note that kronecker_delta[flux_dirn][1] = { 1 if flux_dirn==2, 0 otherwise }.\n",
" Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UY] + P_plus_half_b2_r*kronecker_delta[flux_dirn][1]\n",
" - smallbr[SMALLBX+offset]*smallb_lowerr[SMALLBY] );\n",
" Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UY] + P_plus_half_b2_l*kronecker_delta[flux_dirn][1]\n",
" - smallbl[SMALLBX+offset]*smallb_lowerl[SMALLBY] );\n",
"\n",
" // S_y =\\alpha\\sqrt{\\gamma}( T^0_y )\n",
" CCTK_REAL st_y_r = alpha_sqrt_gamma*( rho0_h_plus_b2_r*u0_r*U_LOWERr[UY] - smallbr[SMALLBT]*smallb_lowerr[SMALLBY] );\n",
" CCTK_REAL st_y_l = alpha_sqrt_gamma*( rho0_h_plus_b2_l*u0_l*U_LOWERl[UY] - smallbl[SMALLBT]*smallb_lowerl[SMALLBY] );\n",
"\n",
" // HLL step for Sy:\n",
" st_y_flux = (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(st_y_r-st_y_l) )/(cmaxL + cminL);\n",
"\n",
" /********** Flux for S_z **********/\n",
" // [S_z flux] = \\alpha \\sqrt{\\gamma} T^m_z, where m is the current flux direction (the m index)\n",
" // Again, offset = 1 for reconstruction in x direction, 2 for y, and 3 for z\n",
" // Note that kronecker_delta[flux_dirn][2] = { 1 if flux_dirn==3, 0 otherwise }.\n",
" Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UZ] + P_plus_half_b2_r*kronecker_delta[flux_dirn][2]\n",
" - smallbr[SMALLBX+offset]*smallb_lowerr[SMALLBZ] );\n",
" Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UZ] + P_plus_half_b2_l*kronecker_delta[flux_dirn][2]\n",
" - smallbl[SMALLBX+offset]*smallb_lowerl[SMALLBZ] );\n",
"\n",
" // S_z =\\alpha\\sqrt{\\gamma}( T^0_z )\n",
" CCTK_REAL st_z_r = alpha_sqrt_gamma*( rho0_h_plus_b2_r*u0_r*U_LOWERr[UZ] - smallbr[SMALLBT]*smallb_lowerr[SMALLBZ] );\n",
" CCTK_REAL st_z_l = alpha_sqrt_gamma*( rho0_h_plus_b2_l*u0_l*U_LOWERl[UZ] - smallbl[SMALLBT]*smallb_lowerl[SMALLBZ] );\n",
"\n",
" // HLL step for Sz:\n",
" st_z_flux = (cminL*Fr + cmaxL*Fl - cminL*cmaxL*(st_z_r-st_z_l) )/(cmaxL + cminL);\n",
"\n",
" cmax = cmaxL;\n",
" cmin = cminL;\n",
"}\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='code_validation'></a>\n",
"\n",
"# Step 10: 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": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Validation test for mhdflux.C: FAILED!\n",
"Diff:\n",
"2c2\n",
"< // Compute the flux for advecting rho_star, tau (Font's energy variable), \n",
"---\n",
"> // Compute the flux for advecting rho_star, tau (Font's energy variable),\n",
"5c5\n",
"< static inline void mhdflux(int i,int j,int k,const int flux_dirn,CCTK_REAL *Ul,CCTK_REAL *Ur, CCTK_REAL *FACEVAL,CCTK_REAL *FACEVAL_LAPSE_PSI4,eos_struct &eos,\n",
"---\n",
"> static inline void mhdflux(int i,int j,int k,const int flux_dirn,CCTK_REAL *Ul,CCTK_REAL *Ur, CCTK_REAL *FACEVAL,CCTK_REAL *FACEVAL_LAPSE_PSI4,eos_struct &eos, CCTK_REAL Gamma_th,\n",
"9d8\n",
"< \n",
"17a17\n",
"> \n",
"20,23c20,24\n",
"< CCTK_REAL P_coldr,eps_coldr,dPcold_drhor=0,eps_thr=0,h_r=0,gamma_coldr;\n",
"< compute_P_cold__eps_cold__dPcold_drho__eps_th__h__gamma_cold(Ur,eos,P_coldr,eps_coldr,dPcold_drhor,eps_thr,h_r,gamma_coldr);\n",
"< CCTK_REAL P_coldl,eps_coldl,dPcold_drhol=0,eps_thl=0,h_l=0,gamma_coldl;\n",
"< compute_P_cold__eps_cold__dPcold_drho__eps_th__h__gamma_cold(Ul,eos,P_coldl,eps_coldl,dPcold_drhol,eps_thl,h_l,gamma_coldl);\n",
"---\n",
"> CCTK_REAL P_coldr,eps_coldr,dPcold_drhor=0,eps_thr=0,h_r=0,Gamma_coldr;\n",
"> compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(Ur,eos,Gamma_th,P_coldr,eps_coldr,dPcold_drhor,eps_thr,h_r,Gamma_coldr);\n",
"> CCTK_REAL P_coldl,eps_coldl,dPcold_drhol=0,eps_thl=0,h_l=0,Gamma_coldl;\n",
"> compute_P_cold__eps_cold__dPcold_drho__eps_th__h__Gamma_cold(Ul,eos,Gamma_th,P_coldl,eps_coldl,dPcold_drhol,eps_thl,h_l,Gamma_coldl);\n",
"> \n",
"31a33\n",
"> \n",
"40c42\n",
"< compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ur,u0_r,ONE_OVER_LAPSE_SQRT_4PI, \n",
"---\n",
"> compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ur,u0_r,ONE_OVER_LAPSE_SQRT_4PI,\n",
"43c45\n",
"< CCTK_REAL U_LOWERr[4] = { 0.0, u_x_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4], \n",
"---\n",
"> CCTK_REAL U_LOWERr[4] = { 0.0, u_x_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4r*u0_r*FACEVAL_LAPSE_PSI4[PSI4],\n",
"50c52\n",
"< compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ul,u0_l,ONE_OVER_LAPSE_SQRT_4PI, \n",
"---\n",
"> compute_smallba_b2_and_u_i_over_u0_psi4(FACEVAL,FACEVAL_LAPSE_PSI4,Ul,u0_l,ONE_OVER_LAPSE_SQRT_4PI,\n",
"53c55\n",
"< CCTK_REAL U_LOWERl[4] = { 0.0, u_x_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4], \n",
"---\n",
"> CCTK_REAL U_LOWERl[4] = { 0.0, u_x_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4], u_y_over_u0_psi4l*u0_l*FACEVAL_LAPSE_PSI4[PSI4],\n",
"56a59\n",
"> \n",
"60c63\n",
"< compute_v02(dPcold_drhor,eos.gamma_th,eps_thr,h_r,smallbr,Ur,v02r);\n",
"---\n",
"> compute_v02(dPcold_drhor,Gamma_th,eps_thr,h_r,smallbr,Ur,v02r);\n",
"62c65\n",
"< compute_v02(dPcold_drhol,eos.gamma_th,eps_thl,h_l,smallbl,Ul,v02l);\n",
"---\n",
"> compute_v02(dPcold_drhol,Gamma_th,eps_thl,h_l,smallbl,Ul,v02l);\n",
"67,68c70,71\n",
"< // next compute c_+ and c_- using a simplified dispersion relation. \n",
"< // Note that, in solving the simplified disp. relation, we overestimate \n",
"---\n",
"> // next compute c_+ and c_- using a simplified dispersion relation.\n",
"> // Note that, in solving the simplified disp. relation, we overestimate\n",
"80a84\n",
"> \n",
"91a96\n",
"> \n",
"122a128\n",
"> \n",
"136c142\n",
"< Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UX] \n",
"---\n",
"> Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UX]\n",
"138c144\n",
"< Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UX] \n",
"---\n",
"> Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UX]\n",
"152c158\n",
"< Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UY] + P_plus_half_b2_r*kronecker_delta[flux_dirn][1] \n",
"---\n",
"> Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UY] + P_plus_half_b2_r*kronecker_delta[flux_dirn][1]\n",
"154c160\n",
"< Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UY] + P_plus_half_b2_l*kronecker_delta[flux_dirn][1] \n",
"---\n",
"> Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UY] + P_plus_half_b2_l*kronecker_delta[flux_dirn][1]\n",
"168c174\n",
"< Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UZ] + P_plus_half_b2_r*kronecker_delta[flux_dirn][2] \n",
"---\n",
"> Fr = alpha_sqrt_gamma*( rho0_h_plus_b2_r*(u0_r*Ur[VX+offset])*U_LOWERr[UZ] + P_plus_half_b2_r*kronecker_delta[flux_dirn][2]\n",
"170c176\n",
"< Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UZ] + P_plus_half_b2_l*kronecker_delta[flux_dirn][2] \n",
"---\n",
"> Fl = alpha_sqrt_gamma*( rho0_h_plus_b2_l*(u0_l*Ul[VX+offset])*U_LOWERl[UZ] + P_plus_half_b2_l*kronecker_delta[flux_dirn][2]\n",
"182a189\n",
"> \n"
]
}
],
"source": [
"# Verify if the code generated by this tutorial module\n",
"# matches the original IllinoisGRMHD source code\n",
"\n",
"# First download the original IllinoisGRMHD source code\n",
"import urllib\n",
"from os import path\n",
"\n",
"original_IGM_file_url = \"https://bitbucket.org/zach_etienne/wvuthorns/raw/5611b2f0b17135538c9d9d17c7da062abe0401b6/IllinoisGRMHD/src/mhdflux.C\"\n",
"original_IGM_file_name = \"mhdflux-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__mhdflux__C = !diff $original_IGM_file_path $outfile_path__mhdflux__C\n",
"\n",
"if Validation__mhdflux__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 mhdflux.C: PASSED!\")\n",
"else:\n",
" # If the validation fails, we keep the original IGM source code file\n",
" print(\"Validation test for mhdflux.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__mhdflux__C:\n",
" print(diff_line)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<a id='latex_pdf_output'></a>\n",
"\n",
"# Step 11: 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__mhdflux.pdf](Tutorial-IllinoisGRMHD__mhdflux.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": 11,
"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__mhdflux.ipynb\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__mhdflux.tex\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__mhdflux.tex\n",
"#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__mhdflux.tex\n",
"!rm -f Tut*.out Tut*.aux Tut*.log"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}