From 83a658b43f6d97d5e321ecc88e269fbd466af1ca Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 18:54:05 -0600 Subject: [PATCH 01/24] Add robustness to calculations --- singularity-eos/eos/eos_davis.hpp | 72 ++++++++++++++++++++++--------- 1 file changed, 52 insertions(+), 20 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 87e2ef06f6..1e6c7b6efb 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -18,6 +18,7 @@ #include #include +#include #include #include #include @@ -42,16 +43,21 @@ class DavisReactants : public EosBase { 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 = (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0; + 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 - (1 + _alpha) / (_Cv0 * ts); } 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 = (1 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1; + 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( @@ -149,7 +162,7 @@ class DavisReactants : public EosBase { static inline unsigned long scratch_size(std::string method, unsigned int nelements) { return 0; } - static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } + static inline unsigned long std::max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"DavisReactants Params: "}; printf("%srho0:%e e0:%e P0:%e\nT0:%e A:%e B:%e\nC:%e G0:%e Z:%e\nalpha:%e " @@ -283,7 +296,7 @@ class DavisProducts : public EosBase { static inline unsigned long scratch_size(std::string method, unsigned int nelements) { return 0; } - static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } + static inline unsigned long std::max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"DavisProducts Params: "}; printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e E0:%e\n", s1, _a, _b, _k, _n, _vc, @@ -300,15 +313,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 +338,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)) / @@ -328,7 +353,7 @@ class DavisProducts : public EosBase { 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, rho); const Real phat = 0.25 * _A * _A / _B * _rho0; const Real b4y = 4.0 * _B * y; @@ -342,7 +367,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 +379,19 @@ 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); + const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); + return _T0 * std::exp(-_Z * y) * std::pow(robust::ratio(_rho0, std::max(rho, 0.)), -_G0 - _Z); } else { - return _T0 * std::pow(_rho0 / rho, -_G0); + return _T0 * std::pow(robust::ratio(_rho0, rho), -_G0); } } 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,7 +402,7 @@ 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); @@ -389,8 +414,8 @@ 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; + return -(psv + (sie - Es(rho)) * rho * (gammav - gamma * std::max(rho, 0.)) - robust::ratio(gamma * std::max(rho, 0.) * esv), + rho); } template @@ -398,6 +423,7 @@ 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)); @@ -418,6 +444,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 +477,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 +500,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))); }; @@ -488,6 +519,7 @@ 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); From 3ce28ef4ea2f8698c6389b1e42d577fe630bbaa2 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 18:54:41 -0600 Subject: [PATCH 02/24] Clang format --- singularity-eos/eos/eos_davis.hpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 1e6c7b6efb..c767f186b1 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -379,12 +379,14 @@ 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 - robust::ratio(1.0, std::max(rho, 0.))) + 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 - robust::ratio(_rho0, std::max(rho, 0.)); - return _T0 * std::exp(-_Z * y) * std::pow(robust::ratio(_rho0, std::max(rho, 0.)), -_G0 - _Z); + return _T0 * std::exp(-_Z * y) * + std::pow(robust::ratio(_rho0, std::max(rho, 0.)), -_G0 - _Z); } else { return _T0 * std::pow(robust::ratio(_rho0, rho), -_G0); } @@ -414,8 +416,9 @@ 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 * std::max(rho, 0.)) - robust::ratio(gamma * std::max(rho, 0.) * esv), - rho); + return -(psv + (sie - Es(rho)) * rho * (gammav - gamma * std::max(rho, 0.)) - + robust::ratio(gamma * std::max(rho, 0.) * esv), + rho); } template From 02b089ee0cd84dac894cc2190817d294250e3df4 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 18:59:22 -0600 Subject: [PATCH 03/24] Add robust_utils.hpp header --- singularity-eos/eos/eos_davis.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index c767f186b1..093125e322 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -22,6 +22,7 @@ #include #include #include +#include #include namespace singularity { From b50971ebe41972387ff2e5d2bf56b715149c08ce Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 19:05:54 -0600 Subject: [PATCH 04/24] Update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6800600af4..12fbe42a13 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,7 @@ - [[PR335]](https://github.com/lanl/singularity-eos/pull/335) Fix missing hermite.hpp in CMake install required for Helmholtz EOS - [[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 +- [[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 From 1feb551c46fdd3781a7aca4e6203dc4f598ea6da Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 19:06:38 -0600 Subject: [PATCH 05/24] Clang format --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 093125e322..80482da087 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -21,8 +21,8 @@ #include #include #include -#include #include +#include #include namespace singularity { From d1714e421c55d945bd7e818df6e800574bf7885c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 20:45:54 -0600 Subject: [PATCH 06/24] Fix some silly mistakes --- singularity-eos/eos/eos_davis.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 80482da087..d3cb8e0eb0 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -163,7 +163,7 @@ class DavisReactants : public EosBase { static inline unsigned long scratch_size(std::string method, unsigned int nelements) { return 0; } - static inline unsigned long std::max_scratch_size(unsigned int nelements) { return 0; } + static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"DavisReactants Params: "}; printf("%srho0:%e e0:%e P0:%e\nT0:%e A:%e B:%e\nC:%e G0:%e Z:%e\nalpha:%e " @@ -417,9 +417,9 @@ 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 * std::max(rho, 0.)) - - robust::ratio(gamma * std::max(rho, 0.) * esv), - rho); + const Real numerator = -(psv + (sie - Es(rho)) * rho * (gammav - gamma * + std::max(rho, 0.)) - gamma * std::max(rho, 0.) * esv); + return robust::ratio(numerator, rho); } template From 62fc6824a4a862c6613787913dd35671d2a0559e Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 20:46:23 -0600 Subject: [PATCH 07/24] Clang format --- singularity-eos/eos/eos_davis.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index d3cb8e0eb0..d6f382b1c4 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -417,8 +417,9 @@ 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; - const Real numerator = -(psv + (sie - Es(rho)) * rho * (gammav - gamma * - std::max(rho, 0.)) - gamma * std::max(rho, 0.) * esv); + const Real numerator = + -(psv + (sie - Es(rho)) * rho * (gammav - gamma * std::max(rho, 0.)) - + gamma * std::max(rho, 0.) * esv); return robust::ratio(numerator, rho); } From 4d78e474d5cf1dde7f36b6932e3baa75feeb7a47 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 20:48:24 -0600 Subject: [PATCH 08/24] Another silly find/replace --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index d6f382b1c4..b66d1dfd13 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -297,7 +297,7 @@ class DavisProducts : public EosBase { static inline unsigned long scratch_size(std::string method, unsigned int nelements) { return 0; } - static inline unsigned long std::max_scratch_size(unsigned int nelements) { return 0; } + static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"DavisProducts Params: "}; printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e E0:%e\n", s1, _a, _b, _k, _n, _vc, From ac2a8bcf21b2256ad19c9900bb661cb191f76d34 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 20:53:28 -0600 Subject: [PATCH 09/24] Last silly mistake hopefully... --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index b66d1dfd13..027138c8da 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -428,7 +428,7 @@ 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") + 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)); From 61edd15006d3b4a6242ee715d96e188b0e57e3d7 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:32:33 -0600 Subject: [PATCH 10/24] Fix missing maxes and remove last negative density branch --- singularity-eos/eos/eos_davis.hpp | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 027138c8da..8dc257272d 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -83,7 +83,7 @@ class DavisReactants : public EosBase { // 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 - (1 + _alpha) / (_Cv0 * ts); + return es - (_Cv0 * ts) / (1 + _alpha); } template PORTABLE_INLINE_FUNCTION Real @@ -354,7 +354,7 @@ class DavisProducts : public EosBase { PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; - const Real y = 1.0 - robust::ratio(_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; @@ -384,13 +384,9 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const { phat / _rho0 * e_s; } PORTABLE_INLINE_FUNCTION Real DavisReactants::Ts(const Real rho) const { - if (rho >= _rho0) { - const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); - return _T0 * std::exp(-_Z * y) * - std::pow(robust::ratio(_rho0, std::max(rho, 0.)), -_G0 - _Z); - } else { - return _T0 * std::pow(robust::ratio(_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) { From 4ec06fec56b4fad04abe495bf5c149a3588cb3e3 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:34:08 -0600 Subject: [PATCH 11/24] Add missing maxes and do a clang format --- singularity-eos/eos/eos_davis.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 8dc257272d..e80aea9839 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -404,8 +404,8 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner 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 * @@ -414,9 +414,9 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner : -phat * 4 * _B * _rho0 * std::exp(b4y); const Real gammav = (rho >= _rho0) ? _Z * _rho0 : 0.0; const Real numerator = - -(psv + (sie - Es(rho)) * rho * (gammav - gamma * std::max(rho, 0.)) - + -(psv + (sie - Es(std::max(rho, 0.))) * rho * (gammav - gamma * std::max(rho, 0.)) - gamma * std::max(rho, 0.) * esv); - return robust::ratio(numerator, rho); + return robust::ratio(numerator, std::max(rho, 0.)); } template From 4e2c4b469440133c5e023fbcc65e71362f2c2c76 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:44:34 -0600 Subject: [PATCH 12/24] Break out dimensionless energy into its own function --- singularity-eos/eos/eos_davis.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index e80aea9839..ca82e25824 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -44,7 +44,7 @@ class DavisReactants : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { const Real es = Es(rho); - const Real power_base = (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0; + const Real power_base = DimlessEdiff(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 @@ -110,7 +110,7 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real power_base = (1 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1; + const Real power_base = DimlessEdiff(sie); if (power_base <= 0) { // Return zero heat capacity instead of an imaginary value return 0.; @@ -180,6 +180,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 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; @@ -352,6 +353,10 @@ class DavisProducts : public EosBase { } }; +PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const real sie) const { + return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0 +} + PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; const Real y = 1.0 - robust::ratio(_rho0, std::max(rho, 0.)); From ad3cb40176631ea3a3f939810db62a9a1858972e Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:44:55 -0600 Subject: [PATCH 13/24] Clang format --- singularity-eos/eos/eos_davis.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index ca82e25824..19d5cb82c4 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -354,8 +354,7 @@ class DavisProducts : public EosBase { }; PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const real sie) const { - return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0 -} + return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0} PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; From 88f7b861428281a6c4e59e75aa7f60880665fba0 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:56:17 -0600 Subject: [PATCH 14/24] Add max to another rho --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 19d5cb82c4..c7f70dd974 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -418,7 +418,7 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner : -phat * 4 * _B * _rho0 * std::exp(b4y); const Real gammav = (rho >= _rho0) ? _Z * _rho0 : 0.0; const Real numerator = - -(psv + (sie - Es(std::max(rho, 0.))) * rho * (gammav - gamma * std::max(rho, 0.)) - + -(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.)); } From c68391ba011145ac3025cbcba46be5045dd13fda Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:03:14 -0600 Subject: [PATCH 15/24] A few more robustness changes --- singularity-eos/eos/eos_davis.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index c7f70dd974..ace6dfb885 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -430,8 +430,8 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu // 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)); + return (Ps(r) + Gamma(r) * std::max(r, 0.) * _Cv0 * Ts(r) / (1 + _alpha) * + (std::pow(robust::ratio(temp, Ts(r)), 1 + _alpha) - 1.0)); }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; @@ -507,7 +507,7 @@ PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperatur 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))); + return (Ps(r) + Gamma(r) * std::max(r, 0.) * _Cv * (temp - Ts(r))); }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; From 5a84cab4fcaae7e2a862b88660d2191a04eca218 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:03:38 -0600 Subject: [PATCH 16/24] Clang format --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index ace6dfb885..b23db2358b 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -431,7 +431,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { return (Ps(r) + Gamma(r) * std::max(r, 0.) * _Cv0 * Ts(r) / (1 + _alpha) * - (std::pow(robust::ratio(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; From 43fb3ad7c8959ba9e183c13dadb0e8865e831b6a Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:07:55 -0600 Subject: [PATCH 17/24] Correct real to Real and switch to error instead of maxing density --- singularity-eos/eos/eos_davis.hpp | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index b23db2358b..87f843a649 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -353,7 +353,7 @@ class DavisProducts : public EosBase { } }; -PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const real sie) const { +PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real sie) const { return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0} PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { @@ -430,7 +430,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu // function of T PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { - return (Ps(r) + Gamma(r) * std::max(r, 0.) * _Cv0 * Ts(r) / (1 + _alpha) * + return (Ps(r) + Gamma(r) * r * _Cv0 * Ts(r) / (1 + _alpha) * (std::pow(robust::ratio(temp, Ts(r)), 1 + _alpha) - 1.0)); }; using RootFinding1D::regula_falsi; @@ -441,6 +441,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") + } sie = InternalEnergyFromDensityTemperature(rho, temp); } @@ -507,7 +511,7 @@ PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperatur 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) * std::max(r, 0.) * _Cv * (temp - Ts(r))); + return (Ps(r) + Gamma(r) * r * _Cv * (temp - Ts(r))); }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; @@ -518,6 +522,10 @@ 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") + } sie = InternalEnergyFromDensityTemperature(rho, temp); } template From e30a855cc18f18bd3214861dc49c9b93a7ad16dc Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:12:47 -0600 Subject: [PATCH 18/24] Add missing density to DimlessEDiff function --- singularity-eos/eos/eos_davis.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 87f843a649..0d3281ee51 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -44,7 +44,7 @@ class DavisReactants : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { const Real es = Es(rho); - const Real power_base = DimlessEdiff(sie); + 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 @@ -110,7 +110,7 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real power_base = DimlessEdiff(sie); + const Real power_base = DimlessEdiff(rho, sie); if (power_base <= 0) { // Return zero heat capacity instead of an imaginary value return 0.; @@ -180,7 +180,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 sie) const; + 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; @@ -353,7 +353,7 @@ class DavisProducts : public EosBase { } }; -PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real sie) const { +PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, const Real sie) const { return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0} PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { From 98541a73e22b5c9aedaac82575e29dcd03b03c79 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:13:22 -0600 Subject: [PATCH 19/24] Clang format --- singularity-eos/eos/eos_davis.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 0d3281ee51..d0813317fc 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -353,7 +353,8 @@ class DavisProducts : public EosBase { } }; -PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, const Real sie) const { +PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, + const Real sie) const { return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0} PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { From 5377650864e40a8b8add4eec7bb29180b0d5eae5 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:15:08 -0600 Subject: [PATCH 20/24] Correct es to Es(rho) --- singularity-eos/eos/eos_davis.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index d0813317fc..bd4cdc1742 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -43,7 +43,6 @@ 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 power_base = DimlessEdiff(rho, sie); if (power_base <= 0) { // This case would result in an imaginary temperature (i.e. negative), but we won't @@ -355,7 +354,7 @@ 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) + 1.0} + 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; From 4f380404272075749d969e01bc9661de71ee09c4 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:20:12 -0600 Subject: [PATCH 21/24] Add missing semicolon --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index bab74467df..137e74542a 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -354,7 +354,7 @@ 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} + 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; From 6708ddada5538e6431e4711c4f6c56dc31e3a58c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:22:16 -0600 Subject: [PATCH 22/24] Add missing bracket --- singularity-eos/eos/eos_davis.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 137e74542a..d53b4b8818 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -354,7 +354,8 @@ 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}; + 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; From c422370b535963ad3ee239420a43f41242b96352 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:23:06 -0600 Subject: [PATCH 23/24] Clang format --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index d53b4b8818..2c2987ee6e 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -354,7 +354,7 @@ 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; + return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1.0; } PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { From ff3901dd0d301607739b9a52ff79361191f4fdf5 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:25:28 -0600 Subject: [PATCH 24/24] Add missing semicolon --- singularity-eos/eos/eos_davis.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 2c2987ee6e..fb34c2a229 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -444,7 +444,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu } if (rho < 0.) { EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: " - "Root find resulted in a negative density") + "Root find resulted in a negative density\n"); } sie = InternalEnergyFromDensityTemperature(rho, temp); } @@ -525,7 +525,7 @@ PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperatur } if (rho < 0.) { EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: " - "Root find resulted in a negative density") + "Root find resulted in a negative density\n"); } sie = InternalEnergyFromDensityTemperature(rho, temp); }