Fix bug in Gruneisen DensityEnergyFromPressureTemperature()

[jhp-lanl]

PR Summary

@Gerbeset wasn't able to use the DensityEnergyFromPressureTemperature() function with ambient conditions and properties for aluminum. This pointed to a bug in the implementation.

As far as I can tell, the bug in question was in the bounds for the root find. The big issue I think was that the upper bound was outside the applicable maximum density of the EOS. Now that we have a maximum density implemented, this can bound the root find appropriately.

EDIT - It turns out that the issue was thermodynamic consistency. For a given temperature, there is a maximum in the pressure since the cold curve energy doesn't increase with increasing density. As a result, at temperatures/energies below the reference curve, there will be a regime where the pressure decreases as a function of density. This is not a phase change, but rather a result of thermodynamic consistency breaking. The root find assumed a monotonic $P(rho, T)$ curve for a given temperature, so it was failing to converge.

See above for a plot of the pressure for various isotherms. This should be monotonic.

This MR introduces a new function that bounds the density to a range where the EOS is thermodynamically stable when performing the P-T lookup.

An open question is whether to apply this domain limiting procedure to all lookups or not.

Less important changes

I also switched from a reference capture lambda to our PORABLE_LAMBDA.

To make the implementation easier to follow (and more general), I decided to just use density-temperature lookups everywhere. As a result, I find the density that satisfies the pressure-temperature point first and then use that density in the energy lookup. This isn't strictly necessary since the EOS assumes that energy is constant along an isotherm in density space (violating thermodynamic consistency), but I figured I'd do it this way anyway to be more correct. I think there should be negligible performance difference if any. Since this is generally an initialization function, I doubt that matters.

I also added a test to make sure that the call is returning correct values.

PR Checklist

• Adds a test for any bugs fixed. Adds tests for new features.
• Format your changes by using the make format command after configuring with cmake.
• N/A Document any new features, update documentation for changes made.
• Make sure the copyright notice on any files you modified is up to date.
• After creating a pull request, note it in the CHANGELOG.md file.
• LANL employees: make sure tests pass both on the github CI and on the Darwin CI

[jhp-lanl]

@chadmeyer and @jonahm-LANL can you review?

[jhp-lanl]

Found another bug... working on it. Will make draft for now

[jhp-lanl]

@chadmeyer and @Yurlungur I had to make a few more changes.

It looks like the primary issue is that the Gruneisen EOS pressure is not monotonic along an isotherm. I didn't realize that subtraction of the reference energy curve will eventually overwhelms the pressure reference curve, causing a maximum in the pressure.

I introduced a function to calculate the location of this maximum and essentially bound the EOS to the region up to this maximum (in P-T space). We could think about imposing this maximum on the other lookups too, but it would add an expensive root find. The P-T lookup is most likely at initialization so it's okay to be a bit more expensive there.

[Yurlungur]

Good catch! I do think we should try to limit root finding generically in equations of state. However, it probably suggests each EOS should have additional introspection.

Before we merge I would like Chad's eyes on this. @chadmeyer

[jhp-lanl]

I'm converting back to draft. I'm encountering a situation where dPres_drho_e is returning a negative value and I'm not quite sure why.

[jhp-lanl]

I'm converting back to draft. I'm encountering a situation where dPres_drho_e is returning a negative value and I'm not quite sure why.

More specifically, the reference pressure in the beginning of the algorithm is of order -1e-13. Something is very weird here maybe in how the EOS parameters are being set up.

[jhp-lanl]

I'm converting back to draft. I'm encountering a situation where dPres_drho_e is returning a negative value and I'm not quite sure why.

More specifically, the reference pressure in the beginning of the algorithm is of order -1e-13. Something is very weird here maybe in how the EOS parameters are being set up.

I figured out the issue... it was a unit conversion issue with the input specification. I think this MR is fine.

[Yurlungur]

@chadmeyer please take a look

[chadmeyer]

Good catch! I do think we should try to limit root finding generically in equations of state. However, it probably suggests each EOS should have additional introspection.

I will observe here that RE of PT is generally going to be a root find for most EOSs, and isn't meant to be called much outside of certain special cases (i.e. not in an inner loop outside of initialization/remapping of an EOS).

[chadmeyer]

I think this seems like a reasonable change, given the issue seen. Is there a more thermodynamically consistent form of the Gruneisen we could implement at some point?

[jhp-lanl]

I think this seems like a reasonable change, given the issue seen. Is there a more thermodynamically consistent form of the Gruneisen we could implement at some point?

@chadmeyer , we actually already have one that @aematts provided in eos_mgusup.hpp (documentation). The only downside of that EOS is that it is limited to $b = 0, S_2 = 0, S_3 = 0$, referencing the parameters in this EOS.

Ann has also provided a power law Mie-Gruneisen that offers more flexibility in the shape of the Hugoniot reference curve, but could not use the Steinberg parameters directly. I think this version also has restrictions on the functional form of the Gruneisen parameter, possibly limiting accuracy off of the principal Hugoniot (re-shocks or rarefactions).

More generally though, I think the only way to make this EOS thermodynamically consistent is to numerically integrate a temperature reference curve (i.e. the temperature along the Hugoniot reference curve). You could imagine doing this a priori and then either fitting it to a function or tabulating it (and maybe storing it in a spiner table). You could also do this in-line, but that would probably be expensive.

[Yurlungur]

It's worth pointing out that this Gruneisen is specifically in the code because it's not thermodynamically consistent but it is consistent with several host code implementations.

[aematts]

I think this seems like a reasonable change, given the issue seen. Is there a more thermodynamically consistent form of the Gruneisen we could implement at some point?

@chadmeyer , we actually already have one that @aematts provided in eos_mgusup.hpp (documentation). The only downside of that EOS is that it is limited to b = 0 , S 2 = 0 , S 3 = 0 , referencing the parameters in this EOS.

Ann has also provided a power law Mie-Gruneisen that offers more flexibility in the shape of the Hugoniot reference curve, but could not use the Steinberg parameters directly. I think this version also has restrictions on the functional form of the Gruneisen parameter, possibly limiting accuracy off of the principal Hugoniot (re-shocks or rarefactions).

More generally though, I think the only way to make this EOS thermodynamically consistent is to numerically integrate a temperature reference curve (i.e. the temperature along the Hugoniot reference curve). You could imagine doing this a priori and then either fitting it to a function or tabulating it (and maybe storing it in a spiner table). You could also do this in-line, but that would probably be expensive.

@jhp-lanl , @chadmeyer , @Yurlungur : It is true that all three of the EOS's I have implemented uses $\Gamma \rho = const$, that is $b_0=0$. It is mainly because then it is fairly easy to derive an entropy for a MG type EOS. I can think about how to modify my derivation to use a Gruneisen type $\Gamma$ relation.

I also want to mention the Extended Vinet EOS that is based on an isotherm (any temperature) and might be better for off-Hugoniot cases. It is a full EOS and not the cold curve form that is commonly named Vinet at LANL.

[jhp-lanl]

I also want to mention the Extended Vinet EOS that is based on an isotherm (any temperature) and might be better for off-Hugoniot cases. It is a full EOS and not the cold curve form that is commonly named Vinet at LANL.

I completely agree. Hugoniot reference curves are in general less flexible, and if calibrations exist, an isothermal reference curve EOS is preferable.

[jhp-lanl]

One last edge case I found: if the reference temperature is absurd (i.e. a unit conversion issue), it's possible for the maximum pressure to occur at a density less than the reference density. In this case, all compression states are unstable, which is bogus. I added an error to help with this case.

Will merge when tests pass.