Ideal EOS Returns Negative Heat Capacity at high T / low Rho

[evbauer]

See top left corner of attached plot from EOS plotter for just the ideal gas+radiation EOS at solar composition, where red represents a bad value (either NaN or negative value returned by ideal EOS).

Ran into this while working with Shelley Cheng on some models with high mass loss rates. Looks likely to be something like a floating point subtraction issue, so maybe @fxt44 can offer some advice on how to patch this sort of thing up.

The red in the lower-right corner of the plot is a separate issue described in issue #565.

[evbauer]

Edit: had a dyslexic moment when I first wrote the title. It's high temperature and low density that's relevant for this issue. The opposite corner is the topic of #565.

[evbauer]

The underlying issue here is that $c_P = c_V \Gamma1 / \chi\rho$ and $\chi\rho \rightarrow 0$ for high temperature and low density because the radiation pressure dominates over the gas pressure. Where we're having problems is the $\chi\rho \lesssim 10^{-16}$ region. Floating point errors occasionally return a negative value there.

Maybe the simplest hack to avoid this is to set a floor of 1d-99 somewhere around here: https://github.com/MESAHub/mesa/blob/acd763a682009d91e6ff6f1023b9cf72d21bf8dc/eos/private/skye_thermodynamics.f90#L65

Would that be good enough or should we consider trying something more sophisticated to avoid the floating point issues in the first place? I'll want to check the other EOS quantities in this region as well to see if anything else is returning bad values.

My main motivation here is just to help models avoid crashing in early Newton iterations that might touch this regime. Obviously stellar profiles shouldn't actually live here, but I always want to squash bad numerics in the input physics wherever I can find them so that we can give MESA a better chance at iterating toward a correct solution.

[fxt44]

for t=1e10 k, rho=1e-12 g/cc, pure hydrogen, i get cp = 1.25e-20 erg/k from a helm-style interpolation and a direct calculation. subtraction of two large and nearly equal numbers in a helm-style eos is always a potential danger in a radiation dominated regime (which never occurs in a direct calculation).

[evbauer]

This reveals an interesting interaction between floating point arithmetic and the fact that Skye/ideal are built on top of auto_diff and the free energy. I was wondering how we could end up with subtraction issues with a simple analytic EOS like ideal, and it turns out it's actually pretty easy to see by writing out everything explicitly for just the radiation term and its implementation.

Of course, radiation pressure has no density dependence, so we should have $\chi_\rho = \partial \ln P /\partial \ln \rho= 0$ everywhere if we use a pure radiation EOS. Instead, we get the attached plot of floating point noise because somehow the EOS implementation ends up trying to subtract two numbers to recover zero. Why? This is a consequence of deriving everything from the free energy using autodiff. The free energy of the radiation term is $$F{\rm rad} = -\frac{aT^4}{3\rho}$$ The pressure is derived from this as $$P{\rm rad} = \rho^2\frac{\partial F{\rm rad}}{\partial \rho} = \rho^2 \frac{aT^4}{3\rho^2}$$ Analytically, of course, that reduces to $P_{\rm rad} = aT^4/3$ like we want, but at least in our current implementation auto_diff sees the above expression when taking derivatives to construct derived quantities for the EOS. Since autodiff doesn't know how to do the cancellation of the $\rho^2$ term in advance, it will end up implementing something like $\chi\rho$ by blind application of the chain rule and product rule: $$\chi\rho \equiv \frac{\rho}{P}\frac{\partial P}{\partial \rho} = \frac{\rho}{P} \left(2 \rho \frac{aT^4}{3\rho^2} - 2\rho^2 \frac{aT^4}{3\rho^3}\right)$$ That's what gives us a subtraction leading to all the noise and bad values in $\chi\rho$. In the full ideal EOS, once the radiation pressure dominates over the ion pressure term, this noise surfaces and gives bad values for things like $\chi_\rho$ and $c_P$.

