broken support for piecewise polytropes in IllinoisGRMHD

[Roland Haas]

Looking at the source code IllinoisGRMHD in the ET seems to support piecewise polytropic EOS, at the very least for a single piece by setting up IllinoisGRMHD’s gamma_th parameter and the first entry in the eos_params_arrays2 grid array (eg by setting ID_converter_ILGRMHD’s Initial_Gamma and Initial_K parameters). However D_converter_ILGRMHD contains an explicit check for a simple Gamma-law (ie no difference between thermal and cold Gamma) in

// P = (\Gamma - 1) rho epsilon
// -> \Gamma = P/(rho epsilon) + 1
double measured_gamma = ( P[index]/(rho_b[index] * eps[index]) + 1.0 );
if(rho_b[index]>rho_b_atm && fabs(gamma_th - measured_gamma)/gamma_th > 1e-2)
  CCTK_VError(VERR_DEF_PARAMS,
              "Expected simple gamma law with gamma_th=%.15e, but found a point with gamma law such that gamma_th=%.15e. error = %e| rb=%e rbatm=%e P=%e\n",
              gamma_th,measured_gamma,(gamma_th-measured_gamma)/gamma_th,rho_b[index],rho_b_atm,P[index] );

‌

There is also no documentation (at all as far as I can tell) on what EOS IllinoisGRMHD supports or how one would go about setting them up.

[Zach Etienne]

The main documentation (the IllinoisGRMHD paper https://arxiv.org/pdf/1501.07276.pdf ) says:

“… the current version supports only single gamma-law EOSs…”

Full support for the “hybrid equation of state” (of which piecewise polytrope is a limiting case) was incorporated into the NRPy+ version of IllinoisGRMHD, which is rather pedagogically documented in many Jupyter notebooks:
https://nbviewer.jupyter.org/github/zachetienne/nrpytutorial/tree/master/IllinoisGRMHD/doc/

We proposed this version for inclusion in the Toolkit but it was rejected due to very minor backwards incompatibility.

[Roland Haas]

The current code has all kinds of things that seems to have support for piecwiese polytropes (at the very least for a single piece) eg the function compute_P_cold__eps_cold__dPcold_drho__eps_th__h__gamma_cold. Are all those incorrectly coded and only work for a Gamma-law?

I had a look at the notebooks, none tell me how to use the code to set up piecewise polytropes. They do go into details on how the piecewise polytropes may work showing first python code and then C code that (presumably) does the same thing. I could not find any information how to use it though. Admittedly I did not read all the 10+ jupyter notebooks so it could be in one the ones I did not read.

‌

[Roland Haas]

Note that the code was rejected for the 2020_05 release since it requires retiring a parameter, which (per policy provided in the ticket #2300) requires deprecation of the parameter one release in advance. After that it can be removed. At this point the deprecation notice would have to be included in 2021_05 and the parameter can be remove in 2021_11 (or the deprecation could have been included in 2020_11 and removed in 2021_05, but no action was taken).

[Roland Haas]

Turns out there is indeed code in IllinoisGRMHD (in the con2prim routines) that make explicit assumptions about the value of IllinoisGRMHD::gamma_th and ID_convert_ILGRMHD::Gamma_Initial (which turns into the gamma_tab array in IllinoisGRMHD).