EOS_Omni::poly_gamma_ini should default to poly_gamma

[Frank Löffler]

At the moment, EOS_Omni::poly_gamma and EOS_Omni::poly_gamma_ini are two parameters for the adiabatic exponent, one used only for initial data, the other for evolutions. Both default to 2.0. This means that if someone, unaware of poly_gamma_ini, changes only poly_gamma the result is usually catastrophic (wrong results or c2p failures). I myself didn't notice today and had to debug my way down to find my mistake.

I propose to make the default of EOS_Omni::poly_gamma_ini to be whatever EOS_Omni::poly_gamma is. One way to achieve that would be to use some special value for EOS_Omni::poly_gamma_ini which would usually not being used to indicate that behavior and to make that default. However, this being an exponent, at least in theory, all real values are allowed. Another approach would be to introduce a new boolean parameter which would allow users to use an initial value for gamma which is different from the one of the initial data. This would be 'saver', but has the drawback of a new parameter. Nevertheless, I favor the second approach.

Keyword:

[Frank Löffler]

I put this up for review to decide about which option to use, will then prepare a patch, and put this up for review again then.

[Bruno Mundim]

the second approach with the boolean parameter looks better...

[Erik Schnetter]

If there is a boolean, then there also needs to be a check that people don't modify the *_ini parameter without modifying the boolean.

It would be simpler to set the default of the *_ini parameter to -1000000 (this will never be used in an actual simulation). In this way, if one sets the *_ini parameter, it will be used, and existing parameter files do not need to be changed.

[Bruno Mundim]

changing the default to -100000000 seems even better!

[Roland Haas]

I am a little bit confused as to what poly_gamma_ini does. Doing a grep for this name on the EOS_Omni sources it seems to be used only in EOS_Omni_Startup.F90 to set up poly_k_cgs (and similar quantities for gl and hybrid eos). poly_k_cgs seems then to be used unmodified during the lifetime of the run. So poly_gamma_ini does not seem to just control values during initial data (which I would find somewhat unexpected for an EOS [which should be just a table of states] to do) but also during runtime. Is this correct?

It seems as if poly_gamma_ini is used to construct a conversion factor from M_sun=c=G=1 units to cgs units (which conversion we seem to do backwards and forwards for the analytic EOS calls).

[Frank Löffler]

Right. I also find it a bit confusing to have cgs units used for constants in the code but then transforming back and forth all the time to cancel them out again. At the moment, a polytrope is calculated as p = K * rho^poly_gamma * magic_number^(poly_gamma_ini - poly_gamma) where magic_number is defined in EOS_Omni_Module.F90 with the comment:

Magic numbers in convenient units, calculated by someone whom we trust to have worked very carefully

I think there are two separate issues here. The first is the question why these units are used in the first place - why converting poly_k to poly_k_cgs if inside the EOS we use p_gf and rho_gf to reverse it again.

The second issue is that at the moment poly_gamma_ini isn't really doing what it is supposed to be doing. It would be correct if poly_k_cgs would be reset to poly_k * rho_gf^poly_gamma / p_gf after it has been used for initial data.

[Erik Schnetter]

Okay; this comment should definitively have been caught by code review. The comment is funny to read (and could stay), but a real explanation needs to follow, including listing the actual constants used, their units, and their source (e.g. Wikipedia).

[Roland Haas]

Going one step further: should we remove poly_gamma_ini altogether as well as (for the Gamma-law, polytropic and possibly hybrid EOS) remove the conversion from/to cgs units that happens in each EOS call? They seem to not do anything since it always converts back and forth and no user visible quantity uses cgs units.

We would still need comments explaining the magic numbers (or replace them by values from thorn Constants if they are the same).

[Frank Löffler]

I don't use poly_gamma_ini and never did. I wrote an email to the users mailing list. I agree with the constants.

[Roland Haas]

It's been a week since Frank send out the email. Not objections were voiced. I'd say lets remove poly_gamma_ini.

The email did not mention removing the conversion so there we should check if there is any possible way that this conversion is visible to the user and if not, I would prefer to remove it as well.

[Frank Löffler]

Removing 'review' until I have a patch.

[Roland Haas]

added patches for EOS_Omni, GRHydro (which used the module) and TOVSolver (for a test case for the hybrid EOS).

[Roland Haas]

Replying to [comment:7 eschnett]:

Okay; this comment should definitively have been caught by code review. The comment is funny to read (and could stay), but a real explanation needs to follow, including listing the actual constants used, their units, and their source (e.g. Wikipedia).

I just did that and attach a diff of the module file. The results change a bit (order 1-e2) which is a bit more than what I would have expected from using a slightly different solar mass.

I took the constant (and the solar mass) from http://asa.usno.navy.mil/SecK/2013/Astronomical_Constants_2013.pdf which says

The IAU 2009 System of Astronomical Constants (1) as published in the Report of the IAU Working Group on Numerical Standards for Fundamental Astronomy (NSFA, 2011), (2) planetary equatorial radii, taken from the report of the IAU WG on Cartographic Coordinates and Rotational Elements: 2009 (2011), and lastly (3) other useful constants. and which is in fact the primary reference on wikipedia (in the solar mass article).

For fun I also looked at thorn Constants and found that (a) it does not give a reference for its numbers and (b) its gravitational constant G=6.6732e-11 differs from the one I used 6.67428d−11 (all SI units). G is just hard to measure it seems. Using Google to look for the number used, it might have come from RNSID (http://www.gravity.phys.uwm.edu/rns/source/rns.v2.0/consts.h) which also does not seem to give a source.

[Frank Löffler]

Constants used values from TOVSolver (which now uses Constants), without changing them, and without source there as well. I suggest to update Constants to use the currently 'accepted' values, and of course to state the source. EOS_Omni could then use Constants as well. The purpose of Constants is not only to have these constants handy, but to have only one value for a specific constant used within a code. If that would not be the case cancelations might actually give surprising results.

[Frank Löffler]

Changing Constants should only affect the screen output within TOVSolver. They are not used for something else.

[Roland Haas]

The numbers is EOS_Omni are used when one uses the tabulated EOS (for the analytic EOS the values should all cancel eventually).

When (and I agree we should) change them, we will need to warn users that their results might change. There are also likely other places in the toolkit where M_sun and/or G,c are used (or carefully hidden in some typed in number).

[Frank Löffler]

There is now a ticket about removing the parameter altogether (#1245).