diff --git a/CHANGELOG.md b/CHANGELOG.md index 465399f4cf..e8677ecf46 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface - [[PR373]](https://github.com/lanl/singularity-eos/pull/373) Initialize cache in `get_sg_eos*` functions +- [[PR374]](https://github.com/lanl/singularity-eos/pull/374) Make the Davis EOS more numerically robust ### Changed (changing behavior/API/variables/...) - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 0fb2217e0d..fb34c2a229 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -18,8 +18,10 @@ #include #include +#include #include #include +#include #include #include @@ -41,17 +43,21 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real es = Es(rho); - const Real tmp = std::pow((1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0, - 1.0 / (1.0 + _alpha)); - if (tmp > 0) return Ts(rho) * tmp; - return Ts(rho) + (sie - es) / _Cv0; // This branch is a negative temperature + const Real power_base = DimlessEdiff(rho, sie); + if (power_base <= 0) { + // This case would result in an imaginary temperature (i.e. negative), but we won't + // allow that so return zero + return 0.; + } + const Real tmp = std::pow(power_base, 1.0 / (1.0 + _alpha)); + return Ts(rho) * tmp; } template PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( const Real rho, const Real temp, Indexer_t &&lambda = static_cast(nullptr)) const { const Real t_s = Ts(rho); + PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); return Es(rho) + _Cv0 * t_s / (1.0 + _alpha) * (std::pow(temp / t_s, 1.0 + _alpha) - 1.0); } @@ -72,8 +78,11 @@ class DavisReactants : public EosBase { template PORTABLE_INLINE_FUNCTION Real MinInternalEnergyFromDensity( const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { - MinInternalEnergyIsNotEnabled("DavisReactants"); - return 0.0; + // Minimum enegy is when the returned temperature is zero. This only happens + // when the base to the exponent is zero (see T(rho, e) equation) + const Real es = Es(rho); + const Real ts = Ts(rho); + return es - (_Cv0 * ts) / (1 + _alpha); } template PORTABLE_INLINE_FUNCTION Real @@ -100,8 +109,12 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - return _Cv0 / std::pow((1 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1, - -_alpha / (1 + _alpha)); + const Real power_base = DimlessEdiff(rho, sie); + if (power_base <= 0) { + // Return zero heat capacity instead of an imaginary value + return 0.; + } + return _Cv0 / std::pow(power_base, -_alpha / (1 + _alpha)); } template PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( @@ -166,6 +179,7 @@ class DavisReactants : public EosBase { // static constexpr const char _eos_type[] = "DavisReactants"; static constexpr unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; + PORTABLE_FORCEINLINE_FUNCTION Real DimlessEdiff(const Real rho, const Real sie) const; PORTABLE_INLINE_FUNCTION Real Ps(const Real rho) const; PORTABLE_INLINE_FUNCTION Real Es(const Real rho) const; PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const; @@ -300,15 +314,24 @@ class DavisProducts : public EosBase { static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; PORTABLE_INLINE_FUNCTION Real F(const Real rho) const { + if (rho <= 0) { + return 0.; + } const Real vvc = 1.0 / (rho * _vc); return 2.0 * _a / (std::pow(vvc, 2 * _n) + 1.0); } PORTABLE_INLINE_FUNCTION Real Ps(const Real rho) const { + if (rho <= 0) { + return 0.; + } const Real vvc = 1 / (rho * _vc); return _pc * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) / std::pow(vvc, _k + _a) * (_k - 1.0 + F(rho)) / (_k - 1.0 + _a); } PORTABLE_INLINE_FUNCTION Real Es(const Real rho) const { + if (rho <= 0) { + return 0.; + } const Real vvc = 1 / (rho * _vc); const Real ec = _pc * _vc / (_k - 1.0 + _a); // const Real de = ecj-(Es(rho0)-_E0); @@ -316,6 +339,9 @@ class DavisProducts : public EosBase { std::pow(vvc, _k - 1.0 + _a); } PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const { + if (rho <= 0) { + return 0.; + } const Real vvc = 1 / (rho * _vc); return std::pow(2.0, -_a * _b / _n) * _pc * _vc / (_Cv * (_k - 1 + _a)) * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n * (1 - _b)) / @@ -326,9 +352,14 @@ class DavisProducts : public EosBase { } }; +PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, + const Real sie) const { + return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1.0; +} + PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; - const Real y = 1.0 - _rho0 / rho; + const Real y = 1.0 - robust::ratio(_rho0, std::max(rho, 0.)); const Real phat = 0.25 * _A * _A / _B * _rho0; const Real b4y = 4.0 * _B * y; @@ -342,7 +373,7 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { } } PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const { - const Real y = 1 - _rho0 / rho; + const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); const Real phat = 0.25 * _A * _A / _B * _rho0; const Real b4y = 4 * _B * y; Real e_s; @@ -354,19 +385,17 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const { } else { e_s = -y - (1.0 - std::exp(b4y)) / (4.0 * _B); } - return _e0 + _P0 * (1.0 / _rho0 - 1.0 / rho) + phat / _rho0 * e_s; + return _e0 + _P0 * (1.0 / _rho0 - robust::ratio(1.0, std::max(rho, 0.))) + + phat / _rho0 * e_s; } PORTABLE_INLINE_FUNCTION Real DavisReactants::Ts(const Real rho) const { - if (rho >= _rho0) { - const Real y = 1 - _rho0 / rho; - return _T0 * std::exp(-_Z * y) * std::pow(_rho0 / rho, -_G0 - _Z); - } else { - return _T0 * std::pow(_rho0 / rho, -_G0); - } + const Real rho0overrho = robust::ratio(_rho0, std::max(rho, 0.)); + const Real y = 1 - rho0overrho; + return _T0 * std::exp(-_Z * y) * std::pow(rho0overrho, -_G0 - _Z); } PORTABLE_INLINE_FUNCTION Real DavisReactants::Gamma(const Real rho) const { if (rho >= _rho0) { - const Real y = 1 - _rho0 / rho; + const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); return _G0 + _Z * y; } else { return _G0; @@ -377,11 +406,11 @@ template PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda) const { using namespace math_utils; - const Real y = 1 - _rho0 / rho; + const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); const Real phat = 0.25 * _A * _A / _B * _rho0; const Real b4y = 4 * _B * y; - const Real gamma = Gamma(rho); - const Real esv = -Ps(rho); + const Real gamma = Gamma(std::max(rho, 0.)); + const Real esv = -Ps(std::max(rho, 0.)); const Real psv = (rho >= _rho0) ? -phat * _rho0 * @@ -389,8 +418,10 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner 3 * y / pow<4>(1 - y) + 4 * pow<2>(y) / pow<5>(1 - y)) : -phat * 4 * _B * _rho0 * std::exp(b4y); const Real gammav = (rho >= _rho0) ? _Z * _rho0 : 0.0; - return -(psv + (sie - Es(rho)) * rho * (gammav - gamma * rho) - gamma * rho * esv) / - rho; + const Real numerator = + -(psv + (sie - Es(rho)) * std::max(rho, 0.) * (gammav - gamma * std::max(rho, 0.)) - + gamma * std::max(rho, 0.) * esv); + return robust::ratio(numerator, std::max(rho, 0.)); } template @@ -398,9 +429,10 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { // First, solve P=P(rho,T) for rho. Note P(rho,e) has an sie-es term, which is only a // function of T + PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { return (Ps(r) + Gamma(r) * r * _Cv0 * Ts(r) / (1 + _alpha) * - (std::pow(temp / Ts(r), 1 + _alpha) - 1.0)); + (std::pow(robust::ratio(temp, Ts(r)), 1 + _alpha) - 1.0)); }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; @@ -410,6 +442,10 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: " "Root find failed to find a solution given P, T\n"); } + if (rho < 0.) { + EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: " + "Root find resulted in a negative density\n"); + } sie = InternalEnergyFromDensityTemperature(rho, temp); } @@ -418,6 +454,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::FillEos(Real &rho, Real &temp, Rea Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda) const { + // BROKEN: This can only handle density-energy inputs! MUST BE CHANGED if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); if (output & thermalqs::temperature) temp = TemperatureFromDensityInternalEnergy(rho, sie); @@ -450,6 +487,9 @@ template PORTABLE_INLINE_FUNCTION Real DavisProducts::BulkModulusFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda) const { using namespace math_utils; + if (rho <= 0) { + return 0.; + } const Real vvc = 1 / (rho * _vc); const Real Fx = -4 * _a * std::pow(vvc, 2 * _n - 1) / pow<2>(1 + std::pow(vvc, 2 * _n)); const Real tmp = std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) / @@ -470,6 +510,7 @@ PORTABLE_INLINE_FUNCTION Real DavisProducts::BulkModulusFromDensityInternalEnerg template PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperature( const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { + PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { return (Ps(r) + Gamma(r) * r * _Cv * (temp - Ts(r))); }; @@ -482,12 +523,17 @@ PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperatur EOS_ERROR("DavisProducts::DensityEnergyFromPressureTemperature: " "Root find failed to find a solution given P, T\n"); } + if (rho < 0.) { + EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: " + "Root find resulted in a negative density\n"); + } sie = InternalEnergyFromDensityTemperature(rho, temp); } template PORTABLE_INLINE_FUNCTION void DavisProducts::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda) const { + // BROKEN: This can only handle density-energy inputs! MUST BE CHANGED if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); if (output & thermalqs::temperature) temp = TemperatureFromDensityInternalEnergy(rho, sie);