{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# Tutorial-IllinoisGRMHD: mhdflux.C\n", "\n", "## Authors: Leo Werneck & Zach Etienne\n", "\n", "**This module is currently under development**\n", "\n", "## This Jupyter notebook provides an in-depth tutorial on computing the MHD flux terms within IllinoisGRMHD. It details the implementation of the Harten-Lax-van Leer (HLL) flux and elaborates on basic quantities, EOS quantities, speed limitations, comoving magnetic fields, characteristic speeds, and MHD flux for each component.\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](#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": [ "\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": [ "\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": [ "\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": [ "\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 contains 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 implements 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. its 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 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": [ "\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": [ "\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 the 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": [ "\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": [ "\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": [ "\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": [ "\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": [ "\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": [ "\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 (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 }