diff --git a/CHANGELOG.md b/CHANGELOG.md index 6800600af4..e8677ecf46 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,9 +20,12 @@ - [[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 +- [[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 +- [[PR372]](https://github.com/lanl/singularity-eos/pull/372) Removed E0 from Davis Products EOS in favor of using the shifted EOS modifier. CHANGES API! ### Infrastructure (changes irrelevant to downstream codes) - [[PR329]](https://github.com/lanl/singularity-eos/pull/329) Move vinet tests into analytic test suite diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 12f555caac..ca42caed77 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1254,17 +1254,44 @@ Here, there are four dimensionless parameters that are settable by the user, :math:`e_\mathrm{C}`, :math:`V_\mathrm{C}` and :math:`T_\mathrm{C}` are tuning parameters with units related to their non-subscripted counterparts. +Note that the energy zero (i.e. the reference energy) for the Davis products EOS +is arbitrary. For the isentrope to properly pass through the CJ state of a +reacting material, the energy release of the reaction needs to be accounted for +properly. If done external to the EOS, an energy source term is required in the +Euler equations. However, a common convention is to specify the reactants and +product EOS in a consistent way such that the reference energy corresponds to +the rest state of the material *before* it reacts. + +The energy at the CJ state can be calculated as + +.. math:: + + e_\mathrm{CJ} = \frac{P_0 + P_\mathrm{CJ}}{2(V_0 - V_\mathrm{CJ})}, + +relative to :math:`e = 0` at the reference state of the *reactants*. Therefore +the energy offset of the products EOS is given by + +.. math:: + + e_0 = e_S(V_\mathrm{CJ}) - e_\mathrm{CJ}. + +Practically, this means :math:`e_0` should be positive for any energetic material. + +To provide the energy offset to the Davis Products EOS, `the energy shift +modifier`_ should be used. Note that the convention there +is that the shift is positive, so :math:`-e_0` should be provided to the shift +modifier. + The constructor for the Davis Products EOS is .. code-block:: cpp DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc, - const Real pc, const Real Cv, const Real E0) + const Real pc, const Real Cv) where ``a`` is :math:`a`, ``b`` is :math:`b`, ``k`` is :math:`k`, ``n`` is :math:`n`, ``vc`` is :math:`V_\mathrm{C}`, ``pc`` is -:math:`P_\mathrm{C}`, ``Cv`` is :math:`C_{V,0}`, and ``E0`` is -:math:`e_\mathrm{C}`. +:math:`P_\mathrm{C}`, ``Cv`` is :math:`C_{V,0}`. Spiner EOS ```````````` diff --git a/doc/sphinx/src/modifiers.rst b/doc/sphinx/src/modifiers.rst index 552b16940b..75e11310af 100644 --- a/doc/sphinx/src/modifiers.rst +++ b/doc/sphinx/src/modifiers.rst @@ -22,6 +22,8 @@ We list below the available modifiers and their constructors. The Shifted EOS ----------------- +.. _modifiers shifted EOS: + The shifted equation of state modifies zero point energy of an underlying model by some shift. So for example, it transforms diff --git a/python/module.cpp b/python/module.cpp index a6b3ad1c79..2b89795552 100644 --- a/python/module.cpp +++ b/python/module.cpp @@ -60,9 +60,9 @@ PYBIND11_MODULE(singularity_eos, m) { eos_class(m, "DavisProducts") .def(py::init()) .def( - py::init(), + py::init(), py::arg("a"), py::arg("b"), py::arg("k"), py::arg("n"), py::arg("vc"), - py::arg("pc"), py::arg("Cv"), py::arg("E0") + py::arg("pc"), py::arg("Cv") ); #ifdef SPINER_USE_HDF diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index cd3c094f8a..d1ca6f89db 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -182,6 +182,7 @@ class PTESolverBase { void Finalize() { for (int m = 0; m < nmat; m++) { temp[m] *= Tnorm; + PORTABLE_REQUIRE(temp[m] >= 0., "Non-positive temperature returned"); u[m] *= uscale; press[m] *= uscale; } diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 87e2ef06f6..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; @@ -177,8 +191,8 @@ class DavisProducts : public EosBase { DavisProducts() = default; PORTABLE_INLINE_FUNCTION DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc, - const Real pc, const Real Cv, const Real E0) - : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv), _E0(E0) {} + const Real pc, const Real Cv) + : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv) {} DavisProducts GetOnDevice() { return *this; } template PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( @@ -286,8 +300,8 @@ class DavisProducts : public EosBase { 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, - _pc, _Cv, _E0); + printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e \n", s1, _a, _b, _k, _n, _vc, _pc, + _Cv); } inline void Finalize() {} static std::string EosType() { return std::string("DavisProducts"); } @@ -295,20 +309,29 @@ class DavisProducts : public EosBase { private: static constexpr Real onethird = 1.0 / 3.0; - Real _a, _b, _k, _n, _vc, _pc, _Cv, _E0; + Real _a, _b, _k, _n, _vc, _pc, _Cv; // static constexpr const char _eos_type[] = "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); @@ -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); diff --git a/singularity-eos/eos/get_sg_eos_p_t.cpp b/singularity-eos/eos/get_sg_eos_p_t.cpp index fb1db64d2d..5845a7f634 100644 --- a/singularity-eos/eos/get_sg_eos_p_t.cpp +++ b/singularity-eos/eos/get_sg_eos_p_t.cpp @@ -39,6 +39,10 @@ void get_sg_eos_p_t(const char *name, int ncell, int nmat, indirection_v &offset // for small loops const int32_t token{tokens.acquire()}; const int32_t tid{small_loop ? iloop : token}; + // need to initialize the scratch before it's used to avoid undefined behavior + for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { + solver_scratch(tid, idx) = 0.0; + } // caching mechanism singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0)); double mass_sum{0.0}; diff --git a/singularity-eos/eos/get_sg_eos_rho_e.cpp b/singularity-eos/eos/get_sg_eos_rho_e.cpp index 43049e31a7..95cb97b2ed 100644 --- a/singularity-eos/eos/get_sg_eos_rho_e.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_e.cpp @@ -39,6 +39,10 @@ void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v, int npte{0}; // initialize values for solver / lookup i_func(i, tid, mass_sum, npte, 0.0, 1.0, 0.0); + // need to initialize the scratch before it's used to avoid undefined behavior + for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { + solver_scratch(tid, idx) = 0.0; + } // get cache from offsets into scratch const int neq = npte + 1; singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) + diff --git a/singularity-eos/eos/get_sg_eos_rho_p.cpp b/singularity-eos/eos/get_sg_eos_rho_p.cpp index dc6da80ecc..4a788b1897 100644 --- a/singularity-eos/eos/get_sg_eos_rho_p.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_p.cpp @@ -41,6 +41,10 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v, // initialize values for solver / lookup i_func(i, tid, mass_sum, npte, 0.0, 0.0, 1.0); Real sie_tot_true{0.0}; + // need to initialize the scratch before it's used to avoid undefined behavior + for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { + solver_scratch(tid, idx) = 0.0; + } const int neq = npte + 1; singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) + neq * (neq + 4) + 2 * npte); diff --git a/singularity-eos/eos/get_sg_eos_rho_t.cpp b/singularity-eos/eos/get_sg_eos_rho_t.cpp index f1b91eed1f..7d79bc1b2c 100644 --- a/singularity-eos/eos/get_sg_eos_rho_t.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_t.cpp @@ -42,6 +42,10 @@ void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v, i_func(i, tid, mass_sum, npte, 1.0, 0.0, 0.0); // calculate pte condition (lookup for 1 mat cell) Real sie_tot_true{0.0}; + // need to initialize the scratch before it's used to avoid undefined behavior + for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { + solver_scratch(tid, idx) = 0.0; + } const int neq = npte; singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) + neq * (neq + 4) + 2 * npte); diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index e000bfe268..c80ab710b5 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -154,23 +154,23 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, - const double pc, const double Cv, const double E0, - int const *const enabled, double *const vals) { + const double pc, const double Cv, int const *const enabled, + double *const vals) { assert(matindex >= 0); - EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); + EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv)); if (enabled[3] == 1) { singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } - EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); + EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv)); eos[matindex] = eos_.GetOnDevice(); return 0; } int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, - const double pc, const double Cv, const double E0) { - return init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, E0, def_en, def_v); + const double pc, const double Cv) { + return init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, def_en, def_v); } int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index 2c0ef227ec..e325601758 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -112,13 +112,13 @@ end function init_sg_JWL interface integer(kind=c_int) function & - init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, E0, & + init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, & sg_mods_enabled, sg_mods_values) & bind(C, name='init_sg_DavisProducts') import integer(c_int), value, intent(in) :: matindex type(c_ptr), value, intent(in) :: eos - real(kind=c_double), value, intent(in) :: a, b, k, n, vc, pc, Cv, E0 + real(kind=c_double), value, intent(in) :: a, b, k, n, vc, pc, Cv type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values end function init_sg_DavisProducts end interface @@ -460,12 +460,12 @@ integer function init_sg_JWL_f(matindex, eos, A, B, R1, R2, w, rho0, Cv, & end function init_sg_JWL_f integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & - Cv, E0, sg_mods_enabled, & + Cv, sg_mods_enabled, & sg_mods_values) & result(err) integer(c_int), value, intent(in) :: matindex type(sg_eos_ary_t), intent(in) :: eos - real(kind=8), value, intent(in) :: a, b, k, n, vc, pc, Cv, E0 + real(kind=8), value, intent(in) :: a, b, k, n, vc, pc, Cv integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars @@ -478,7 +478,7 @@ integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values err = init_sg_DavisProducts(matindex-1, eos%ptr, a, b, k, n, vc, pc, Cv, & - E0, c_loc(sg_mods_enabled_use), & + c_loc(sg_mods_enabled_use), & c_loc(sg_mods_values_use)) end function init_sg_DavisProducts_f diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index 9a7ddd43f7..9a31d10f54 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -39,8 +39,8 @@ int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const doubl int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, - const double pc, const double Cv, const double E0, - int const *const enabled, double *const vals); + const double pc, const double Cv, int const *const enabled, + double *const vals); int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, const double e0, const double P0, const double T0, @@ -134,7 +134,7 @@ int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const doubl int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, - const double pc, const double Cv, const double E0); + const double pc, const double Cv); int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, const double e0, const double P0, const double T0, diff --git a/test/pte_test_utils.hpp b/test/pte_test_utils.hpp index d61b09747e..af2e68ad94 100644 --- a/test/pte_test_utils.hpp +++ b/test/pte_test_utils.hpp @@ -59,7 +59,7 @@ inline void set_eos(T *eos) { singularity::EOS dr = singularity::DavisReactants( 1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, 0.56, 0.0, 0.4265, 0.001074e10); singularity::EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419, - 3.2e10, 0.001072e10, 0.0); + 3.2e10, 0.001072e10); eos[0] = gr.GetOnDevice(); eos[1] = dr.GetOnDevice(); eos[2] = dp.GetOnDevice(); diff --git a/test/python_bindings.py b/test/python_bindings.py index a0ff479ad7..951e90523e 100644 --- a/test/python_bindings.py +++ b/test/python_bindings.py @@ -87,7 +87,7 @@ def testDavisReactants(self): eos = singularity_eos.DavisReactants(1,1,1,1,1,1,1,1,1,1,1) def testDavisProducts(self): - eos = singularity_eos.DavisProducts(1,1,1,1,1,1,1,1) + eos = singularity_eos.DavisProducts(1,1,1,1,1,1,1) class Modifiers(unittest.TestCase, EOSTestBase): def setUp(self): diff --git a/test/test_f_iface.f90 b/test/test_f_iface.f90 index 52a92e95cc..7b845a2119 100644 --- a/test/test_f_iface.f90 +++ b/test/test_f_iface.f90 @@ -43,7 +43,7 @@ program test_sg_fortran_interface 4.6d0, 0.34d0, 0.56d0,0.d0, 0.4265d0, 0.001074d10) mat = mat + 1 res = init_sg_DavisProducts_f(mat, eos, 0.798311d0, 0.58d0, 1.35d0, 2.66182d0,& - 0.75419d0, 3.2d10, 0.001072d10, 0.d0) + 0.75419d0, 3.2d10, 0.001072d10) ! cleanup res = finalize_sg_eos_f(nmat, eos)