Skip to content

Commit

Permalink
Merge pull request #374 from lanl/jhp/DavisRobust
Browse files Browse the repository at this point in the history 
Various robustness additions to the Davis EOS
  • Loading branch information
@jhp-lanl
jhp-lanl authored May 14, 2024
2 parents a213659 + ff3901d commit 5229919
Show file tree
Hide file tree
Showing 2 changed files with 72 additions and 25 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
96 changes: 71 additions & 25 deletions singularity-eos/eos/eos_davis.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,10 @@
#include <cmath>
#include <cstdio>

#include <ports-of-call/portable_errors.hpp>
#include <singularity-eos/base/constants.hpp>
#include <singularity-eos/base/math_utils.hpp>
#include <singularity-eos/base/robust_utils.hpp>
#include <singularity-eos/base/root-finding-1d/root_finding.hpp>
#include <singularity-eos/eos/eos_base.hpp>

Expand All @@ -41,17 +43,21 @@ class DavisReactants : public EosBase<DavisReactants> {
PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
const Real rho, const Real sie,
Indexer_t &&lambda = static_cast<Real *>(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 <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(
const Real rho, const Real temp,
Indexer_t &&lambda = static_cast<Real *>(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);
}
Expand All @@ -72,8 +78,11 @@ class DavisReactants : public EosBase<DavisReactants> {
template <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION Real MinInternalEnergyFromDensity(
const Real rho, Indexer_t &&lambda = static_cast<Real *>(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 <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION Real
Expand All @@ -100,8 +109,12 @@ class DavisReactants : public EosBase<DavisReactants> {
PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(
const Real rho, const Real sie,
Indexer_t &&lambda = static_cast<Real *>(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 <typename Indexer_t = Real *>
PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature(
Expand Down Expand Up @@ -166,6 +179,7 @@ class DavisReactants : public EosBase<DavisReactants> {
// 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;
Expand Down Expand Up @@ -300,22 +314,34 @@ class DavisProducts : public EosBase<DavisProducts> {
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);
return ec * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) /
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)) /
Expand All @@ -326,9 +352,14 @@ class DavisProducts : public EosBase<DavisProducts> {
}
};

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;

Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -377,30 +406,33 @@ template <typename Indexer_t>
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 *
(4 * _B * (1 + b4y + 0.5 * (pow<2>(b4y) + _C / 3 * pow<3>(b4y))) +
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 <typename Indexer_t>
PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperature(
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;
Expand All @@ -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);
}

Expand All @@ -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);
Expand Down Expand Up @@ -450,6 +487,9 @@ template <typename Indexer_t>
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) /
Expand All @@ -470,6 +510,7 @@ PORTABLE_INLINE_FUNCTION Real DavisProducts::BulkModulusFromDensityInternalEnerg
template <typename Indexer_t>
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)));
};
Expand All @@ -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 <typename Indexer_t>
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);
Expand Down

0 comments on commit 5229919

Please sign in to comment.