I don't have an immediate suggestion for the best path toward fixing this at the moment, but I think this diagnosis should at least inform a start on the solution. I'll spend some time considering what kind of solution this motivates.

[fxt44]

very good! now you see why subtractions when starting from a helmholtz freee energy can be an issue.

how about doing radiation analytically and adding it in at the end? both helm (standalone) and more direct approaches do this ...

[evbauer]

how about doing radiation analytically and adding it in at the end? both helm (standalone) and more direct approaches do this ...

Yes, was just thinking something along these lines... Initial PR incoming. This makes me also want to think through if there are any other obvious places this kind of issue might pop up in Skye and ideal.

[fxt44]

should you want to save some work (or maybe it just adds some work), here are all the fundamental derivatatives:

% more eosfxt_radiation.f90

! pressure in erg/cm**3
        prad    = asol * third * temp * temp * temp * temp

! first derivatives
        dpraddd = 0.0d0
        dpraddt = 4.0d0 * prad * tempinv
        dpradda = 0.0d0
        dpraddz = 0.0d0

! second derivatives
        dpradddd = 0.0d0
        dpradddt = 0.0d0
        dpraddda = 0.0d0
        dpradddz = 0.0d0
        dpraddtt = 4.0d0 * asol * temp * temp
        dpraddta = 0.0d0
        dpraddtz = 0.0d0
        dpraddaa = 0.0d0
        dpraddaz = 0.0d0
        dpraddzz = 0.0d0

! energy in erg/gr
        erad    = 3.0d0 * prad * deninv

! first derivatives
        deraddd = -erad * deninv
        deraddt = 3.0d0 * dpraddt * deninv
        deradda = 0.0d0
        deraddz = 0.0d0

! second derivatives
        deradddd = -2.0d0 * deraddd * deninv
        deradddt = -deraddt * deninv
        deraddda = 0.0d0
        deradddz = 0.0d0
        deraddtt = 3.0d0 * dpraddtt * deninv
        deraddta = 0.0d0
        deraddtz = 0.0d0
        deraddaa = 0.0d0
        deraddaz = 0.0d0
        deraddzz = 0.0d0

! entropy in erg/g/kelvin
! from the first law e + P/rho = T*s + mu_rad * n_rad [erg/g/K]
! with zero chemical potential for photons
        srad    = (prad*deninv + erad) * tempinv

! first derivatives
        dsraddd = ((dpraddd - prad*deninv)*deninv + deraddd) * tempinv
        dsraddt = (dpraddt*deninv + deraddt - srad)  * tempinv
        dsradda = 0.0d0
        dsraddz = 0.0d0

! second derivatives
        dsradddd = ((dpradddd &
                    - (2.0d0*dpraddd - 2.0d0*prad*deninv)*deninv)*deninv &
                     + deradddd) * tempinv
        dsradddt = ((dpradddt - dpraddt*deninv)*deninv &
                    + deradddt - dsraddd) * tempinv
        dsraddda = 0.0d0
        dsradddz = 0.0d0
        dsraddtt = ((dpraddtt*deninv+deraddtt-dsraddt)-dsraddt)*tempinv
        dsraddta = 0.0d0
        dsraddtz = 0.0d0
        dsraddaa = 0.0d0
        dsraddaz = 0.0d0
        dsraddzz = 0.0d0

[fxt44]

where deninv = 1/rho and tempinv = 1/T (do those divisions only once).

[fxt44]

other obvious places this kind of issue might pop up in Skye and ideal.

anytime radiation dominates, many other quantities become nonsense because Prad, Erad, and Srad are large. for example, when ions or electrons set the correct derivatives they get completely swamped by the radiation terms.

[evbauer]

The fix in #591 seems to work nicely for $\chi_\rho$ and $c_P$ on ideal. Compare to plot at the top of this issue.

[fxt44]

bravo! left a PR comment on erad and srad.