Adding an equation of state that includes radiation pressure

[shibersag]

Currently, Octo-Tiger implements two types of equation of states (EoSs):

1. Ideal - Only gas (matter) contributes to the pressure (i.e. thermal pressure). The gas is assumed to be an ideal gas with some adiabatic index \gamma. If this EoS is used with \epsilon=1, it can also function as a polytropic EoS.
2. WD - zero-temperature white dwarf EoS. In addition to thermal pressure it includes pressure from electrons in a degenerate state.

The first EoS is useful when simulating low mass main-sequence stars and the second is useful when simulating white dwarfs.

To model more massive stars in which the pressure from radiation is not negligible we need to incorporate into Octo-Tiger a new equation of state that includes radiation pressure.

This would bring Octo-Tiger to the realm of massive stars and would widen the possible scenarios we can simulate. Massive stars evolve to compact objects, like neutron stars or black holes, which most of the recently detected gravitational waves originated from, and as such there is currently a growing interest in binaries with massive stars.

The main difficulty in adding a radiation pressure term to the EoS is that the temperature and pressure cannot be easily derived from the internal energy but instead need to be solved numerically preferably by the Newton-Raphson method.

[shibersag]

I wrote a first implementation of this equation of state (IPR, ideal plus radiation). It is currently on eos_ideal_rad branch. @dmarce1 , @G-071 , would you mind taking a look? I'm currently running a simulation of massive binary with this implementation. I'll see how this simulation evolves.

[shibersag]

I would just say that with the current implementation there is a problem with the silo output. Since the pressure is no longer an analytic expression of the energy, the pressure and temperature that Octo-tiger outputs in the silo files is incorrect. @dmarce1 Should we add another field to store the correct pressure (or temperature)?

[G-071]

@shibersag Is there a good scenario we can use to test it (i.e the new EoS) with our Jenkins Pipeline? Once you and @dmarce1 judge the implementation to be correct, we should immediately set up the tests to check all further PRs to avoid it breaking when we start optimizing the performance/add more features.

When you open the PR and give me the scenario parameters I can set up such a test and the respective Jenkins Pipeline for continuous testing (note: the Jenkins machines we are using are quite powerful but it would be nice to have some smaller tests than also run on your average development machine).

[shibersag]

@G-071 That is a good idea. Once we will have the first prototype, we can set-up a test. I'm not sure though how small this test can be. I will have to see. BTW, do we have any test for eos=WD right now?

[G-071]

@shibersag Only a very basic one where I changed the rotating_star to this eos. It is sensitive enough to be triggered by the bug we fixed in #387, but I wouldn't mind getting a better/larger test scenario for that!

[shibersag]

@G-071 , @dmarce1, The ipr eos is quite ready. Please have a look on the eos_ideal_rad branch and have your comments. I ended up using the tau variable field to store the temperature in order for the output silo files to have all the necessary correct data. In case we will find a way to incorporate the dual-energy formalism within the ipr eos this will have to be changed.

I also added a floor density runtime parameter used for the SCF in this branch. This is unrelated directly to the ipr eos but it is something that I find useful to have, and which I already used in some of my massive star simulations. Massive stars span a large range in densities so the hard coded floor value wasn't low enough in compare to the surface density of the stars.

I will do some more final tests and then will send a pull request.

[diehlpk]

@shibersag Please do a pull request, so we can give you comments there.

[diehlpk]

refs #414

[diehlpk]

@shibersag can this be closed, if so please do?

[shibersag]

I need to do some more testing before closing it

[diehlpk]

@shibersag Can we close the issue?

[diehlpk]

@shibersag Can we close the issue?

[shibersag]

@diehlpk ok, let's closed it. But we should make a note that this feature is currently still experimental and hasn't been fully tested yet.