This module is currently under development
IllinoisGRMHD
solves the equations of General Relativistic MagnetoHydroDynamics (GRMHD) using a high-resolution shock capturing scheme. It is a rewrite of the Illinois Numerical Relativity (ILNR) group's GRMHD code, and generates results that agree to roundoff error with that original code. Its feature set coincides with the features of the ILNR group's recent code (ca. 2009-2014), which was used in their modeling of the following systems:
IllinoisGRMHD
is particularly good at modeling GRMHD flows into black holes without the need for excision. Its HARM-based conservative-to-primitive solver has also been modified to check the physicality of conservative variables prior to primitive inversion and moves them into the physical range if they become unphysical.
Currently, IllinoisGRMHD consists of
IllinoisGRMHD
evolves the vector potential $A_{\mu}$ (on staggered grids) instead of the magnetic fields ($B^i$) directly, to guarantee that the magnetic fields will remain divergenceless even at AMR boundaries. On uniform resolution grids, this vector potential formulation produces results equivalent to those generated using the standard, staggered flux-CT scheme. This scheme is based on that of Del Zanna et al. (2003).
The basic equations solved by IllinoisGRMHD
are
where $g_{\mu\nu}$ is the ADM 4-metric, $R_{\mu\nu}$ and $R$ are the Ricci tensor and scalar, respectively, constructed from $g_{\mu\nu}$, $\nabla_{\mu}$ is the covariant derivative compatible with $g_{\mu\nu}$, $g\equiv\det\left(g_{\mu\nu}\right)$, $\rho_0$ is the rest-mass density of the fluid, $u^{\mu}$ is the four-velocity of the fluid, $F^{\mu\nu}$ is the Faraday tensor and $F^{*\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\rho\sigma}F_{\mu\nu}$ its dual ($\epsilon^{\mu\nu\rho\sigma}$ is the Levi-Civita symbol).
The final equation that must be solved is the equation of state (EOS) for the matter. Currently IllinoisGRMHD
implements a hybrid EOS of the form
where P is the pressure, $\epsilon$ the specific internal energy (the subscript cold indicate the cold components of these quantities), and $\Gamma_{\rm th}$ is a constant parameter which determines the conversion efficiency of kinetic to thermal energy at shocks.
In this step, we will write down the equations used by IllinoisGRMHD
in the form they are implemented. To give an example, the GRMHD equations are written in the conservative form
where $\boldsymbol{C} = \left\{\rho_{\star},\tilde\tau,\tilde{S}_{i},\tilde{B}^{i}\right\}$ is the vector of conservative variables, $\boldsymbol{F}$ is the flux vector, and $\boldsymbol{S}$ the vector of source terms (explicit expressions for the components of $\boldsymbol{C}$, $\boldsymbol{F}$, and $\boldsymbol{S}$ can be found below).
IllinoisGRMHD
solves Einstein's field equations (with $G_{\rm N}=1=c$) in the presence of matter sources,
where $G^{\mu\nu} \equiv R^{\mu\nu} + \frac{1}{2}g^{\mu\nu}R$ is the Einstein tensor and the total energy-momentum tensor, $T^{\mu\nu}$, is the sum of the matter and electromagnetic energy-momentum tensors:
$$ T^{\mu\nu} = T^{\mu\nu}_{\rm matter} + T^{\mu\nu}_{\rm EM}\ . $$Einstein's field equations are solved using the BSSN formalism (see this tutorial module for an overview).
The conservation of the baryon number equation is written as (cf. eqs. (6), (7), and (18) in Etienne et al.)
$$ \boxed{\partial_{t}\rho_{\star} + \partial_{j}\left(\rho_{\star}v^{j}\right) = 0}\ , $$where $\rho_{\star}\equiv \alpha \sqrt{\gamma}\rho_{0}u^{0}$ and $v^{i}\equiv u^{i}/u^{0}$.
In the ideal MHD limit, we can write down the total energy-momentum tensor as (cf. eq. (8) in Etienne et al. and the discussion before it)
$$ T^{\mu \nu} = (\rho_0 h +b^2) u^{\mu} u^{\nu} + \left( P + \frac{b^2}{2}\right) g^{\mu \nu} - b^{\mu} b^{\nu} $$The spatial components of the energy-momentum conservation equation give (cf. eq. (18) of Etienne et al. and Eqs. (35) and (36) of Duez et al.)
$$ \boxed{\partial_{t}\tilde{S}_{i} + \partial_{j}\left(\alpha\sqrt{\gamma} T^{j}_{\ i}\right) = \frac{1}{2}\sqrt{\gamma} T^{\alpha\beta}\partial_{i}g_{\alpha\beta}}\ , $$where
$$ \tilde{S}_{i} = \sqrt{\gamma}S_{i} = \alpha \sqrt{\gamma} T^{0}_{\ i} = \left(\rho_{\star}h + \alpha u^{0}\sqrt{\gamma}b^{2}\right)u_{i} - \alpha \sqrt{\gamma} b^{0} b^{i}\ . $$The time component of the energy-momentum conservation equation gives (cf. eq. (18) of Etienne et al.)
$$ \boxed{\partial_{t}\tilde{\tau} + \partial_{i}\left(\alpha^{2}\sqrt{\gamma}T^{0i} - \rho_{\star}v^{i}\right) = s}\ , $$where
\begin{align} \tilde{\tau} &= \sqrt{\gamma}n_{\mu}n_{\nu} - \rho_{\star} = \alpha^{2}\sqrt{\gamma}T^{00} - \rho_{\star}\ ,\\ s &= -\alpha \sqrt{\gamma} T^{\mu\nu} \nabla_{\nu} n_{\mu} = \alpha\sqrt{\gamma}\left[\left(T^{00}\beta^{i}\beta^{j} + 2T^{0i}\beta^{j} + T^{ij}\right)K_{ij} - \left(T^{00}\beta^{i} + T^{0i}\right)\partial_{i}\alpha\right]\ , \end{align}with $n_{\mu} = \left(\alpha,0,0,0\right)$ being the normal vector and $K_{ij}$ the extrinsic curvature.
From the spatial components of the dual of Maxwell's equations,
$$ \nabla_{\nu} F^{*\mu\nu}=\frac{1}{\sqrt{-g}}\partial_{\nu}\left(\sqrt{-g}F^{*\mu\nu}\right)=0\ , $$we get the magnetic induction equation, which in conservative form may be written as (cf. eq. (12) Etienne et al.)
$$ \partial_{t}\tilde{B}^{i} + \partial_{j}\left(v^{j}\tilde{B}^{i} - v^{i}\tilde{B}^{j}\right) = 0\ , $$where $\tilde{B}^{i} = \sqrt{\gamma}B^{i}$. We must also guarantee that no magnetic monopoles form, via the constraint
$$ \partial_{i}\tilde{B}^{i} = 0\ , $$which is the time component of the dual of Maxwell's equations. Satisfying this constraint equation while evolving the magnetic field forward in time via the evolution equation above turns out to be a nontrivial endeavor, particularly on AMR grids. Instead, we choose to evolve the magnetic 4-vector potential $\mathcal{A}_{\mu}$ instead of the magnetic fields directly, so that
\begin{align} \mathcal{A}_{\mu} &= \Phi n_{\mu} + A_{\mu}\ ,\\ \tilde{B}^{i} &= \tilde{\epsilon}^{ijk}\partial_{j}A_{k}\ , \end{align}where $A_{\mu}$ is purely spatial ($A_{\mu}n^{\mu} = 0$) and $\Phi$ is the EM scalar potential.
Special finite difference operators for the vector potential are defined in IllinoisGRMHD so that the divergence of a curl is zero to roundoff error, which implies that the divergence of $\tilde{B}^{i}$ (as defined above) is zero and the condition $\partial_{i}\tilde{B}^{i}=0$ is satisfied automatically, even on AMR grids.
In terms of $A_{i}$, the induction equation becomes
$$ \boxed{\partial_t A_i = \tilde\epsilon_{ijk} v^j \tilde{B}^k - \partial_i (\alpha \Phi - \beta^j A_j)}\ . $$The final equation comes from choosing the covariant version of the "generalized Lorenz gauge condition" that was introduced by the Illinois Relativity group,
$$ \nabla_{\mu}\mathcal{A}^{\mu} = \xi n_{\mu} \mathcal{A}^{\mu}\ , $$where $\xi$ is a parameter with dimensions 1/Length, chosen so that the CFL condition remains satisfied. This gauge choice yields the additional evolution equation
$$ \boxed{\partial_t [\sqrt{\gamma} \Phi] + \partial_j (\alpha \sqrt{\gamma} A^j - \beta^j [\sqrt{\gamma}\Phi]) = -\xi\alpha \sqrt{\gamma}\Phi}\ . $$Note that in IllinoisGRMHD
the evolved variable is $\sqrt{\gamma}\Phi$, not $\Phi$.
IllinoisGRMHD
currently implements a hybrid EOS of the form
The function $\epsilon_{\rm cold}$ is related to $P_{\rm cold}$ by the first law of thermodynamics,
$$ \epsilon_{\rm cold}(\rho_{0}) = \int \frac{P_{\rm cold}(\rho_{0})}{\rho_{0}^{2}}d\rho_{0}\ . $$Currently, all functions within IllinoisGRMHD
support piecewise-defined $P_{\rm cold}(\rho_{0})$ (the so-called "piecewise polytrope" EOS) with up to nine different polytropic indices.
In typical runs of the code, particularly on the ones presented in the release paper, the $\Gamma$-law EOS, $P = (\Gamma - 1)\rho_{0}\epsilon$, is adopted. This corresponds to setting $P_{\rm cold} = (\Gamma - 1)\rho_{0}\epsilon_{\rm cold}$ in the boxed equation above, which is equivalent to $P_{\rm cold} = \kappa\rho_{0}^{\Gamma}$ (with constant $\kappa$), and $\Gamma_{\rm th} = \Gamma$. In the absence of shocks, $\epsilon = \epsilon_{\rm cold}$ so that $P = P_{\rm cold}$.
IllinoisGRMHD
has also performed successful binary neutron star mergers using the widely adopted piecewise polytropic EOSs of Read et al.. In the EOS low-level functions tutorial, you will find special Python
functions designed to help the user who wishes to use Read et al. EOSs or any other piecewise polytropic EOS.
#!jupyter nbconvert --to latex --template ../../latex_nrpy_style.tplx Tutorial-IllinoisGRMHD__Overview.ipynb
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Overview.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Overview.tex
#!pdflatex -interaction=batchmode Tutorial-IllinoisGRMHD__Overview.tex
!rm -f Tut*.out Tut*.aux Tut*.log