diff --git a/.github/workflows/tests_minimal.yml b/.github/workflows/tests_minimal.yml index 781802b3e9f..52f2e70e361 100644 --- a/.github/workflows/tests_minimal.yml +++ b/.github/workflows/tests_minimal.yml @@ -36,4 +36,4 @@ jobs: .. make -j4 make install - make test + ctest --output-on-failure diff --git a/.github/workflows/tests_minimal_kokkos.yml b/.github/workflows/tests_minimal_kokkos.yml index b9666138452..5038f44b4b1 100644 --- a/.github/workflows/tests_minimal_kokkos.yml +++ b/.github/workflows/tests_minimal_kokkos.yml @@ -36,4 +36,4 @@ jobs: .. make -j4 make install - make test + ctest --output-on-failure diff --git a/CHANGELOG.md b/CHANGELOG.md index 7e75ccc12dc..cffd00fc4f3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ - [[PR444]](https://github.com/lanl/singularity-eos/pull/444) Add Z split modifier and electron ideal gas EOS ### Fixed (Repair bugs, etc) +- [[PR449]](https://github.com/lanl/singularity-eos/pull/449) Ensure that DensityEnergyFromPressureTemperature works for all equations of state and is properly tested - [[PR439]](https://github.com/lanl/singularity-eos/pull/439) Add mean atomic mass and number to EOS API - [[PR437]](https://github.com/lanl/singularity-eos/pull/437) Fix segfault on HIP, clean up warnings, add strict sanitizer test diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index e9bed5ddfb8..ae1db924a5a 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1143,18 +1143,70 @@ quantities as outputs. Methods Used for Mixed Cell Closures -------------------------------------- -Several methods were developed in support of mixed cell closures. In particular: +Several methods were developed in support of mixed cell closures. In particular the function .. cpp:function:: Real MinimumDensity() const; -and +the function .. cpp:function:: Real MinimumTemperature() const; +and the function + +.. cpp:function:: Real MaximumDensity() const; + provide bounds for valid inputs into a table, which can be used by a -root finder to meaningful bound the root search. Similarly, +root finder to meaningful bound the root search. + +.. warning:: + + For unbounded equations of state, ``MinimumDensity`` and + ``MinimumTemperature`` will return zero, while ``MaximumDensity`` + will return a very large finite number. Which number you get, + however, is not guaranteed. You may wish to apply more sensible + bounds in your own code. + +Similarly, .. cpp:function:: Real RhoPmin(const Real temp) const; returns the density at which pressure is minimized for a given -temperature. This is again useful for root finds. +temperature. The function + +.. cpp:function:: Real MinimumPressure() const; + +provides the minimum pressure an equation of state supports, which may +be the most negative tension state. The function + +.. cpp:function:: Real MaximumPressureAtTemperature(const Real temp) const; + +provides a maximum possible pressure an equation of state supports at +a given temperature. (Most models are unbounded in pressure.) This is +again useful for root finds. + +The function + +.. code-block:: cpp + + template + void DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Indexer_t &&lambda, Real &rho, Real &sie) const; + +is designed for working in Pressure-Temperature space. Given a +pressure ``press`` and temperature ``temp``, it sets a density ``rho`` +and specific internal energy ``sie``. The ``lambda`` is optional and +defaults to a ``nullptr``. + +Typically this operation requires a root find. You may pass in an +initial guess for the density ``rho`` in-place and most EOS models +will use it. + +.. warning:: + + Pressure is not necessarily monotone in density and it may be double + valued. Thus you are not guaranteed to find the correct root and the + value of your initial guess may determine correctness. The fact that + ``rho`` may be used as an initial guess means you **must** pass in + an initialized variable, even if it is zero-initialized. Do not pass + uninitialized memory into this function! + diff --git a/eospac-wrapper/eospac_wrapper.cpp b/eospac-wrapper/eospac_wrapper.cpp index 3a2fdb8ee84..24ec3c61c09 100644 --- a/eospac-wrapper/eospac_wrapper.cpp +++ b/eospac-wrapper/eospac_wrapper.cpp @@ -25,9 +25,9 @@ namespace EospacWrapper { void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) { - constexpr int NT = 2; + constexpr int NT = 3; EOS_INTEGER tableHandle[NT]; - EOS_INTEGER tableType[NT] = {EOS_Info, EOS_Ut_DT}; + EOS_INTEGER tableType[NT] = {EOS_Info, EOS_Ut_DT, EOS_Pt_DT}; EOS_INTEGER commentsHandle[1]; EOS_INTEGER commentsType[1] = {EOS_Comment}; @@ -48,7 +48,8 @@ void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) { } } - eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Info", "EOS_Ut_DT"}, eospacWarn); + eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Info", "EOS_Ut_DT", "EOS_Pt_DT"}, + eospacWarn); for (int i = 0; i < numInfoTables; i++) { eosSafeTableInfo(&(tableHandle[i]), NI[i], infoItems[i].data(), infoVals[i].data(), @@ -66,6 +67,8 @@ void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) { metadata.TMax = temperatureFromSesame(infoVals[1][3]); metadata.sieMin = sieFromSesame(infoVals[1][4]); metadata.sieMax = sieFromSesame(infoVals[1][5]); + metadata.PMin = pressureFromSesame(infoVals[2][4]); + metadata.PMax = pressureFromSesame(infoVals[2][5]); metadata.numRho = static_cast(infoVals[1][6]); metadata.numT = static_cast(infoVals[1][7]); metadata.rhoConversionFactor = infoVals[1][8]; diff --git a/eospac-wrapper/eospac_wrapper.hpp b/eospac-wrapper/eospac_wrapper.hpp index a9c285d6567..169c6b1e459 100644 --- a/eospac-wrapper/eospac_wrapper.hpp +++ b/eospac-wrapper/eospac_wrapper.hpp @@ -59,6 +59,7 @@ class SesameMetadata { double rhoMin, rhoMax; double TMin, TMax; double sieMin, sieMax; + double PMin, PMax; double rhoConversionFactor; double TConversionFactor; double sieConversionFactor; diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 5ca9ed928e1..7f685dd88e2 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -16,11 +16,14 @@ #define _SINGULARITY_EOS_EOS_EOS_BASE_ #include +#include #include #include #include #include +#include +#include #include namespace singularity { @@ -781,6 +784,23 @@ class EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return 0; } + // Report maximum value of density. Default is unbounded. + // JMM: Should we use actual infinity, the largest real, or just a + // big number? For comparisons, actual infinity is better. It also + // has the advantage of being projective with modifiers that modify + // the max. On the other hand, it's more fraught if someone tries to + // put it into a formula without guarding against it. + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { return 1e100; } + + // These are for the PT space PTE solver to bound the iterations in + // a safe range. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return 0; } + // Gruneisen EOS's often have a maximum density, which implies a maximum pressure. + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } + PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { return 0.0; } @@ -836,6 +856,58 @@ class EosBase { PORTABLE_ALWAYS_THROW_OR_ABORT(msg); } + // JMM: This method is often going to be overloaded for special cases. + template + PORTABLE_INLINE_FUNCTION void + DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Indexer_t &&lambda, Real &rho, Real &sie) const { + using RootFinding1D::findRoot; // more robust but slower. Better default. + using RootFinding1D::Status; + + // Pressure is not monotone in density at low densities, which can + // prevent convergence. We want to approach tension from above, + // not below. Choose close to, but above, normal density for a + // metal like copper. + constexpr Real DEFAULT_RHO_GUESS = 12; + + CRTP copy = *(static_cast(this)); + + // P(rho) not monotone. When relevant, bound rhopmin. + Real rhomin = std::max(copy.RhoPmin(temp), copy.MinimumDensity()); + Real rhomax = copy.MaximumDensity(); + PORTABLE_REQUIRE(rhomax > rhomin, "max bound > min bound"); + + auto PofRT = [&](const Real r) { + return copy.PressureFromDensityTemperature(r, temp, lambda); + }; + Real rhoguess = rho; // use input density + if ((rhoguess <= rhomin) || (rhoguess >= rhomax)) { // avoid edge effects + if ((rhomin < DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS < rhomax)) { + rhoguess = DEFAULT_RHO_GUESS; + } else { + rhoguess = 0.5 * (rhomin + rhomax); + } + } + auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), + robust::EPS(), rho); + // JMM: This needs to not fail and instead return something sane. + // If root find failed to converge, density will at least be + // within bounds. + if (status != Status::SUCCESS) { + PORTABLE_WARN("DensityEnergyFromPressureTemperature failed to find root\n"); + } + sie = copy.InternalEnergyFromDensityTemperature(rho, temp, lambda); + return; + } + PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, + const Real temp, + Real &rho, + Real &sie) const { + CRTP copy = *(static_cast(this)); + copy.DensityEnergyFromPressureTemperature(press, temp, static_cast(nullptr), + rho, sie); + } + // Serialization /* The methodology here is there are *three* size methods all EOS's provide: diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index 763cb536955..0e6ef06e0cd 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -1104,6 +1104,7 @@ class EOSPAC : public EosBase { } PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rho_min_; } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return temp_min_; } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return press_min_; } private: static constexpr const unsigned long _preferred_input = @@ -1128,6 +1129,7 @@ class EOSPAC : public EosBase { Real dvdt_ref_ = 1; Real rho_min_ = 0; Real temp_min_ = 0; + Real press_min_ = 0; // TODO(JMM): Is the fact that EOS_INTEGER isn't a size_t a // problem? Could it ever realistically overflow? EOS_INTEGER shared_size_, packed_size_; @@ -1199,6 +1201,8 @@ inline EOSPAC::EOSPAC(const int matid, bool invert_at_setup, Real insert_data, rho_ref_ = m.normalDensity; rho_min_ = m.rhoMin; temp_min_ = m.TMin; + press_min_ = m.PMin; + // use std::max to hydrogen, in case of bad table AZbar_.Abar = std::max(1.0, m.meanAtomicMass); AZbar_.Zbar = std::max(1.0, m.meanAtomicNumber); diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index a2892e9a2af..d343705b847 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -170,6 +170,17 @@ class Gruneisen : public EosBase { s1, _C0, _s1, _s2, _s3, _G0, _b, _rho0, _T0, _P0, _Cv, _rho_max); _AZbar.PrintParams(); } + + // Report minimum values of density and temperature + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { return _rho_max; } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return PressureFromDensityInternalEnergy(0, 0); } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureAtTemperature(const Real T) const { + return MaxStableDensityAtTemperature(T); + } + template PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp index b6ade5ac19e..512f41ff31d 100644 --- a/singularity-eos/eos/eos_helmholtz.hpp +++ b/singularity-eos/eos/eos_helmholtz.hpp @@ -264,6 +264,12 @@ class HelmElectrons { std::size_t numTemp() const { return NTEMP; } PORTABLE_FORCEINLINE_FUNCTION std::size_t numRho() const { return NRHO; } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumDensity() const { return rhoMin(); } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumTemperature() const { return TMin_; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { return rhoMax(); } private: inline void InitDataFile_(const std::string &filename); @@ -676,16 +682,38 @@ class Helmholtz : public EosBase { ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, Real &dpde, Real &dvdt, Indexer_t &&lambda = static_cast(nullptr)) const { - // JMM: I'm not sure what to put here or if it matters. Some - // reference state, maybe stellar denity, would be appropriate. - PORTABLE_ALWAYS_ABORT("Stub"); + // JMM: Conditions for an oxygen burning shell in a stellar + // core. Not sure if that's the best choice. + rho = 1e7; + temp = 1.5e9; + FillEos(rho, temp, sie, press, cv, bmod, + thermalqs::specific_internal_energy | thermalqs::pressure | + thermalqs::specific_heat | thermalqs::bulk_modulus, + lambda); + dpde = GruneisenParamFromDensityTemperature(rho, temp, lambda) * rho; + dvdt = 0; } template PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - // This is only used for mixed cell closures. Stubbing it out for now. - PORTABLE_ALWAYS_ABORT("Stub"); + using RootFinding1D::regula_falsi; + using RootFinding1D::Status; + PORTABLE_REQUIRE(temp > = 0, "Non-negative temperature required"); + auto PofRT = [&](const Real r) { + return PressureFromDensityTemperature(r, temp, lambda); + }; + auto status = regula_falsi(PofRT, press, 1e7, electrons_.rhoMin(), + electrons_.rhoMax(), 1.0e-8, 1.0e-8, rho); + if (status != Status::SUCCESS) { + PORTABLE_THROW_OR_ABORT("Helmholtz::DensityEnergyFromPressureTemperature: " + "Root find failed to find a solution given P, T\n"); + } + if (rho < 0.) { + PORTABLE_THROW_OR_ABORT("Helmholtz::DensityEnergyFromPressureTemperature: " + "Root find produced negative energy solution\n"); + } + sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); } static std::string EosType() { return std::string("Helmholtz"); } static std::string EosPyType() { return EosType(); } diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index c315c1619aa..b85cd9175fd 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -125,6 +125,27 @@ class MGUsup : public EosBase { FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda = static_cast(nullptr)) const; + + // Since reference isotherm scales as eta = 1 -_rho0/rho, EOS + // diverges for rho << rho0. eta^2 is used, so... + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumDensity() const { return _RHOMINFAC * _rho0; } + + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { + if (_s > 1) { + return 0.99 * robust::ratio(_s * _rho0, _s - 1); + } else { // for s <= 1, no maximum, but we need to pick something. + return 1e3 * _rho0; + } // note that s < 0 implies unphysical shock derivative. + } + // Hugoniot pressure ill defined at reference density. On one side, + // negative. On the other positive. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return -1e100; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } + template PORTABLE_INLINE_FUNCTION void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, @@ -148,11 +169,6 @@ class MGUsup : public EosBase { _AZbar.PrintParams(); printf("\n\n"); } - // Density/Energy from P/T not unique, if used will give error - template - PORTABLE_INLINE_FUNCTION void - DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Indexer_t &&lambda, Real &rho, Real &sie) const; inline void Finalize() {} static std::string EosType() { return std::string("MGUsup"); } static std::string EosPyType() { return EosType(); } @@ -160,6 +176,7 @@ class MGUsup : public EosBase { private: static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; + Real _RHOMINFAC = std::sqrt(robust::EPS()); Real _rho0, _T0, _Cs, _s, _G0, _Cv0, _E0, _S0; MeanAtomicProperties _AZbar; }; @@ -200,6 +217,8 @@ PORTABLE_INLINE_FUNCTION Real MGUsup::HugTemperatureFromDensity(Real rho) const Real eta = 1.0 - robust::ratio(_rho0, rho); Real f1 = 1.0 - _s * eta; if (f1 <= 0.0) { + printf("f1, eta, rho, rho0, s = %.14e %.14e %.14e %.14e %.14e\n", f1, eta, rho, _rho0, + _s); PORTABLE_ALWAYS_THROW_OR_ABORT("MGUsup model parameters s and rho0 together with rho " "give a negative argument for a logarithm."); } @@ -350,13 +369,6 @@ PORTABLE_INLINE_FUNCTION Real MGUsup::BulkModulusFromDensityInternalEnergy( value = robust::ratio(_rho0, rho) * value; return value; } -// AEM: Give error since function is not well defined -template -PORTABLE_INLINE_FUNCTION void MGUsup::DensityEnergyFromPressureTemperature( - const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - EOS_ERROR("MGUsup::DensityEnergyFromPressureTemperature: " - "Not implemented.\n"); -} // AEM: We should add entropy and Gruneissen parameters here so that it is complete // If we add also alpha and BT, those should also be in here. template diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp index 7e400b89a7c..76da7607d07 100644 --- a/singularity-eos/eos/eos_powermg.hpp +++ b/singularity-eos/eos/eos_powermg.hpp @@ -157,11 +157,27 @@ class PowerMG : public EosBase { } printf("\n\n"); } - // Density/Energy from P/T not unique, if used will give error - template - PORTABLE_INLINE_FUNCTION void - DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Indexer_t &&lambda, Real &rho, Real &sie) const; + + // In principle, PowerMG should be valid for all densities. In + // practice, however, since the EOS is parametrized in terms of the + // compression, eta = 1 - rho0/rho, things will fail when rho is + // much less than rho0. Worse, since this is a power law in + // compression, it potentially depends on eta^M where M is the + // maximum index in the power series. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumDensity() const { + // JMM: I think formally this should be the mth root of machine + // epsilon times rho0, but for m = 20 or something that's of order + // 1. Things seem reasonably well behaved with this bound. They do + // NOT seem well behaved for, e.g., 10*rho0*machine epsilon. + return 1e-4 * _rho0; + } + // Essentially unbounded... I think. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return -1e100; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } + inline void Finalize() {} static std::string EosType() { return std::string("PowerMG"); } static std::string EosPyType() { return EosType(); } @@ -431,13 +447,6 @@ PORTABLE_INLINE_FUNCTION Real PowerMG::EntropyFromDensityInternalEnergy( } return value; } -// AEM: Give error since function is not well defined -template -PORTABLE_INLINE_FUNCTION void PowerMG::DensityEnergyFromPressureTemperature( - const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - EOS_ERROR("PowerMG::DensityEnergyFromPressureTemperature: " - "Not implemented.\n"); -} // AEM: We should add entropy and Gruneissen parameters here so that it is complete // If we add also alpha and BT, those should also be in here. template diff --git a/singularity-eos/eos/eos_sap_polynomial.hpp b/singularity-eos/eos/eos_sap_polynomial.hpp index 7d361190c9c..c7358c48fa5 100644 --- a/singularity-eos/eos/eos_sap_polynomial.hpp +++ b/singularity-eos/eos/eos_sap_polynomial.hpp @@ -162,6 +162,12 @@ class SAP_Polynomial : public EosBase { return rho / _rho0 - 1; } + // Essentially unbounded... I think. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return -1e100; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } + template PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, @@ -209,14 +215,6 @@ class SAP_Polynomial : public EosBase { printf(" b3 = %g\n", _b3); _AZbar.PrintParams(); } - template - PORTABLE_INLINE_FUNCTION void - DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Indexer_t &&lambda, Real &rho, Real &sie) const { - PORTABLE_WARN("This function is a stub for an incomplete EoS."); - sie = 0.0; - rho = 0.0; - } inline void Finalize() {} static std::string EosType() { return std::string("SAP_Polynomial"); } static std::string EosPyType() { return EosType(); } diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 7a9bb2e2853..6af2e644d1b 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -204,6 +204,10 @@ class SpinerEOSDependsRhoT : public EosBase { } PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return T_(lTMin_); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return PMin_; } + PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return _n_lambda; } PORTABLE_INLINE_FUNCTION @@ -305,6 +309,7 @@ class SpinerEOSDependsRhoT : public EosBase { Real lRhoMin_, lRhoMax_, rhoMax_; Real lRhoMinSearch_; Real lTMin_, lTMax_, TMax_; + Real PMin_; Real rhoNormal_, TNormal_, sieNormal_, PNormal_; Real CvNormal_, bModNormal_, dPdENormal_, dVdTNormal_; Real lRhoOffset_, lTOffset_; // offsets must be non-negative @@ -472,6 +477,13 @@ class SpinerEOSDependsRhoSie : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return TMin(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return PMin_; } + PORTABLE_INLINE_FUNCTION + Real RhoPmin(const Real temp) const { + return rho_at_pmin_.interpToReal(toLog_(temp, lTOffset_)); + } PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return _n_lambda; } @@ -515,22 +527,24 @@ class SpinerEOSDependsRhoSie : public EosBase { DataBox sie_; // depends on (rho,T) DataBox T_; // depends on (rho, sie) + DataBox rho_at_pmin_; SP5Tables dependsRhoT_; SP5Tables dependsRhoSie_; - int numRho_; + int numRho_, numT_; Real rhoNormal_, TNormal_, sieNormal_, PNormal_; Real CvNormal_, bModNormal_, dPdENormal_, dVdTNormal_; Real lRhoMin_, lRhoMax_, rhoMax_; DataBox PlRhoMax_, dPdRhoMax_; - + Real PMin_; Real lRhoOffset_, lTOffset_, lEOffset_; // offsets must be non-negative #define DBLIST \ - &sie_, &T_, &(dependsRhoT_.P), &(dependsRhoT_.bMod), &(dependsRhoT_.dPdRho), \ - &(dependsRhoT_.dPdE), &(dependsRhoT_.dTdRho), &(dependsRhoT_.dTdE), \ - &(dependsRhoT_.dEdRho), &(dependsRhoSie_.P), &(dependsRhoSie_.bMod), \ - &(dependsRhoSie_.dPdRho), &(dependsRhoSie_.dPdE), &(dependsRhoSie_.dTdRho), \ - &(dependsRhoSie_.dTdE), &(dependsRhoSie_.dEdRho), &PlRhoMax_, &dPdRhoMax_ + &sie_, &T_, &rho_at_pmin_, &(dependsRhoT_.P), &(dependsRhoT_.bMod), \ + &(dependsRhoT_.dPdRho), &(dependsRhoT_.dPdE), &(dependsRhoT_.dTdRho), \ + &(dependsRhoT_.dTdE), &(dependsRhoT_.dEdRho), &(dependsRhoSie_.P), \ + &(dependsRhoSie_.bMod), &(dependsRhoSie_.dPdRho), &(dependsRhoSie_.dPdE), \ + &(dependsRhoSie_.dTdRho), &(dependsRhoSie_.dTdE), &(dependsRhoSie_.dEdRho), \ + &PlRhoMax_, &dPdRhoMax_ std::vector GetDataBoxPointers_() const { return std::vector{DBLIST}; } @@ -771,11 +785,11 @@ inline herr_t SpinerEOSDependsRhoT::loadDataboxes_(const std::string &matid_str, rho_at_pmin_.resize(numT_); rho_at_pmin_.setRange(0, P_.range(0)); for (int i = 0; i < numT_; i++) { - Real pmin = std::numeric_limits::max(); + PMin_ = std::numeric_limits::max(); int jmax = -1; for (int j = 0; j < numRho_; j++) { - if (P_(j, i) < pmin) { - pmin = P_(j, i); + if (P_(j, i) < PMin_) { + PMin_ = P_(j, i); jmax = j; } } @@ -1595,6 +1609,7 @@ herr_t SpinerEOSDependsRhoSie::loadDataboxes_(const std::string &matid_str, hid_ // Metadata for root finding extrapolation numRho_ = sie_.dim(2); + numT_ = sie_.dim(1); lRhoMin_ = sie_.range(1).min(); lRhoMax_ = sie_.range(1).max(); rhoMax_ = fromLog_(lRhoMax_, lRhoOffset_); @@ -1603,6 +1618,22 @@ herr_t SpinerEOSDependsRhoSie::loadDataboxes_(const std::string &matid_str, hid_ PlRhoMax_ = dependsRhoT_.P.slice(numRho_ - 1); dPdRhoMax_ = dependsRhoT_.dPdRho.slice(numRho_ - 1); + // fill in minimum pressure as a function of temperature + rho_at_pmin_.resize(numT_); + rho_at_pmin_.setRange(0, sie_.range(0)); + for (int i = 0; i < numT_; i++) { + PMin_ = std::numeric_limits::max(); + int jmax = -1; + for (int j = 0; j < numRho_; j++) { + if (dependsRhoT_.P(j, i) < PMin_) { + PMin_ = dependsRhoT_.P(j, i); + jmax = j; + } + } + if (jmax < 0) printf("Failed to find minimum pressure.\n"); + rho_at_pmin_(i) = fromLog_(dependsRhoT_.P.range(1).x(jmax), lRhoOffset_); + } + // reference state Real lRhoNormal = toLog_(rhoNormal_, lRhoOffset_); // if rho normal not on the table, set it to the middle diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 4ff7788dda2..d7c5ecb3190 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -237,6 +237,7 @@ class StellarCollapse : public EosBase { } PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return TMin(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return _n_lambda; } inline RootFinding1D::Status rootStatus() const { return status_; } @@ -653,7 +654,26 @@ PORTABLE_INLINE_FUNCTION Real StellarCollapse::GruneisenParamFromDensityInternal template PORTABLE_INLINE_FUNCTION void StellarCollapse::DensityEnergyFromPressureTemperature( const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - EOS_ERROR("StellarCollapse::DensityEnergyFromPRessureTemperature is a stub"); + using RootFinding1D::regula_falsi; + using RootFinding1D::Status; + checkLambda_(lambda); + Real lrguess = lRho_(rho); + Real lT = lT_(temp); + Real lP = P2lP_(press); + Real Ye = lambda[Lambda::Ye]; + + if ((lrguess < lRhoMin_) || (lrguess > lRhoMax_)) { + lrguess = lRho_(rhoNormal_); + } + // Since EOS is tabulated in logrho, use that for root find. + auto lPofRT = [&](Real lR) { return lP_.interpToReal(Ye, lT, lR); }; + auto status = regula_falsi(lPofRT, lP, lrguess, lRhoMin_, lRhoMax_, ROOT_THRESH, + ROOT_THRESH, lrguess); + + Real lE = lE_.interpToReal(Ye, lT, lrguess); + rho = rho_(lrguess); + sie = le2e_(lE); + lambda[Lambda::lT] = lT; } template diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index 31dcb7867a3..7b48a2e7e04 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -317,6 +317,16 @@ class Variant { }, eos_); } + PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, + const Real temp, + Real &rho, + Real &sie) const { + return mpark::visit( + [&press, &temp, &rho, &sie](const auto &eos) { + return eos.DensityEnergyFromPressureTemperature(press, temp, rho, sie); + }, + eos_); + } PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { @@ -333,6 +343,23 @@ class Variant { return mpark::visit([](const auto &eos) { return eos.MinimumTemperature(); }, eos_); } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { + return mpark::visit([](const auto &eos) { return eos.MaximumDensity(); }, eos_); + } + + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { + return mpark::visit([](const auto &eos) { return eos.MinimumPressure(); }, eos_); + } + + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureAtTemperature(const Real temp) const { + return mpark::visit( + [&temp](const auto &eos) { return eos.MaximumPressureAtTemperature(temp); }, + eos_); + } + // Atomic mass/atomic number functions PORTABLE_INLINE_FUNCTION Real MeanAtomicMass() const { diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp index 8ea913e4049..67d3e722159 100644 --- a/singularity-eos/eos/eos_vinet.hpp +++ b/singularity-eos/eos/eos_vinet.hpp @@ -131,6 +131,18 @@ class Vinet : public EosBase { ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, Real &dpde, Real &dvdt, Indexer_t &&lambda = static_cast(nullptr)) const; + + // EOS depends on (rho0/rho)^(1/3) so this is a reasonable bound + // before the EOS will become ill-behaved. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumDensity() const { return std::cbrt(10 * robust::EPS()) * _rho0; } + + // Essentially unbounded... I think. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return -1e100; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } + // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here SG_ADD_BASE_CLASS_USINGS(Vinet) @@ -152,11 +164,6 @@ class Vinet : public EosBase { } printf("\n\n"); } - // Density/Energy from P/T not unique, if used will give error - template - PORTABLE_INLINE_FUNCTION void - DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Indexer_t &&lambda, Real &rho, Real &sie) const; inline void Finalize() {} static std::string EosType() { return std::string("Vinet"); } static std::string EosPyType() { return EosType(); } @@ -375,13 +382,6 @@ PORTABLE_INLINE_FUNCTION Real Vinet::BulkModulusFromDensityInternalEnergy( Vinet_F_DT_func(rho, temp, output); return output[7] * output[7] * rho; } -// AEM: Give error since function is not well defined -template -PORTABLE_INLINE_FUNCTION void Vinet::DensityEnergyFromPressureTemperature( - const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - EOS_ERROR("Vinet::DensityEnergyFromPressureTemperature: " - "Not implemented.\n"); -} // AEM: We should add entropy and Gruneissen parameters here so that it is complete // If we add also alpha and BT, those should also be in here. template diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 2e7519f7509..853f1b2f944 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -248,6 +248,15 @@ class UnitSystem : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return inv_temp_unit_ * t_.MinimumTemperature(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return inv_rho_unit_ * t_.MaximumDensity(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return inv_press_unit_ * t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return inv_press_unit_ * t_.MaximumPressureAtTemperature(temp_unit_ * temp); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index fb072f7d7d6..67a5f5185aa 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -243,6 +243,22 @@ class BilinearRampEOS : public EosBase> { return; } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { + return t_.MinimumDensity(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { + return t_.MinimumTemperature(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return t_.MaximumDensity(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return t_.MaximumPressureAtTemperature(temp); + } + template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( const Real rho, const Real temperature, diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index f497e359246..d6cbe4d76c7 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -155,6 +155,15 @@ class RelativisticEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return t_.MinimumTemperature(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return t_.MaximumDensity(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return t_.MaximumPressureAtTemperature(temp); + } static constexpr unsigned long PreferredInput() { return T::PreferredInput(); } diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index 60d522026fe..d48e6d9504a 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -341,6 +341,15 @@ class ScaledEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return t_.MinimumTemperature(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return inv_scale_ * t_.MaximumDensity(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return t_.MaximumPressureAtTemperature(temp); + } PORTABLE_INLINE_FUNCTION Real MeanAtomicMass() const { return inv_scale_ * t_.MeanAtomicMass(); } diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index f42003a4121..933f655366b 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -351,6 +351,15 @@ class ShiftedEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return t_.MinimumTemperature(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return t_.MaximumDensity(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return t_.MaximumPressureAtTemperature(temp); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp index eb9d68313cf..9b3505a3b71 100644 --- a/singularity-eos/eos/modifiers/zsplit_eos.hpp +++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp @@ -212,6 +212,16 @@ class ZSplit : public EosBase> { // dvdt, dpde unchanged } + template + PORTABLE_FUNCTION void + DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Indexer_t &&lambda, Real &rho, Real &sie) const { + const Real scale = GetScale_(lambda); + const Real iscale = GetInvScale_(lambda); + t_.DensityEnergyFromPressureTemperature(iscale * press, temp, lambda, rho, sie); + sie *= scale; + } + PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return NL + t_.nlambda(); } static constexpr unsigned long PreferredInput() { return T::PreferredInput(); } @@ -229,6 +239,15 @@ class ZSplit : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return t_.MinimumTemperature(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return t_.MaximumDensity(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return t_.MaximumPressureAtTemperature(temp); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index b5bf7bde90f..3b73e574815 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -34,6 +34,7 @@ #include #include +#include inline std::string demangle(const char *name) { @@ -143,6 +144,38 @@ inline void compare_two_eoss(const E1 &test_e, const E2 &ref_e) { return; } +template +PORTABLE_INLINE_FUNCTION bool +CheckRhoSieFromPT(EOS eos, Real rho, Real T, + Indexer_t &&lambda = static_cast(nullptr)) { + const Real P = eos.PressureFromDensityTemperature(rho, T, lambda); + const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T, lambda); + Real rtest = 12; // set these to something + Real etest = 1; + eos.DensityEnergyFromPressureTemperature(P, T, lambda, rtest, etest); + Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda); + Real residual = P_test - P; + Real frac_residual = + std::min(std::abs(singularity::robust::ratio(residual, P)), std::abs(residual)); + bool results_good = (isClose(rho, rtest, 1e-8) && isClose(sie, etest, 1e-8)) + // This is as good as it will get sometimes. + || (std::abs(frac_residual) <= 1e-8); + if (!results_good) { + printf("RhoSie of PT failure!\n" + "\trho_true = %.14e\n" + "\tsie_true = %.14e\n" + "\tP = %.14e\n" + "\tT = %.14e\n" + "\trho = %.14e\n" + "\tsie = %.14e\n" + "\tP_test = %.14e\n" + "\tresidual = %.14e\n" + "\tfracres = %.14e\n", + rho, sie, P, T, rtest, etest, P_test, residual, frac_residual); + } + return results_good; +} + // Macro that checks for an exception or is a no-op depending on // whether or not a non-serial backend is supplied #ifdef PORTABILITY_STRATEGY_NONE diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp index 77735704505..da6bc18771d 100644 --- a/test/test_eos_carnahan_starling.cpp +++ b/test/test_eos_carnahan_starling.cpp @@ -161,6 +161,18 @@ SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") { "Energy"); } } + + WHEN("We check RhoSie from PT") { + int nwrong = 0; + portableReduce( + "Check RhoSieFromPT", 0, 1, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += + !CheckRhoSieFromPT(eos, 7.8290736890381501e-03, 1.5320999999999999e+03); + }, + nwrong); + REQUIRE(nwrong == 0); + } } } } diff --git a/test/test_eos_helmholtz.cpp b/test/test_eos_helmholtz.cpp index 01a2aa11038..037835cacaf 100644 --- a/test/test_eos_helmholtz.cpp +++ b/test/test_eos_helmholtz.cpp @@ -125,13 +125,16 @@ SCENARIO("Helmholtz equation of state - Table interpolation (tgiven)", "[Helmhol eos.BulkModulusFromDensityTemperature(rho_in[i], temp_in[j], lambda); Real gruen = eos.GruneisenParamFromDensityTemperature(rho_in[i], temp_in[j], lambda); - if (!isClose(ein, ein_ref[k], 1e-10)) nwrong += 1; if (!isClose(press, press_ref[k], 1e-10)) nwrong += 1; if (!isClose(cv, cv_ref[k], 1e-6)) nwrong += 1; if (!isClose(bulkmod, bulkmod_ref[k], 1e-8)) nwrong += 1; if (!isClose(gruen, gruen_ref[k], 1e-6)) nwrong += 1; + // RhoSie of PT + nwrong += !CheckRhoSieFromPT(eos, rho_in[i], temp_in[i], lambda); + + // Deserialized EOS ein = eos_2.InternalEnergyFromDensityTemperature(rho_in[i], temp_in[j], lambda); press = eos_2.PressureFromDensityTemperature(rho_in[i], temp_in[j], lambda); diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp index b24a94bb3fa..58a35149298 100644 --- a/test/test_eos_ideal.cpp +++ b/test/test_eos_ideal.cpp @@ -106,10 +106,10 @@ SCENARIO("Ideal gas mean atomic properties", portableReduce( "Check mean atomic number", 0, N, PORTABLE_LAMBDA(const int i, int &nw) { - double rho = i; - double T = 100.0 * i; - double Ab_eval = device_eos.MeanAtomicMassFromDensityTemperature(rho, T); - double Zb_eval = device_eos.MeanAtomicNumberFromDensityTemperature(rho, T); + Real rho = i; + Real T = 100.0 * i; + Real Ab_eval = device_eos.MeanAtomicMassFromDensityTemperature(rho, T); + Real Zb_eval = device_eos.MeanAtomicNumberFromDensityTemperature(rho, T); nw += !(isClose(Ab_eval, Abar, 1e-12)) + !(isClose(Zb_eval, Zbar, 1e-12)); }, nwrong); @@ -120,6 +120,29 @@ SCENARIO("Ideal gas mean atomic properties", } } +SCENARIO("Ideal gas density energy from prssure temperature", + "[IdealGas][DensityEnergyFromPressureTemperature]") { + constexpr Real Cv = 5.0; + constexpr Real gm1 = 0.4; + GIVEN("An ideal gas") { + EOS host_eos = IdealGas(gm1, Cv); + auto device_eos = host_eos.GetOnDevice(); + WHEN("We compute density and energy from pressure and temperature") { + constexpr int N = 100; + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 1, N, + PORTABLE_LAMBDA(const int i, int &nw) { + Real rho = i; + Real T = 100.0 * i; + nw += !CheckRhoSieFromPT(device_eos, rho, T); + }, + nwrong); + THEN("There are no errors") { REQUIRE(nwrong == 0); } + } + } +} + // A non-standard pattern where we do a reduction class CheckPofRE { public: @@ -283,5 +306,18 @@ SCENARIO("Ideal electron gas", "[IdealGas][IdealEelctrons]") { REQUIRE(nwrong == 0); } } + + WHEN("We compute density and energy from pressure and temperature") { + constexpr int N = 100; + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 1, N, + PORTABLE_LAMBDA(const int i, int &nw) { + Real ll[1] = {static_cast(i)}; + nw += !CheckRhoSieFromPT(eos, rho, T, ll); + }, + nwrong); + THEN("There are no errors") { REQUIRE(nwrong == 0); } + } } } diff --git a/test/test_eos_mgusup.cpp b/test/test_eos_mgusup.cpp index df35715eac5..162709cc0cb 100644 --- a/test/test_eos_mgusup.cpp +++ b/test/test_eos_mgusup.cpp @@ -408,6 +408,17 @@ SCENARIO("MGUsup EOS rho T", "[VectorEOS][MGUsupEOS]") { } } + WHEN("a [rho,sie](P, T) lookup is performed") { + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += !CheckRhoSieFromPT(eos, v_density[i], v_temperature[i]); + }, + nwrong); + THEN("There are no errors") { REQUIRE(nwrong == 0); } + } + portableFor( "entropy", 0, num, PORTABLE_LAMBDA(const int i) { v_entropy[i] = diff --git a/test/test_eos_noble_abel.cpp b/test/test_eos_noble_abel.cpp index eba52ef49bb..c7f3d75b4b0 100644 --- a/test/test_eos_noble_abel.cpp +++ b/test/test_eos_noble_abel.cpp @@ -125,6 +125,18 @@ SCENARIO("NobleAbel1", "[NobleAbel][NobleAbel1]") { } } + WHEN("A [rho,e](P, T) loopup is performed") { + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += + !CheckRhoSieFromPT(eos, 7.2347087589269467e-03, 4.0000000000000000e+03); + }, + nwrong); + THEN("There are no errors") { REQUIRE(nwrong == 0); } + } + WHEN("A B_S(rho, e) lookup is performed") { eos.BulkModulusFromDensityInternalEnergy(v_density, v_energy, v_bulkmodulus, num); #ifdef PORTABILITY_STRATEGY_KOKKOS diff --git a/test/test_eos_powermg.cpp b/test/test_eos_powermg.cpp index d6565a11105..2c0e7e7eace 100644 --- a/test/test_eos_powermg.cpp +++ b/test/test_eos_powermg.cpp @@ -271,6 +271,17 @@ SCENARIO("PowerMG EOS rho sie", "[VectorEOS][PowerMGEOS]") { } } + WHEN("a [rho,sie](P, T) lookup is performed") { + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 0, 1, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += !CheckRhoSieFromPT(eos, 7.285, 7.78960970661031411e+02); + }, + nwrong); + THEN("There are no errors") { REQUIRE(nwrong == 0); } + } + WHEN("A S(rho, e) lookup is performed") { portableFor( "Test entropy", 0, num, PORTABLE_LAMBDA(const int i) { diff --git a/test/test_eos_stellar_collapse.cpp b/test/test_eos_stellar_collapse.cpp index 53524019bdf..7087e4398cb 100644 --- a/test/test_eos_stellar_collapse.cpp +++ b/test/test_eos_stellar_collapse.cpp @@ -278,6 +278,9 @@ SCENARIO("Stellar Collapse EOS", "[StellarCollapse]") { if (!isClose(b1, b2)) { nwrong_d() += 1; } + if (!CheckRhoSieFromPT(sc_d, R, T, lambda)) { + nwrong_d() += 1; + } }); #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::deep_copy(nwrong_h, nwrong_d); diff --git a/test/test_eos_stiff.cpp b/test/test_eos_stiff.cpp index a34be2a2949..3f741daf4b5 100644 --- a/test/test_eos_stiff.cpp +++ b/test/test_eos_stiff.cpp @@ -583,6 +583,16 @@ SCENARIO("Test Stiff Gas Entropy Calls", "[StiffGas][StiffGas5]") { array_compare(num, density, energy, h_entropy, entropy_true, "Density", "Energy"); } + THEN("[rho, e](P, T) also works") { + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += !CheckRhoSieFromPT(eos, v_density[i], v_local_temp[i]); + }, + nwrong); + AND_THEN("There are no errors") { REQUIRE(nwrong == 0); } + } } } } diff --git a/test/test_eos_vinet.cpp b/test/test_eos_vinet.cpp index 0c7f15d9960..916f41aabbc 100644 --- a/test/test_eos_vinet.cpp +++ b/test/test_eos_vinet.cpp @@ -346,6 +346,16 @@ SCENARIO("Vinet EOS rho T", "[VectorEOS][VinetEOS]") { array_compare(num, density, temperature, h_pressure, pressure_true, "Density", "Temperature"); } + THEN("[rho, e](P, T) also works") { + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += !CheckRhoSieFromPT(eos, v_density[i], v_temperature[i]); + }, + nwrong); + AND_THEN("There are no errors") { REQUIRE(nwrong == 0); } + } } portableFor( diff --git a/test/test_eos_zsplit.cpp b/test/test_eos_zsplit.cpp index aeaff8c18c6..f4eecd35f64 100644 --- a/test/test_eos_zsplit.cpp +++ b/test/test_eos_zsplit.cpp @@ -88,6 +88,11 @@ SCENARIO("ZSplit of Ideal Gas", "[ZSplit][IdealGas][IdealElectrons]") { printf("Bad sie Sum! %d %.14e %.14e %.14e\n", i, et, ei, ee); nw += 1; } + + if (Z > 0) { + nw += !CheckRhoSieFromPT(eos_zi, rho, temp, lambda); + nw += !CheckRhoSieFromPT(eos_ze, rho, temp, lambda); + } }, nwrong); THEN("Pe + Pi = Pt") { REQUIRE(nwrong == 0); } @@ -119,6 +124,9 @@ SCENARIO("ZSplit of Ideal Gas", "[ZSplit][IdealGas][IdealElectrons]") { printf("Bad sie ionized! %d %.14e %.14e\n", i, e_fi, e_ie); nw += 1; } + + nw += !CheckRhoSieFromPT(eos_ze, rho, T, lambda); + nw += !CheckRhoSieFromPT(eos_ie, rho, T, lambda); }, nwrong); THEN("Pe_fi == PE_ie") { REQUIRE(nwrong == 0); }