{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "\n", "# [Polytropic TOV](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation) Initial Data\n", "\n", "## Authors: Phil Chang, Zach Etienne, & Leo Werneck\n", "### Formatting improvements courtesy Brandon Clark\n", "\n", "## This module sets up initial data for a [TOV](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation) star in *spherical, isotropic coordinates*\n", "\n", "**Notebook Status:** Validated \n", "\n", "**Validation Notes:** This module has been validated to exhibit convergence to zero of the Hamiltonian constraint violation at the expected order to the exact solution (see [start-to-finish TOV module](Tutorial-Start_to_Finish-BSSNCurvilinear-Setting_up_TOV_initial_data.ipynb) for full test). Note that convergence at the surface of the star is lower order due to the sharp drop to zero in $T^{\\mu\\nu}$.\n", "\n", "### NRPy+ Source Code for this module: [TOV/TOV_Solver.py](../edit/TOV/TOV_Solver.py)\n", "\n", "[comment]: <> (Introduction: TODO)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Table of Contents\n", "$$\\label{toc}$$\n", "\n", "This notebook is organized as follows:\n", "\n", "1. [Step 1](#initializenrpy): Initialize core Python/NRPy+ modules\n", "1. [Step 2](#tov): The TOV Equations\n", "1. [Step 3](#code_validation): Code Validation against `TOV.TOV_Solver` NRPy+ module\n", "1. [Step 4](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 1: Initialize core Python/NRPy+ modules \\[Back to [top](#toc)\\]\n", "$$\\label{initializenrpy}$$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:06:26.375685Z", "iopub.status.busy": "2021-03-07T17:06:26.374015Z", "iopub.status.idle": "2021-03-07T17:06:26.609669Z", "shell.execute_reply": "2021-03-07T17:06:26.610190Z" } }, "outputs": [], "source": [ "# Step 1: Import needed Python/NRPy+ modules\n", "import numpy as np # NumPy: A numerical methods module for Python\n", "import scipy.integrate as si # SciPy: Python module for mathematics, science, and engineering applications\n", "import math, sys # Standard Python modules for math; multiplatform OS-level functions\n", "import TOV.Polytropic_EOSs as ppeos # NRPy+: Piecewise polytrope equation of state support" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 2: The TOV equations \\[Back to [top](#toc)\\]\n", "$$\\label{tov}$$\n", "\n", "The [TOV line element](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation) in terms of the *Schwarzschild coordinate* $r$ is written (in the $-+++$ form):\n", "$$\n", "ds^2 = - c^2 e^\\nu dt^2 + \\left(1 - \\frac{2Gm}{rc^2}\\right)^{-1} dr^2 + r^2 d\\Omega^2,\n", "$$\n", "where $m(r)$ is the mass-energy enclosed at a given $r$, and is equal to the total star's mass outside the stellar radius $r=R$.\n", "\n", "In terms of the *isotropic coordinate* $\\bar{r}$ with $G=c=1$ (i.e., the coordinate system and units we'd prefer to use), the ($-+++$ form) line element is written:\n", "$$\n", "ds^2 = - e^{\\nu} dt^2 + e^{4\\phi} \\left(d\\bar{r}^2 + \\bar{r}^2 d\\Omega^2\\right),\n", "$$\n", "where $\\phi$ here is the *conformal factor*.\n", "\n", "Setting components of the above line element equal to one another, we get (in $G=c=1$ units):\n", "\n", "\\begin{align}\n", "r^2 &= e^{4\\phi} \\bar{r}^2 \\implies e^{4\\phi} = \\frac{r^2}{\\bar{r}^2} \\\\\n", "\\left(1 - \\frac{2m}{r}\\right)^{-1} dr^2 &= e^{4\\phi} d\\bar{r}^2 \\\\\n", "\\implies \\frac{d\\bar{r}(r)}{dr} &= \\left(1 - \\frac{2m}{r} \\right)^{-1/2} \\frac{\\bar{r}(r)}{r}.\n", "\\end{align}\n", "\n", "The TOV equations provide radial ODEs for the pressure and $\\nu$ (from [the Wikipedia article on the TOV solution](https://en.wikipedia.org/wiki/Tolman%E2%80%93Oppenheimer%E2%80%93Volkoff_equation)):\n", "\n", "\\begin{align}\n", "\\frac{dP}{dr} &= - \\frac{1}{r} \\left( \\frac{\\rho + P}{2} \\right) \\left(\\frac{2 m}{r} + 8 \\pi r^2 P\\right) \\left(1 - \\frac{2 m}{r}\\right)^{-1} \\\\\n", "\\frac{d \\nu}{d r} &= \\frac{1}{r}\\left(1 - \\frac{2 m}{r}\\right)^{-1} \\left(\\frac{2 m}{r} + 8 \\pi r^2 P\\right) \\\\\n", "\\end{align}\n", "\n", "Assuming a polytropic equation of state, which relates the pressure $P$ to the baryonic rest-mass density $\\rho_B$,\n", "\n", "$$\n", "P(\\rho_B) = K \\rho_B^\\Gamma,\n", "$$\n", "the specific internal energy will be given by\n", "$$\n", "\\epsilon = \\frac{P}{\\rho_B (\\Gamma - 1)},\n", "$$\n", "\n", "so the total mass-energy density $\\rho$ is given by\n", "$$\n", "\\rho = \\rho_B (1 + \\epsilon).\n", "$$\n", "\n", "Given this, the mass-energy $m(r)$ density is the solution to the ODE:\n", "$$\n", "\\frac{dm(r)}{dr} = 4\\pi r^2 \\rho(r)\n", "$$\n", "\n", "Thus the full set of ODEs that need to be solved is given by\n", "\n", "$$\n", "\\boxed{\n", "\\begin{matrix}\n", "\\frac{dP}{dr} &=& - \\frac{1}{r} \\left( \\frac{\\rho + P}{2} \\right) \\left(\\frac{2 m}{r} + 8 \\pi r^2 P\\right) \\left(1 - \\frac{2 m}{r}\\right)^{-1} \\\\\n", "\\frac{d \\nu}{d r} &=& \\frac{1}{r}\\left(1 - \\frac{2 m}{r}\\right)^{-1} \\left(\\frac{2 m}{r} + 8 \\pi r^2 P\\right) \\\\\n", "\\frac{m(r)}{dr} &=& 4\\pi r^2 \\rho(r) \\\\\n", "\\frac{d\\bar{r}(r)}{dr} &=& \\left(1 - \\frac{2m}{r} \\right)^{-1/2} \\frac{\\bar{r}(r)}{r}\n", "\\end{matrix}\n", "}\\ .\n", "$$\n", "\n", "The following code solves these equations, and was largely written by Phil Chang." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:06:26.774454Z", "iopub.status.busy": "2021-03-07T17:06:26.639256Z", "iopub.status.idle": "2021-03-07T17:06:26.791454Z", "shell.execute_reply": "2021-03-07T17:06:26.791961Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1256 1256 1256 1256 1256 1256\n", "Just generated a TOV star with R_Schw = 9.566044579232513e-01 , M = 1.405030336771405e-01 , M/R_Schw = 1.468768334847266e-01 .\n" ] } ], "source": [ "# Step 2: The TOV equations\n", "\n", "## TOV SOLVER FOR SINGLE AND PIECEWISE POLYTROPES\n", "## Authors: Phil Chang, Zachariah B. Etienne, Leo Werneck\n", "\n", "# Full documentation for this module may be found in the NRPy+ tutorial Jupyter notebook:\n", "# Tutorial-Start_to_Finish-BSSNCurvilinear-Setting_up_TOV_initial_data.ipynb\n", "\n", "\n", "# Inputs:\n", "# * Output data file name\n", "# * rho_baryon_central, the central density of the TOV star.\n", "# * n, the polytropic equation of state index. n=1 models cold, degenerate neutron star matter.\n", "# * K_Polytrope, the polytropic constant.\n", "# * Verbose output toggle (default = True)\n", "\n", "# Output: An initial data file (default file name = \"outputTOVpolytrope.txt\") that well\n", "# samples the (spherically symmetric) solution both inside and outside the star.\n", "# It is up to the initial data module to perform the 1D interpolation to generate\n", "# the solution at arbitrary radius. The file has the following columns:\n", "# Column 1: Schwarzschild radius\n", "# Column 2: rho(r), *total* mass-energy density (as opposed to baryonic rest-mass density)\n", "# Column 3: P(r), Pressure\n", "# Column 4: m(r), mass enclosed\n", "# Column 5: e^{nu(r)}, g_{tt}(r)\n", "# Column 6: e^{4 phi(r)}, conformal factor g_{rr}(r)\n", "# Column 7: rbar(r), Isotropic radius\n", "\n", "# rbar refers to the isotropic radius, and\n", "# R_Schw refers to the Schwarzschild radius\n", "\n", "def TOV_Solver(eos,\n", " outfile = \"outputTOVpolytrope.txt\",\n", " rho_baryon_central = 0.129285,\n", " verbose = True,\n", " return_M_and_RSchw = False,\n", " accuracy = \"medium\",\n", " integrator_type = \"default\",\n", " no_output_File = False,\n", " pressure_renormalization=1): # reset the pressure to stellar oscillations studies\n", "\n", " def TOV_rhs(r_Schw, y) :\n", " # In \\tilde units\n", " #\n", " P = y[0]\n", " m = y[1]\n", " # nu = y[2] # nu is not needed as input into TOV_rhs\n", " rbar = y[3]\n", "\n", " # Compute rho_b and eps_cold, to be used below\n", " # to compute rho_(total)\n", " rho_baryon, eps_cold = ppeos.Polytrope_EOS__compute_rhob_and_eps_cold_from_P_cold(eos,P)\n", "\n", "# with open(\"rhob_P_cold_and_eps_cold.dat\",\"a+\") as file:\n", "# file.write(str(r_Schw).format(\"%.15e\")+\" \"+str(rho_baryon).format(\"%.15e\")+\" \"+str(P).format(\"%.15e\")+\" \"+str(eps_cold).format(\"%.15e\")+\"\\n\")\n", "\n", " # Compute rho, the *total* mass-energy density:\n", " # .------------------------------.\n", " # | rho = (1 + eps)*rho_(baryon) |\n", " # .------------------------------.\n", " # with eps = eps_cold, for the initial data.\n", " rho = (1.0 + eps_cold)*rho_baryon\n", "\n", "# m = 4*math.pi/3. * rho*r_Schw**3\n", " if( r_Schw < 1e-4 or m <= 0.):\n", " # From https://github.com/natj/tov/blob/master/tov.py#L33:\n", " # dPdr = -cgs.G*(eden + P/cgs.c**2)*(m + 4.0*pi*r**3*P/cgs.c**2)\n", " # dPdr = dPdr/(r*(r - 2.0*cgs.G*m/cgs.c**2))\n", " dPdrSchw = -(rho + P)*(4.*math.pi/3.*r_Schw*rho + 4.*math.pi*r_Schw*P)/(1.-8.*math.pi*rho*r_Schw*r_Schw)\n", " drbardrSchw = 1./(1. - 8.*math.pi*rho*r_Schw*r_Schw)**0.5\n", " else:\n", " dPdrSchw = -(rho + P)*(m + 4.*math.pi*r_Schw**3*P)/(r_Schw*r_Schw*(1.-2.*m/r_Schw))\n", " drbardrSchw = 1./(1. - 2.*m/r_Schw)**0.5*rbar/r_Schw\n", "\n", " dmdrSchw = 4.*math.pi*r_Schw*r_Schw*rho\n", " dnudrSchw = -2./(P + rho)*dPdrSchw\n", " return [dPdrSchw, dmdrSchw, dnudrSchw, drbardrSchw]\n", "\n", " def integrateStar( eos, P, dumpData = False ):\n", "\n", " if accuracy == \"medium\":\n", " min_step_size = 1e-5\n", " max_step_size = 1e-2\n", " integrator = 'dop853'\n", " elif accuracy == \"low\":\n", " min_step_size = 1e-3\n", " max_step_size = 1e-1\n", " integrator = 'dopri5'\n", " elif accuracy == \"verylow\":\n", " min_step_size = 1e-1\n", " max_step_size = 5e-1\n", " integrator = 'dopri5'\n", " elif accuracy == \"high\":\n", " min_step_size = 1e-5\n", " max_step_size = 1e-5\n", " integrator = 'dop853'\n", " elif accuracy == \"veryhigh\":\n", " min_step_size = 1e-7\n", " max_step_size = 1e-6\n", " integrator = 'dop853'\n", " else:\n", " print(\"Unknown accuracy option: \"+str(accuracy))\n", "\n", " if integrator_type == \"default\":\n", " pass\n", " else:\n", " integrator = integrator_type\n", "\n", " integrator = si.ode(TOV_rhs).set_integrator(integrator)#,rtol=1e-4,atol=1e-4)\n", "# integrator = si.ode(TOV_rhs).set_integrator('dopri5',rtol=1e-4)\n", " y0 = [P, 0., 0., 0.]\n", " r_Schw = 0.\n", " integrator.set_initial_value(y0,r_Schw)\n", " dr_Schw = min_step_size\n", " P = y0[0]\n", "\n", " PArr = []\n", " r_SchwArr = []\n", " mArr = []\n", " nuArr = []\n", " rbarArr = []\n", "\n", " while integrator.successful() and P > 1e-19*y0[0] :\n", " P, m, nu, rbar = integrator.integrate(r_Schw + dr_Schw)\n", " # Update the value of r_Schw to the latest integrated value\n", " r_Schw += dr_Schw\n", "\n", " dPdrSchw, dmdrSchw, dnudrSchw, drbardrSchw = TOV_rhs( r_Schw, [P,m,nu,rbar])\n", " dr_Schw = 0.1*min(abs(P/dPdrSchw), abs(m/dmdrSchw))\n", " dr_Schw = min(dr_Schw, max_step_size)\n", " PArr.append(P)\n", " r_SchwArr.append(r_Schw)\n", " mArr.append(m)\n", " nuArr.append(nu)\n", " rbarArr.append(rbar)\n", "\n", " M = mArr[-1]\n", " R_Schw = r_SchwArr[-1]\n", "\n", " if no_output_File == True:\n", " return R_Schw, M\n", "\n", " # Apply integration constant to ensure rbar is continuous across TOV surface\n", " for ii in range(len(rbarArr)):\n", " rbarArr[ii] *= 0.5*(np.sqrt(R_Schw*(R_Schw - 2.0*M)) + R_Schw - M) / rbarArr[-1]\n", "\n", " nuArr_np = np.array(nuArr)\n", " # Rescale solution to nu so that it satisfies BC: exp(nu(R))=exp(nutilde-nu(r=R)) * (1 - 2m(R)/R)\n", " # Thus, nu(R) = (nutilde - nu(r=R)) + log(1 - 2*m(R)/R)\n", " nuArr_np = nuArr_np - nuArr_np[-1] + math.log(1.-2.*mArr[-1]/r_SchwArr[-1])\n", "\n", " r_SchwArrExtend_np = 10.**(np.arange(0.01,5.0,0.01))*r_SchwArr[-1]\n", "\n", " r_SchwArr.extend(r_SchwArrExtend_np)\n", " mArr.extend(r_SchwArrExtend_np*0. + M)\n", " PArr.extend(r_SchwArrExtend_np*0.)\n", " exp2phiArr_np = np.append( np.exp(nuArr_np), 1. - 2.*M/r_SchwArrExtend_np)\n", " nuArr.extend(np.log(1. - 2.*M/r_SchwArrExtend_np))\n", " rbarArr.extend( 0.5*(np.sqrt(r_SchwArrExtend_np**2 - 2.*M*r_SchwArrExtend_np) + r_SchwArrExtend_np - M) )\n", "\n", " #phiArr_np = np.append( np.exp(nuArr_np), 1. - 2.*M/r_SchwArrExtend_np)\n", "\n", " # Appending to a Python array does what one would reasonably expect.\n", " # Appending to a numpy array allocates space for a new array with size+1,\n", " # then copies the data over... over and over... super inefficient.\n", " r_SchwArr_np = np.array(r_SchwArr)\n", " PArr_np = np.array(PArr)\n", " rho_baryonArr_np = np.array(PArr) # This is just to initialize the array\n", " for j in range(len(PArr_np)):\n", " # Compute rho_b from P\n", " rho_baryonArr_np[j] = ppeos.Polytrope_EOS__compute_rhob_from_P_cold(eos,PArr_np[j])\n", "\n", " mArr_np = np.array(mArr)\n", " rbarArr_np = np.array(rbarArr)\n", " confFactor_exp4phi_np = (r_SchwArr_np/rbarArr_np)**2\n", "\n", " # Compute the *total* mass-energy density (as opposed to the *baryonic* mass density)\n", " rhoArr_np = []\n", " for i in range(len(PArr)):\n", " rho_baryon, eps_cold = ppeos.Polytrope_EOS__compute_rhob_and_eps_cold_from_P_cold(eos,PArr[i])\n", " rho = (1.0 + eps_cold ) * rho_baryon\n", " rhoArr_np.append(rho)\n", "\n", " if verbose:\n", " print(len(r_SchwArr_np),len(rhoArr_np),len(rho_baryonArr_np),len(PArr_np),len(mArr_np),len(exp2phiArr_np))\n", "\n", " PArr_np *= pressure_renormalization # set for pressure renormalization studies\n", "\n", " # Special thanks to Leonardo Werneck for pointing out this issue with zip()\n", " if sys.version_info[0] < 3:\n", " np.savetxt(outfile, zip(r_SchwArr_np,rhoArr_np,rho_baryonArr_np,PArr_np,mArr_np,exp2phiArr_np,confFactor_exp4phi_np,rbarArr_np),\n", " fmt=\"%.15e\")\n", " else:\n", " np.savetxt(outfile, list(zip(r_SchwArr_np,rhoArr_np,rho_baryonArr_np,PArr_np,mArr_np,exp2phiArr_np,confFactor_exp4phi_np,rbarArr_np)),\n", " fmt=\"%.15e\")\n", "\n", " return R_Schw, M\n", "\n", " # Set initial condition from rho_baryon_central\n", " P_initial_condition = ppeos.Polytrope_EOS__compute_P_cold_from_rhob(eos, rho_baryon_central)\n", "\n", " # Integrate the initial condition\n", " R_Schw_TOV, M_TOV = integrateStar(eos, P_initial_condition, True)\n", " if verbose:\n", " print(\"Just generated a TOV star with R_Schw = %.15e , M = %.15e , M/R_Schw = %.15e .\" %(R_Schw_TOV,M_TOV,(M_TOV / R_Schw_TOV)))\n", "\n", " if return_M_and_RSchw:\n", " return M_TOV, R_Schw_TOV\n", "\n", "############################\n", "# Single polytrope example #\n", "############################\n", "# Set neos = 1 (single polytrope)\n", "neos = 1\n", "\n", "# Set rho_poly_tab (not needed for a single polytrope)\n", "rho_poly_tab = []\n", "\n", "# Set Gamma_poly_tab\n", "Gamma_poly_tab = [2.0]\n", "\n", "# Set K_poly_tab0\n", "K_poly_tab0 = 1. # ZACH NOTES: CHANGED FROM 100.\n", "\n", "# Set the eos quantities\n", "eos = ppeos.set_up_EOS_parameters__complete_set_of_input_variables(neos,rho_poly_tab,Gamma_poly_tab,K_poly_tab0)\n", "\n", "# Set initial condition (Pressure computed from central density)\n", "rho_baryon_central = 0.129285\n", "\n", "M_TOV, R_Schw_TOV = TOV_Solver(eos,outfile=\"outputTOVpolytrope.txt\",rho_baryon_central=0.129285,verbose = True,\n", " return_M_and_RSchw=True,accuracy=\"medium\",integrator_type=\"default\",\n", " pressure_renormalization=1.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 3: Code Validation \\[Back to [top](#toc)\\]\n", "$$\\label{code_validation}$$\n", "\n", "Here, as a code validation check, we verify agreement in the SymPy expressions for these TOV initial data between\n", "\n", "1. this tutorial and \n", "2. the NRPy+ [TOV.TOV_Solver](../edit/TOV/TOV_Solver.py) module." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:06:26.937141Z", "iopub.status.busy": "2021-03-07T17:06:26.797807Z", "iopub.status.idle": "2021-03-07T17:06:26.954030Z", "shell.execute_reply": "2021-03-07T17:06:26.954518Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1256 1256 1256 1256 1256 1256\n", "Just generated a TOV star with\n", "* M = 1.405030336771405e-01 ,\n", "* R_Schw = 9.566044579232513e-01 ,\n", "* R_iso = 8.100085557410308e-01 ,\n", "* M/R_Schw = 1.468768334847266e-01 \n", "\n", "TOV initial data test PASSED.\n" ] } ], "source": [ "# Step 3: Code Validation against TOV.TOV_Solver module\n", "\n", "import filecmp\n", "\n", "import TOV.TOV_Solver as TOV\n", "TOV.TOV_Solver(eos,\n", " outfile=\"outputTOVpolytrope-validation.txt\",\n", " rho_baryon_central=0.129285,\n", " verbose = True,\n", " accuracy=\"medium\",\n", " integrator_type=\"default\",\n", " no_output_File = False)\n", "\n", "if filecmp.cmp('outputTOVpolytrope.txt',\n", " 'outputTOVpolytrope-validation.txt') == False:\n", " print(\"ERROR: TOV initial data test FAILED!\")\n", " sys.exit(1)\n", "else:\n", " print(\"TOV initial data test PASSED.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "# Step 4: 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-ADM_Initial_Data-TOV](Tutorial-ADM_Initial_Data-TOV.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": 4, "metadata": { "execution": { "iopub.execute_input": "2021-03-07T17:06:26.960083Z", "iopub.status.busy": "2021-03-07T17:06:26.958715Z", "iopub.status.idle": "2021-03-07T17:06:30.249332Z", "shell.execute_reply": "2021-03-07T17:06:30.250235Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Created Tutorial-ADM_Initial_Data-TOV.tex, and compiled LaTeX file to PDF\n", " file Tutorial-ADM_Initial_Data-TOV.pdf\n" ] } ], "source": [ "import cmdline_helper as cmd # NRPy+: Multi-platform Python command-line interface\n", "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-ADM_Initial_Data-TOV\")" ] } ], "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.2" } }, "nbformat": 4, "nbformat_minor": 2 }