From 9758631f10819d37090868d41b4674f3128b34d3 Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Thu, 28 Aug 2025 12:42:04 -0600 Subject: [PATCH 01/25] First attempt at implementing an EOS haven't tested compiling yet. --- singularity-eos/eos/eos_simple_macaw.hpp | 262 +++++++++++++++++++++++ singularity-eos/eos/eos_type_lists.hpp | 2 +- test/test_eos_simple_macaw.cpp | 63 ++++++ 3 files changed, 326 insertions(+), 1 deletion(-) create mode 100644 singularity-eos/eos/eos_simple_macaw.hpp create mode 100644 test/test_eos_simple_macaw.cpp diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp new file mode 100644 index 00000000000..b47441117c7 --- /dev/null +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -0,0 +1,262 @@ +//------------------------------------------------------------------------------ +// © 2021-2025. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#ifndef _SINGULARITY_EOS_EOS_EOS_SIMPLE_MACAW_HPP_ +#define _SINGULARITY_EOS_EOS_EOS_SIMPLE_MACAW_HPP_ + +// stdlib +#include +#include +#include + +// Ports-of-call +#include + +// Base stuff +#include +#include +#include +#include + +namespace singularity { + +using namespace eos_base; + +class SimpleMACAW : public EosBase { + public: + SimpleMACAW() = default; + PORTABLE_INLINE_FUNCTION + SimpleMACAW(Real A, Real B, Real Cvinf, Real v0, Real T0, Real Gc, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) + : _A(A), _B(B), _Cvinf(Cvinf), _v0(v0), _T0(T0), _Gc(Gc), _AZbar(AZbar) { + CheckParams(); + } + SimpleMACAW GetOnDevice() { return *this; } + template + PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real x = std::pow( rho * _v0, _Gc); + const Real sie_term = (sie - SieColdCurve(1.0 / rho)) / _Cvinf; + // It's a quadratic equation to solve for temperature from sie and rho + const Real b = -sie_term; + const Real c = -x * sie_term; + return 0.5 * (-b + std::sqrt(b * b - 4.0 * c)); + } + PORTABLE_INLINE_FUNCTION void CheckParams() const { + PORTABLE_ALWAYS_REQUIRE(_A > 0, "Parameter, 'A', must be positive"); + PORTABLE_ALWAYS_REQUIRE(_B > 0, "Parameter, 'B', must be positive"); + PORTABLE_ALWAYS_REQUIRE(_v0 > 0, "Reference specific volume, 'v0', must be positive"); + PORTABLE_ALWAYS_REQUIRE(_T0 > 0, "Reference temperature, 'T0', must be positive"); + PORTABLE_ALWAYS_REQUIRE(_Cvinf > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); + PORTABLE_ALWAYS_REQUIRE(_Gc > 0 && _Gc < 1, "Gruneisen coefficient, 'Gc', must be in the interval (0,1)"); + _AZbar.CheckParams(); + } + template + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real v = 1.0 / rho; + return rho * (SieColdCurve(v) + SieThermalPortion(v, temperature)); + } + template + PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real v = 1.0 / rho; + return PressureColdCurve(v) + PressureThermalPortion(v, temperature); + } + template + PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + const v = 1.0 / rho; + return PressureColdCurve(v) + _Gc * rho * (sie - SieColdCurve(v)); + } + + // Entropy member functions + template + PORTABLE_INLINE_FUNCTION Real + EntropyFromDensityTemperature(const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real tau = robust::ratio(temperature, TemperatureScale(rho)); + return _Cvinf * (robust::ratio(tau, 1.0 + tau) + std::log(1.0 + tau)); + } + template + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real T = TemperatureFromDensityInternalEnergy(rho, sie); + return EntropyFromDensityTemperature(rho, T); + } + + // Cold curve, thermal portion, and other helper functions + template + PORTABLE_INLINE_FUNCTION Real SieColdCurve( + const Real v, Indexer_t &&lambda = static_cast(nullptr)) const { + const Real x = robust::ratio(v, _v0); + return _A * _v0 * (std::pow(x, -_B) + _B * x - (_B + 1.0)); + } + template + PORTABLE_INLINE_FUNCTION Real SieThermalPortion( + const Real v, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real x = robust::ratio(v, _v0); + return _Cvinf * T * (1.0 - _T0 / (_T0 + T * std::pow(x, _Gc))); + } + template + PORTABLE_INLINE_FUNCTION Real PressureColdCurve( + const Real v, Indexer_t &&lambda = static_cast(nullptr)) const { + const Real x = robust::ratio(v, _v0); + return _A * _B * (std::pow(x, -(_B + 1.0)) - 1.0); + } + template + PORTABLE_INLINE_FUNCTION Real PressureThermalPortion( + const Real v, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return rho * _Cvinf * _Gc * pow<2>(T) / (_T0 * std::pow(rho * _v0, _Gc) + T); + } + template + PORTABLE_INLINE_FUNCTION Real TemperatureScale( + const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { + return _T0 * std::pow(rho * _v0, _Gc); + } + + // Thermodynamic derivatives + template + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real tau = robust::ratio(temperature, TemperatureScale(rho)); + const Real numerator = pow<2>(tau) + 2.0 * tau; + const Real denominator = pow<2>(tau + 1.0); + return _Cvinf * robust::ratio(numerator / denominator); + } + template + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real T = TemperatureFromDensityInternalEnergy(rho, sie); + return SpecificHeatFromDensityTemperature(rho, T); + } + template + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + return std::max(robust::SMALL(), _gm1 * (_gm1 + 1.0) * rho * _Cv * temperature); + } + template + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0) + + _Gc * (_Gc + 1.0) * rho * (sie - SieColdCurve(1.0 / rho)); + return std::max(robust::SMALL(), _gm1 * (_gm1 + 1.0) * (rho * (sie - _qq) - _Pinf)); + } + template + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + return _Gc; + } + template + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + return _Gc; + } + template + PORTABLE_INLINE_FUNCTION void + FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, + const unsigned long output, + Indexer_t &&lambda = static_cast(nullptr)) const; + //template + //PORTABLE_INLINE_FUNCTION void + //ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, + // Real &bmod, Real &dpde, Real &dvdt, + // Indexer_t &&lambda = static_cast(nullptr)) const { + // // use STP: 1 atmosphere, room temperature + // rho = _rho0; + // temp = _T0; + // sie = _sie0; + // press = _P0; + // cv = _Cv; + // bmod = _bmod0; + // dpde = _dpde0; + //} + // 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(SimpleMACAW) + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) + + static constexpr unsigned long PreferredInput() { return _preferred_input; } + PORTABLE_INLINE_FUNCTION void PrintParams() const { + printf("Simple MACAW Parameters:\nA = %g\nB = %g\nCvinf = %g\nGc = " + "%g\n", + _A, _B, _Cvinf, _Gc); + _AZbar.PrintParams(); + } + template + PORTABLE_INLINE_FUNCTION void + DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Indexer_t &&lambda, Real &rho, Real &sie) const { + sie = std::max( + _qq, + robust::ratio(press + (_gm1 + 1.0) * _Pinf, press + _Pinf) * _Cv * temp + _qq); + rho = std::max(robust::SMALL(), robust::ratio(press + _Pinf, _gm1 * _Cv * temp)); + } + inline void Finalize() {} + static std::string EosType() { return std::string("SimpleMACAW"); } + static std::string EosPyType() { return EosType(); } + + private: + Real _A, _B, _Cvinf, _Gc; + // reference values + Real _v0, _T0; + MeanAtomicProperties _AZbar; + static constexpr const unsigned long _preferred_input = + thermalqs::density | thermalqs::specific_internal_energy; +}; + +template +PORTABLE_INLINE_FUNCTION void +StiffGas::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, + const unsigned long output, Indexer_t &&lambda) const { + if (output & thermalqs::density && output & thermalqs::specific_internal_energy) { + if (output & thermalqs::pressure || output & thermalqs::temperature) { + UNDEFINED_ERROR; + } + DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie); + } + if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) { + if (output & thermalqs::density || output & thermalqs::temperature) { + UNDEFINED_ERROR; + } + sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); + } + if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) { + sie = robust::ratio(press + (_gm1 + 1.0) * _Pinf, _gm1 * rho) + _qq; + } + if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::temperature) + temp = TemperatureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::bulk_modulus) + bmod = BulkModulusFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::specific_heat) + cv = SpecificHeatFromDensityInternalEnergy(rho, sie); +} + +} // namespace singularity + +#endif // _SINGULARITY_EOS_EOS_EOS_SIMPLE_MACAW_HPP_ diff --git a/singularity-eos/eos/eos_type_lists.hpp b/singularity-eos/eos/eos_type_lists.hpp index 30aaced18fa..848f80edda0 100644 --- a/singularity-eos/eos/eos_type_lists.hpp +++ b/singularity-eos/eos/eos_type_lists.hpp @@ -49,7 +49,7 @@ using singularity::variadic_utils::transform_variadic_list; // all eos's static constexpr const auto full_eos_list = - tl +#include +#include +#include +#include +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +#include +#include +#include + +using singularity::SimpleMACAW; +using EOS = singularity::Variant; + +PORTABLE_INLINE_FUNCTION Real QuadFormulaMinus(Real a, Real b, Real c) { + return (-b - std::sqrt(b * b - 4 * a * c)) / (2 * a); +} + +constexpr Real REAL_TOL = std::numeric_limits::epsilon() * 1e4; // 1e-12 for double + +// Only run exception checking when we aren't offloading to the GPUs +SCENARIO("Zero temperature on cold curve", "[SimpleMACAWEOS][Temperature]") { + GIVEN("Parameters for a Simple MACAW EOS") { + // Unit conversions + //constexpr Real cm = 1.; + //constexpr Real us = 1e-06; + //constexpr Real Mbcc_per_g = 1e12; + // Gruneisen parameters for copper + constexpr Real A = 7.3; + constexpr Real B = 3.9; + constexpr Real Cvinf = 0.000389; + constexpr Real v0 = 1. / 8.952; + constexpr Real T0 = 150.; + constexpr Real Gc = 0.5; + // Create the EOS + EOS host_eos = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); + EOS eos = host_eos.GetOnDevice(); + Real rho = 0.5; + for (int i = 0; i < 10; i++) { + rho += rho + i; // cylce through a variety of densities + THEN("The temperature evaluated on the cold curve should produce zero") { + REQUIRE(eos.TemperatureFromDensityInternalEnergy(rho, eos.SieColdCurve(1.0 / rho))); + } + } + } +} From 53659f85219ad703c6893ffd29bc2ed7d14a28b4 Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Thu, 28 Aug 2025 15:13:16 -0600 Subject: [PATCH 02/25] Got it compiling, need to test and add more tests. --- singularity-eos/eos/eos_models.hpp | 1 + singularity-eos/eos/eos_simple_macaw.hpp | 49 +++++++++++++++--------- singularity-eos/eos/singularity_eos.cpp | 20 ++++++++++ test/CMakeLists.txt | 1 + 4 files changed, 52 insertions(+), 19 deletions(-) diff --git a/singularity-eos/eos/eos_models.hpp b/singularity-eos/eos/eos_models.hpp index f0ce27612ca..f2142666cdc 100644 --- a/singularity-eos/eos/eos_models.hpp +++ b/singularity-eos/eos/eos_models.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index b47441117c7..95ad8b62ed2 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -74,16 +74,23 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real v = 1.0 / rho; - return PressureColdCurve(v) + PressureThermalPortion(v, temperature); + return PressureColdCurve(1.0 / rho) + PressureThermalPortion(rho, temperature); } template PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - const v = 1.0 / rho; + const Real v = 1.0 / rho; return PressureColdCurve(v) + _Gc * rho * (sie - SieColdCurve(v)); } + template + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityPressure( + const Real rho, const Real P, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real v = robust::ratio(1.0, rho); + return SieColdCurve(v) + robust::ratio(v * (P - PressureColdCurve(v)), _Gc); + } + // Entropy member functions template @@ -123,9 +130,11 @@ class SimpleMACAW : public EosBase { } template PORTABLE_INLINE_FUNCTION Real PressureThermalPortion( - const Real v, const Real T, + const Real rho, const Real T, Indexer_t &&lambda = static_cast(nullptr)) const { - return rho * _Cvinf * _Gc * pow<2>(T) / (_T0 * std::pow(rho * _v0, _Gc) + T); + const Real numerator = rho * _Cvinf * _Gc * math_utils::pow<2>(T); + const Real denominator = _T0 * std::pow(rho * _v0, _Gc) + T; + return robust::ratio(numerator, denominator); } template PORTABLE_INLINE_FUNCTION Real TemperatureScale( @@ -133,15 +142,15 @@ class SimpleMACAW : public EosBase { return _T0 * std::pow(rho * _v0, _Gc); } - // Thermodynamic derivatives + // Specific heat capacity at constant volume template PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { const Real tau = robust::ratio(temperature, TemperatureScale(rho)); - const Real numerator = pow<2>(tau) + 2.0 * tau; - const Real denominator = pow<2>(tau + 1.0); - return _Cvinf * robust::ratio(numerator / denominator); + const Real numerator = math_utils::pow<2>(tau) + 2.0 * tau; + const Real denominator = math_utils::pow<2>(tau + 1.0); + return _Cvinf * robust::ratio(numerator, denominator); } template PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( @@ -150,20 +159,23 @@ class SimpleMACAW : public EosBase { const Real T = TemperatureFromDensityInternalEnergy(rho, sie); return SpecificHeatFromDensityTemperature(rho, T); } + // Isentropic bulk modulus template PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return std::max(robust::SMALL(), _gm1 * (_gm1 + 1.0) * rho * _Cv * temperature); + const Real sie = InternalEnergyFromDensityTemperature(rho, temperature); + return BulkModulusFromDensityInternalEnergy(rho, sie); } template PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0) - + _Gc * (_Gc + 1.0) * rho * (sie - SieColdCurve(1.0 / rho)); - return std::max(robust::SMALL(), _gm1 * (_gm1 + 1.0) * (rho * (sie - _qq) - _Pinf)); + const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0); + const Real term2 = _Gc * (_Gc + 1.0) * rho * (sie - SieColdCurve(1.0 / rho)); + return term1 + term2; } + // Gruneisen parameter template PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( const Real rho, const Real temperature, @@ -211,10 +223,9 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - sie = std::max( - _qq, - robust::ratio(press + (_gm1 + 1.0) * _Pinf, press + _Pinf) * _Cv * temp + _qq); - rho = std::max(robust::SMALL(), robust::ratio(press + _Pinf, _gm1 * _Cv * temp)); + // Solving for rho requires a Newton method... + sie = 0.0; + rho = 0.0; } inline void Finalize() {} static std::string EosType() { return std::string("SimpleMACAW"); } @@ -231,7 +242,7 @@ class SimpleMACAW : public EosBase { template PORTABLE_INLINE_FUNCTION void -StiffGas::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, +SimpleMACAW::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda) const { if (output & thermalqs::density && output & thermalqs::specific_internal_energy) { if (output & thermalqs::pressure || output & thermalqs::temperature) { @@ -246,7 +257,7 @@ StiffGas::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); } if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) { - sie = robust::ratio(press + (_gm1 + 1.0) * _Pinf, _gm1 * rho) + _qq; + sie = InternalEnergyFromDensityPressure(rho, press); } if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); if (output & thermalqs::temperature) diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index 78e23aa917a..936260ec359 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -124,6 +124,26 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, return init_sg_JWL(matindex, eos, A, B, R1, R2, w, rho0, Cv, def_en, def_v); } +int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, + const double Cinf, const double v0, const double T0, const double Gc, + int const *const enabled, double *const vals) { + assert(matindex >= 0); + EOS eosi = SGAPPLYMODSIMPLE(SimpleMACAW(A, B, Cvinf, v0, T0, Gc)); + if (enabled[3] == 1) { + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); + } + EOS eos_ = SGAPPLYMOD(SimpleMACAW(A, B, Cvinf, v0, T0, Gc)); + eos[matindex] = eos_.GetOnDevice(); + return 0; +} + +int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, + const double Cvinf, const double v0, const double T0, const double Gc) { + return init_sg_JWL(matindex, eos, A, B, Cvinf, v0, T0, Gc, def_en, def_v); +} + + 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, int const *const enabled, diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 903983d8051..cf857633975 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -20,6 +20,7 @@ add_executable( test_eos_sap_polynomial.cpp test_eos_noble_abel.cpp test_eos_stiff.cpp + test_eos_simple_macaw.cpp test_eos_carnahan_starling.cpp test_eos_vinet.cpp test_eos_mgusup.cpp From edf6a779619a680b2eafe3dc4b3eb1688077434e Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Sun, 31 Aug 2025 11:43:38 -0600 Subject: [PATCH 03/25] Working on iterative method for computing density from pressure and temperature. --- singularity-eos/eos/eos_simple_macaw.hpp | 44 ++++++++++++++++-------- singularity-eos/eos/singularity_eos.cpp | 4 +-- singularity-eos/eos/singularity_eos.f90 | 36 +++++++++++++++++++ singularity-eos/eos/singularity_eos.hpp | 4 +++ test/test_eos_simple_macaw.cpp | 6 +--- 5 files changed, 73 insertions(+), 21 deletions(-) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 95ad8b62ed2..73132c7f385 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -193,20 +193,21 @@ class SimpleMACAW : 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; - //template - //PORTABLE_INLINE_FUNCTION void - //ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, - // Real &bmod, Real &dpde, Real &dvdt, - // Indexer_t &&lambda = static_cast(nullptr)) const { - // // use STP: 1 atmosphere, room temperature - // rho = _rho0; - // temp = _T0; - // sie = _sie0; - // press = _P0; - // cv = _Cv; - // bmod = _bmod0; - // dpde = _dpde0; - //} + template + PORTABLE_INLINE_FUNCTION void + ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, + Real &bmod, Real &dpde, Real &dvdt, + Indexer_t &&lambda = static_cast(nullptr)) const { + // use STP: 1 atmosphere, room temperature + rho = 1.0 / _v0; + temp = _T0; + // TODO: Need to use the Newton solver for the density and energy here... + sie = 0.0; + press = PressureFromDensityTemperature(rho, temp); + cv = SpecificHeatFromDensityInternalEnergy(rho, sie); + bmod = BulkdModulusFromDensityInternalEnergy(rho, sie); + dpde = _Gc * rho; + } // 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(SimpleMACAW) @@ -223,6 +224,21 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { + // Setup lambda function for rho + auto f = [](const Real x){ + const Real term1 = _A * _B * std::pow(_v0 * x, _B+1.0) - press - _A * _B; + const Real term2 = (_T0 * std::pow(_v0 * x, _Gc) + temp) / x; + const Real term3 = _Cvinf * _Gc * math_utils::pow<2>(temp); + return term1 * term2 + term3; + } + + RootCounts root_info; + + regula_falsi(f, 0.0 /*target*/, 1.0 /*guess*/, + 1.0e-10 /*left bracket*/, 1.0e+6 /*right bracket*/, + 1.0e-6 /*x? tol*/, 1.0e-6 /*y? tol*/, + rho, root_info, true); + // Solving for rho requires a Newton method... sie = 0.0; rho = 0.0; diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index 936260ec359..867ea6e1c60 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -125,7 +125,7 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, } int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, - const double Cinf, const double v0, const double T0, const double Gc, + const double Cvinf, const double v0, const double T0, const double Gc, int const *const enabled, double *const vals) { assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(SimpleMACAW(A, B, Cvinf, v0, T0, Gc)); @@ -140,7 +140,7 @@ int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const doub int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, const double Cvinf, const double v0, const double T0, const double Gc) { - return init_sg_JWL(matindex, eos, A, B, Cvinf, v0, T0, Gc, def_en, def_v); + return init_sg_simple_MACAW(matindex, eos, A, B, Cvinf, v0, T0, Gc, def_en, def_v); } diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index e13e3d93eb0..24007874d42 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -35,6 +35,7 @@ module singularity_eos init_sg_IdealGas_f,& init_sg_Gruneisen_f,& init_sg_JWL_f,& + init_sg_SimpleMACAW_f,& init_sg_DavisProducts_f,& init_sg_DavisReactants_f,& init_sg_NobleAbel_f,& @@ -112,6 +113,19 @@ end function init_sg_Gruneisen end function init_sg_JWL end interface + interface + integer(kind=c_int) function & + init_sg_simple_MACAW(matindex, eos, A, B, Cvinf, v0, T0, Gc, sg_mods_enabled, & + sg_mods_values) & + bind(C, name='init_sg_simple_macaw') + import + integer(c_int), value, intent(in) :: matindex + type(c_ptr), value, intent(in) :: eos + real(kind=c_double), value, intent(in) :: A, B, Cvinf, v0, T0, Gc + type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values + end function init_sg_simple_MACAW + end interface + interface integer(kind=c_int) function & init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, & @@ -503,6 +517,28 @@ integer function init_sg_JWL_f(matindex, eos, A, B, R1, R2, w, rho0, Cv, & c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) end function init_sg_JWL_f + integer function init_sg_SimpleMACAW_f(matindex, eos, A, B, Cvinf, v0, T0, Gc, & + 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, Cvinf, v0, T0, Gc + 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 + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + + err = init_sg_simple_macaw(matindex-1, eos%ptr, A, B, Cvinf, v0, T0, Gc, & + c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) + end function init_sg_JWL_f + + integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & Cv, sg_mods_enabled, & sg_mods_values) & diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index d374eb40dd1..42e822e1a83 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -34,6 +34,10 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, const double R1, const double R2, const double w, const double rho0, const double Cv, int const *const enabled, double *const vals); +int init_sg_simple_MACAW(const int matindex, EOS *eos, const double A, const double B, + const double Cvinf, const double v0, const double T0, const double Gc, + int const *const enabled, double *const vals); + int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const double s1, const double s2, const double s3, const double G0, const double b, const double rho0, const double T0, const double P0, diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index 681204a8f8b..f7bc9f0ecf0 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -29,10 +29,6 @@ using singularity::SimpleMACAW; using EOS = singularity::Variant; -PORTABLE_INLINE_FUNCTION Real QuadFormulaMinus(Real a, Real b, Real c) { - return (-b - std::sqrt(b * b - 4 * a * c)) / (2 * a); -} - constexpr Real REAL_TOL = std::numeric_limits::epsilon() * 1e4; // 1e-12 for double // Only run exception checking when we aren't offloading to the GPUs @@ -51,7 +47,7 @@ SCENARIO("Zero temperature on cold curve", "[SimpleMACAWEOS][Temperature]") { constexpr Real Gc = 0.5; // Create the EOS EOS host_eos = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); - EOS eos = host_eos.GetOnDevice(); + //EOS eos = host_eos.GetOnDevice(); Real rho = 0.5; for (int i = 0; i < 10; i++) { rho += rho + i; // cylce through a variety of densities From 95b759965f04ca1322f6c316ed6e35f83551cc4d Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Thu, 4 Sep 2025 22:06:23 -0600 Subject: [PATCH 04/25] Got the nonlinear solve working for density from pressure temp. A few more things to do. --- CMakeLists.txt | 2 +- singularity-eos/eos/eos_simple_macaw.hpp | 63 +++++++++++++++--------- singularity-eos/eos/singularity_eos.cpp | 2 +- singularity-eos/eos/singularity_eos.f90 | 10 ++-- singularity-eos/eos/singularity_eos.hpp | 5 +- test/test_eos_simple_macaw.cpp | 48 ++++++++++++------ 6 files changed, 83 insertions(+), 47 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 4a8168b91fb..e17e789dd16 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -596,7 +596,7 @@ target_compile_options( ) if (SINGULARITY_STRICT_WARNINGS) target_compile_options(singularity-eos_Interface INTERFACE - -Wall -Werror -Wno-unknown-pragmas) + -Wall -Werror -Wno-unknown-pragmas) endif() if(TARGET singularity-eos_Library) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 73132c7f385..135b214ddea 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -47,12 +47,10 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real x = std::pow( rho * _v0, _Gc); - const Real sie_term = (sie - SieColdCurve(1.0 / rho)) / _Cvinf; - // It's a quadratic equation to solve for temperature from sie and rho - const Real b = -sie_term; - const Real c = -x * sie_term; - return 0.5 * (-b + std::sqrt(b * b - 4.0 * c)); + /* Equation (18) */ + const Real Delta_e = sie - SieColdCurve(1.0 / rho); + const Real discriminant = Delta_e * (Delta_e + 4.0 * _Cvinf * _T0 * std::pow(rho * _v0, _Gc)); + return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * _Cvinf); } PORTABLE_INLINE_FUNCTION void CheckParams() const { PORTABLE_ALWAYS_REQUIRE(_A > 0, "Parameter, 'A', must be positive"); @@ -91,7 +89,6 @@ class SimpleMACAW : public EosBase { return SieColdCurve(v) + robust::ratio(v * (P - PressureColdCurve(v)), _Gc); } - // Entropy member functions template PORTABLE_INLINE_FUNCTION Real @@ -125,8 +122,8 @@ class SimpleMACAW : public EosBase { template PORTABLE_INLINE_FUNCTION Real PressureColdCurve( const Real v, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real x = robust::ratio(v, _v0); - return _A * _B * (std::pow(x, -(_B + 1.0)) - 1.0); + const Real ratio = robust::ratio(v, _v0); + return _A * _B * (std::pow(ratio, -(_B + 1.0)) - 1.0); } template PORTABLE_INLINE_FUNCTION Real PressureThermalPortion( @@ -175,6 +172,19 @@ class SimpleMACAW : public EosBase { const Real term2 = _Gc * (_Gc + 1.0) * rho * (sie - SieColdCurve(1.0 / rho)); return term1 + term2; } + // Isothermal bulk modulus + template + PORTABLE_INLINE_FUNCTION Real IsothermalBulkModulusFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + /* As mentioned in the paper, from the constraints on the parameters, + * the isothermal bulk modulus is guaranteed to be positive. */ + const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0); + const Real numerator = rho * _T0 * (1.0 - _Gc) * std::pow(rho * _v0, 2.0 * _Gc) + temperature; + const Real denominator = _T0 * std::pow(rho * _v0, _Gc) + temperature; + return term1 + _Cvinf * _Gc * temperature * temperature * numerator / (denominator * denominator); + } + // Gruneisen parameter template PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( @@ -199,13 +209,13 @@ class SimpleMACAW : public EosBase { Real &bmod, Real &dpde, Real &dvdt, Indexer_t &&lambda = static_cast(nullptr)) const { // use STP: 1 atmosphere, room temperature + // TODO: Fix this subroutine... rho = 1.0 / _v0; temp = _T0; - // TODO: Need to use the Newton solver for the density and energy here... sie = 0.0; press = PressureFromDensityTemperature(rho, temp); cv = SpecificHeatFromDensityInternalEnergy(rho, sie); - bmod = BulkdModulusFromDensityInternalEnergy(rho, sie); + bmod = BulkModulusFromDensityInternalEnergy(rho, sie); dpde = _Gc * rho; } // Generic functions provided by the base class. These contain e.g. the vector @@ -224,24 +234,29 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { + /* Since the isothermal bulk modulus is always positive (assuming parameter constraints), + * this guarantees the function to solve is monotonic hence we always have a unique solution. */ + // Setup lambda function for rho - auto f = [](const Real x){ - const Real term1 = _A * _B * std::pow(_v0 * x, _B+1.0) - press - _A * _B; - const Real term2 = (_T0 * std::pow(_v0 * x, _Gc) + temp) / x; - const Real term3 = _Cvinf * _Gc * math_utils::pow<2>(temp); + auto f = [=](const Real x /* density */){ + const Real term1 = _A * _B * (std::pow(_v0 * x, _B + 1.0) - 1.0) - press; + const Real term2 = _T0 * std::pow(_v0 * x, _Gc) + temp; + const Real term3 = _Cvinf * _Gc * temp * temp * x; return term1 * term2 + term3; - } + }; - RootCounts root_info; + const RootFinding1D::RootCounts root_info; - regula_falsi(f, 0.0 /*target*/, 1.0 /*guess*/, - 1.0e-10 /*left bracket*/, 1.0e+6 /*right bracket*/, - 1.0e-6 /*x? tol*/, 1.0e-6 /*y? tol*/, - rho, root_info, true); + // Finding the density on pressure cold curve is a guaranteed upper bound + const Real rho_high = std::pow(press / (_A * _B) + 1.0, 1.0 / (_B + 1.0)) / _v0; + const Real rho_low = 0.0; // zero density is valid for `f` defined above. - // Solving for rho requires a Newton method... - sie = 0.0; - rho = 0.0; + regula_falsi(f, 0.0 /*target*/, 1.0 / _v0 /*guess*/, + rho_low /*left bracket*/, rho_high /*right bracket*/, + 1.0e-9 /*x? tol*/, 1.0e-9 /*y? tol*/, + rho, &root_info, true); + + sie = InternalEnergyFromDensityTemperature(rho, temp); } inline void Finalize() {} static std::string EosType() { return std::string("SimpleMACAW"); } diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index 867ea6e1c60..c4d5418a4ee 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -140,7 +140,7 @@ int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const doub int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, const double Cvinf, const double v0, const double T0, const double Gc) { - return init_sg_simple_MACAW(matindex, eos, A, B, Cvinf, v0, T0, Gc, def_en, def_v); + return init_sg_SimpleMACAW(matindex, eos, A, B, Cvinf, v0, T0, Gc, def_en, def_v); } diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index 24007874d42..71a556206b2 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -115,15 +115,15 @@ end function init_sg_JWL interface integer(kind=c_int) function & - init_sg_simple_MACAW(matindex, eos, A, B, Cvinf, v0, T0, Gc, sg_mods_enabled, & + init_sg_SimpleMACAW(matindex, eos, A, B, Cvinf, v0, T0, Gc, sg_mods_enabled, & sg_mods_values) & - bind(C, name='init_sg_simple_macaw') + bind(C, name='init_sg_SimpleMACAW') import integer(c_int), value, intent(in) :: matindex type(c_ptr), value, intent(in) :: eos real(kind=c_double), value, intent(in) :: A, B, Cvinf, v0, T0, Gc type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values - end function init_sg_simple_MACAW + end function init_sg_SimpleMACAW end interface interface @@ -534,9 +534,9 @@ integer function init_sg_SimpleMACAW_f(matindex, eos, A, B, Cvinf, v0, T0, Gc, & if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values - err = init_sg_simple_macaw(matindex-1, eos%ptr, A, B, Cvinf, v0, T0, Gc, & + err = init_sg_SimpleMACAW(matindex-1, eos%ptr, A, B, Cvinf, v0, T0, Gc, & c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) - end function init_sg_JWL_f + end function init_sg_SimpleMACAW_f integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index 42e822e1a83..d196f01564e 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -34,7 +34,7 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, const double R1, const double R2, const double w, const double rho0, const double Cv, int const *const enabled, double *const vals); -int init_sg_simple_MACAW(const int matindex, EOS *eos, const double A, const double B, +int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, const double Cvinf, const double v0, const double T0, const double Gc, int const *const enabled, double *const vals); @@ -148,6 +148,9 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, const double R1, const double R2, const double w, const double rho0, const double Cv); +int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, + const double Cvinf, const double v0, const double T0, const double Gc); + int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const double s1, const double s2, const double s3, const double G0, const double b, const double rho0, const double T0, const double P0, diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index f7bc9f0ecf0..df8e7628a52 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -31,29 +31,47 @@ using EOS = singularity::Variant; constexpr Real REAL_TOL = std::numeric_limits::epsilon() * 1e4; // 1e-12 for double +// Parameters for copper +constexpr Real A = 7.3; +constexpr Real B = 3.9; +constexpr Real Cvinf = 0.000389; +constexpr Real v0 = 1. / 8.952; +constexpr Real T0 = 150.; +constexpr Real Gc = 0.5; + +// Create the EOS +auto host_eos = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); +auto eos = host_eos.GetOnDevice(); + // Only run exception checking when we aren't offloading to the GPUs SCENARIO("Zero temperature on cold curve", "[SimpleMACAWEOS][Temperature]") { GIVEN("Parameters for a Simple MACAW EOS") { - // Unit conversions - //constexpr Real cm = 1.; - //constexpr Real us = 1e-06; - //constexpr Real Mbcc_per_g = 1e12; - // Gruneisen parameters for copper - constexpr Real A = 7.3; - constexpr Real B = 3.9; - constexpr Real Cvinf = 0.000389; - constexpr Real v0 = 1. / 8.952; - constexpr Real T0 = 150.; - constexpr Real Gc = 0.5; - // Create the EOS - EOS host_eos = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); - //EOS eos = host_eos.GetOnDevice(); Real rho = 0.5; for (int i = 0; i < 10; i++) { rho += rho + i; // cylce through a variety of densities THEN("The temperature evaluated on the cold curve should produce zero") { - REQUIRE(eos.TemperatureFromDensityInternalEnergy(rho, eos.SieColdCurve(1.0 / rho))); + Real v = 1.0 / rho; + Real e = eos.SieColdCurve(v); + REQUIRE(eos.TemperatureFromDensityInternalEnergy(rho, e) == 0.0); } } } } + +// Only run exception checking when we aren't offloading to the GPUs +SCENARIO("Inversion equivalence", "[SimpleMACAWEOS][Temperature]") { + GIVEN("Initial density and specific internal energy") { + const Real rho_0 = 0.3; + const Real sie_0 = 2.7; + const Real P = eos.PressureFromDensityInternalEnergy(rho_0, sie_0); + const Real T = eos.TemperatureFromDensityInternalEnergy(rho_0, sie_0); + Real rho, sie; + Real* lambda; + eos.DensityEnergyFromPressureTemperature(P, T, lambda, rho, sie); + sie = eos.InternalEnergyFromDensityPressure(rho, P); + THEN("The inversion back to rho and sie from P and T should be identital.") { + REQUIRE(isClose(rho, rho_0, REAL_TOL * rho_0)); + REQUIRE(isClose(sie, sie_0, REAL_TOL * sie_0)); + } + } +} From 5cb1f734899966ce4861dd8d4099cd020538884c Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Mon, 8 Sep 2025 12:26:08 -0600 Subject: [PATCH 05/25] Cleaned up the code some more and added another test. --- singularity-eos/eos/eos_simple_macaw.hpp | 60 ++++++++++++++---------- test/test_eos_simple_macaw.cpp | 18 ++++++- 2 files changed, 52 insertions(+), 26 deletions(-) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 135b214ddea..1ad6709907d 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -33,6 +33,8 @@ namespace singularity { using namespace eos_base; +/* The details of this equation of state can be found here: + * https://www.osti.gov/biblio/2479474 */ class SimpleMACAW : public EosBase { public: SimpleMACAW() = default; @@ -48,7 +50,7 @@ class SimpleMACAW : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (18) */ - const Real Delta_e = sie - SieColdCurve(1.0 / rho); + const Real Delta_e = sie - SieColdCurve(rho); const Real discriminant = Delta_e * (Delta_e + 4.0 * _Cvinf * _T0 * std::pow(rho * _v0, _Gc)); return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * _Cvinf); } @@ -58,35 +60,32 @@ class SimpleMACAW : public EosBase { PORTABLE_ALWAYS_REQUIRE(_v0 > 0, "Reference specific volume, 'v0', must be positive"); PORTABLE_ALWAYS_REQUIRE(_T0 > 0, "Reference temperature, 'T0', must be positive"); PORTABLE_ALWAYS_REQUIRE(_Cvinf > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); - PORTABLE_ALWAYS_REQUIRE(_Gc > 0 && _Gc < 1, "Gruneisen coefficient, 'Gc', must be in the interval (0,1)"); + PORTABLE_ALWAYS_REQUIRE(_Gc > 0 && _Gc < 1, "Gruneisen parameter, 'Gc', must be in the interval (0,1)"); _AZbar.CheckParams(); } template PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real v = 1.0 / rho; - return rho * (SieColdCurve(v) + SieThermalPortion(v, temperature)); + return rho * (SieColdCurve(rho) + SieThermalPortion(rho, temperature)); } template PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return PressureColdCurve(1.0 / rho) + PressureThermalPortion(rho, temperature); + return PressureColdCurve(rho) + PressureThermalPortion(rho, temperature); } template PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real v = 1.0 / rho; - return PressureColdCurve(v) + _Gc * rho * (sie - SieColdCurve(v)); + return PressureColdCurve(rho) + _Gc * rho * (sie - SieColdCurve(rho)); } template PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityPressure( const Real rho, const Real P, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real v = robust::ratio(1.0, rho); - return SieColdCurve(v) + robust::ratio(v * (P - PressureColdCurve(v)), _Gc); + return SieColdCurve(rho) + robust::ratio(P - PressureColdCurve(rho), _Gc * rho); } // Entropy member functions @@ -106,24 +105,27 @@ class SimpleMACAW : public EosBase { } // Cold curve, thermal portion, and other helper functions + + /* Note: The cold curves for the simple MACAW are presented as functions of `v` rather than `rho`. + * We keep them as functions of `rho` for computational efficiency. */ template PORTABLE_INLINE_FUNCTION Real SieColdCurve( - const Real v, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real x = robust::ratio(v, _v0); - return _A * _v0 * (std::pow(x, -_B) + _B * x - (_B + 1.0)); + const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { + const Real ratio = rho * _v0; + return _A * _v0 * (std::pow(ratio, _B) + robust::ratio(_B, ratio) - (_B + 1.)); } template PORTABLE_INLINE_FUNCTION Real SieThermalPortion( - const Real v, const Real T, + const Real rho, const Real T, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real x = robust::ratio(v, _v0); - return _Cvinf * T * (1.0 - _T0 / (_T0 + T * std::pow(x, _Gc))); + const Real ratio = rho * _v0; + return _Cvinf * T * (1.0 - robust::ratio(_T0, _T0 + T * std::pow(ratio, -_Gc))); } template PORTABLE_INLINE_FUNCTION Real PressureColdCurve( - const Real v, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real ratio = robust::ratio(v, _v0); - return _A * _B * (std::pow(ratio, -(_B + 1.0)) - 1.0); + const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { + const Real ratio = rho * _v0; + return _A * _B * (std::pow(ratio, _B + 1.0) - 1.0); } template PORTABLE_INLINE_FUNCTION Real PressureThermalPortion( @@ -169,7 +171,7 @@ class SimpleMACAW : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0); - const Real term2 = _Gc * (_Gc + 1.0) * rho * (sie - SieColdCurve(1.0 / rho)); + const Real term2 = _Gc * (_Gc + 1.0) * rho * (sie - SieColdCurve(rho)); return term1 + term2; } // Isothermal bulk modulus @@ -182,7 +184,17 @@ class SimpleMACAW : public EosBase { const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0); const Real numerator = rho * _T0 * (1.0 - _Gc) * std::pow(rho * _v0, 2.0 * _Gc) + temperature; const Real denominator = _T0 * std::pow(rho * _v0, _Gc) + temperature; - return term1 + _Cvinf * _Gc * temperature * temperature * numerator / (denominator * denominator); + return term1 + _Cvinf * _Gc * temperature * temperature * robust::ratio(numerator, denominator * denominator); + } + // Specific heat capacity at constant pressure + template + PORTABLE_INLINE_FUNCTION Real ConstantPressureSpecificHeatFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real BT = IsothermalBulkModulusFromDensityTemperature(rho, temperature); + const Real Bs = BulkModulusFromDensityTemperature(rho, temperature); + const Real cv = SpecificHeatFromDensityTemperature(rho, temperature); + return robust::ratio(Bs * cv, BT); /* General thermodynamic identity */ } // Gruneisen parameter @@ -208,15 +220,15 @@ class SimpleMACAW : 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 { - // use STP: 1 atmosphere, room temperature - // TODO: Fix this subroutine... + // v_0 and T_0 are the only user selected reference state parameters rho = 1.0 / _v0; temp = _T0; - sie = 0.0; press = PressureFromDensityTemperature(rho, temp); + sie = InternalEnergyFromDensityTemperature(rho, temp); + cv = SpecificHeatFromDensityInternalEnergy(rho, sie); bmod = BulkModulusFromDensityInternalEnergy(rho, sie); - dpde = _Gc * rho; + dpde = rho * GruneisenParamFromDensityTemperature(rho, temp); } // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index df8e7628a52..4a6c54034b7 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -58,8 +58,7 @@ SCENARIO("Zero temperature on cold curve", "[SimpleMACAWEOS][Temperature]") { } } -// Only run exception checking when we aren't offloading to the GPUs -SCENARIO("Inversion equivalence", "[SimpleMACAWEOS][Temperature]") { +SCENARIO("Inversion equivalence", "[SimpleMACAWEOS][PT-RhoSie Equivalence]") { GIVEN("Initial density and specific internal energy") { const Real rho_0 = 0.3; const Real sie_0 = 2.7; @@ -75,3 +74,18 @@ SCENARIO("Inversion equivalence", "[SimpleMACAWEOS][Temperature]") { } } } + +SCENARIO("Thermodynamic consistency", "[SimpleMACAWEOS][Thermo Consistency]") { + GIVEN("Initial density and temperature") { + const Real rho_0 = 3.; + const Real T_0 = 178.; + const Real Bs = eos.BulkModulusFromDensityTemperature(rho_0, T_0); + const Real BT = eos.IsothermalBulkModulusFromDensityTemperature(rho_0, T_0); + const Real cv = eos.SpecificHeatFromDensityTemperature(rho_0, T_0); + const Real cp = eos.ConstantPressureSpecificHeatFromDensityTemperature(rho_0, T_0); + THEN("Thermodynamic consistency is: Bs >= BT and cp >= cv.") { + REQUIRE(Bs >= BT); + REQUIRE(cp >= cv); + } + } +} From d4fdcd165e9f33224b200bf2705f821570947cbf Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Mon, 8 Sep 2025 15:33:38 -0600 Subject: [PATCH 06/25] Added documentation for the simple MACAW some more comments. --- doc/sphinx/src/models.rst | 54 ++++++++++++++++++++++++ singularity-eos/eos/eos_simple_macaw.hpp | 24 ++++++++--- 2 files changed, 72 insertions(+), 6 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 0bc2cf1350f..3848b1f9340 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -21,6 +21,8 @@ .. _PowerMG: https://www.osti.gov/biblio/1762624 +.. _SimpleMACAW: https://www.osti.gov/biblio/2479474 + EOS Models =========== @@ -1320,6 +1322,58 @@ This constructor also optionally accepts `MeanAtomicProperties` for the atomic mass and number as a final optional parameter. +.. _MACAW EOS: https://pubs.aip.org/aip/jap/article/134/12/125102/2912256 + +Simple MACAW +```````````` +The `Simple MACAW EOS <_SimpleMACAW>`_ is a simplified version of the `MACAW EOS `_ +and is thermodynamically consistent. It is constructed from a the Helmholtz +free energy using a Murnaghan cold curve and a simplified thermal component +from the MACAW EOS. + +Fundamentally, the equation of state can be written in Mie-Gruneisen form (with constant Gruneisen parameter) as: + +.. math:: + + P(v, e) = P_{\text{cold}}(v) + \Gamma_c \rho (e - e_{\text{cold}}(v)) + +where the cold curves are given by: + +.. math:: + + e_{\text{cold}}(v) = A v_0 \Big[ \Big( \frac{v}{v_0} \Big)^{-B} + \Big( \frac{v}{v_0} \Big) B - (B+1) \Big] + +and + +.. math:: + + p_{\text{cold}}(v) := -e'_{\text{cold}}(v) = AB \Big( \Big( \frac{v}{v_0} \Big)^{-(B+1)} - 1 \Big) + +The specific heat capacity at constant volume for this model is given by, + +.. math:: + + C_v(v, \tau) = C^\infty_v \frac{\tau^2 + 2\tau}{(\tau + 1)^2} \quad \text{where } \tau = \frac{T}{\theta(v)} \quad \text{ and } \quad \theta(v) := T_0 \Big( \frac{v}{v_c} \Big)^{-\Gamma_c} + +Note it obeys the expected physical behavior; that, :math:`\lim_{T\to 0^+} C_v(v,\tau(v,T)) = 0` and +:math:`\lim_{T\to\infty} C_v(v,\tau(v,T) = C^\infty_v < \infty` (Dulong-Petit law). + +The constructor for the Simple MACAW EOS is + +.. code-block:: cpp + + SimpleMACAW(const Real A, const Real B, const Real Cvinf, const Real v0, + const Real T0, const Real Gc) + +where ``A`` is :math:`A`, ``B`` is :math:`B`, ``Cvinf`` is :math:`C^\infty_v`, +``v0`` is :math:`v_0`, ``T0`` is :math:`T_0`, ``Gc`` is :math:`\Gamma_c`. + +In order to maintain thermodynamic consistency, a sufficient set of constraints +is given by :math:`A > 0`, :math:`B > 0`, :math:`C^\infty_v > 0`, :math:`v_0 > +0`, :math:`T_0 > 0`, and :math:`\Gamma_c > 0`. If one wants to maintain +thermodynamic consistency, we also require :math:`\Gamma_c \in (0,1)`. + + JWL EOS `````````` diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 1ad6709907d..e1f465ffecd 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -33,8 +33,12 @@ namespace singularity { using namespace eos_base; -/* The details of this equation of state can be found here: - * https://www.osti.gov/biblio/2479474 */ +/* Most information regarding this equation of state can be found here: + * https://www.osti.gov/biblio/2479474 + * + * Note that all equations are given as functions of the density ($\rho$) rather than functions + * of the specific volume ($v$). The reason for this is it simplifies the computational complexity. + */ class SimpleMACAW : public EosBase { public: SimpleMACAW() = default; @@ -55,10 +59,10 @@ class SimpleMACAW : public EosBase { return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * _Cvinf); } PORTABLE_INLINE_FUNCTION void CheckParams() const { - PORTABLE_ALWAYS_REQUIRE(_A > 0, "Parameter, 'A', must be positive"); - PORTABLE_ALWAYS_REQUIRE(_B > 0, "Parameter, 'B', must be positive"); + PORTABLE_ALWAYS_REQUIRE(_A >= 0, "Parameter, 'A', must be non-negative"); + PORTABLE_ALWAYS_REQUIRE(_B >= 0, "Parameter, 'B', must be non-negative"); PORTABLE_ALWAYS_REQUIRE(_v0 > 0, "Reference specific volume, 'v0', must be positive"); - PORTABLE_ALWAYS_REQUIRE(_T0 > 0, "Reference temperature, 'T0', must be positive"); + PORTABLE_ALWAYS_REQUIRE(_T0 >= 0, "Reference temperature, 'T0', must be non-negative"); PORTABLE_ALWAYS_REQUIRE(_Cvinf > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); PORTABLE_ALWAYS_REQUIRE(_Gc > 0 && _Gc < 1, "Gruneisen parameter, 'Gc', must be in the interval (0,1)"); _AZbar.CheckParams(); @@ -67,18 +71,21 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (16) */ return rho * (SieColdCurve(rho) + SieThermalPortion(rho, temperature)); } template PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (15) */ return PressureColdCurve(rho) + PressureThermalPortion(rho, temperature); } template PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (19) */ return PressureColdCurve(rho) + _Gc * rho * (sie - SieColdCurve(rho)); } template @@ -93,6 +100,7 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (8) */ const Real tau = robust::ratio(temperature, TemperatureScale(rho)); return _Cvinf * (robust::ratio(tau, 1.0 + tau) + std::log(1.0 + tau)); } @@ -138,6 +146,7 @@ class SimpleMACAW : public EosBase { template PORTABLE_INLINE_FUNCTION Real TemperatureScale( const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (13) */ return _T0 * std::pow(rho * _v0, _Gc); } @@ -146,6 +155,7 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (6) */ const Real tau = robust::ratio(temperature, TemperatureScale(rho)); const Real numerator = math_utils::pow<2>(tau) + 2.0 * tau; const Real denominator = math_utils::pow<2>(tau + 1.0); @@ -170,6 +180,7 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (21) (multiplied by $\rho$ to get bulk mod) */ const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0); const Real term2 = _Gc * (_Gc + 1.0) * rho * (sie - SieColdCurve(rho)); return term1 + term2; @@ -179,7 +190,8 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION Real IsothermalBulkModulusFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - /* As mentioned in the paper, from the constraints on the parameters, + /* Equation (22) (multiplied by $\rho$ to get bulk mod) + * As mentioned in the paper, from the constraints on the parameters, * the isothermal bulk modulus is guaranteed to be positive. */ const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0); const Real numerator = rho * _T0 * (1.0 - _Gc) * std::pow(rho * _v0, 2.0 * _Gc) + temperature; From bd82905efcd84f9e99c6be0ef6eccb6a63865a30 Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Tue, 9 Sep 2025 10:02:38 -0600 Subject: [PATCH 07/25] Reordered the parameters to fix a clang sanitizer issue. --- CMakeLists.txt | 2 +- singularity-eos/eos/eos_simple_macaw.hpp | 10 ++++------ 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e17e789dd16..4a8168b91fb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -596,7 +596,7 @@ target_compile_options( ) if (SINGULARITY_STRICT_WARNINGS) target_compile_options(singularity-eos_Interface INTERFACE - -Wall -Werror -Wno-unknown-pragmas) + -Wall -Werror -Wno-unknown-pragmas) endif() if(TARGET singularity-eos_Library) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index e1f465ffecd..29efe516133 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -61,9 +61,9 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION void CheckParams() const { PORTABLE_ALWAYS_REQUIRE(_A >= 0, "Parameter, 'A', must be non-negative"); PORTABLE_ALWAYS_REQUIRE(_B >= 0, "Parameter, 'B', must be non-negative"); + PORTABLE_ALWAYS_REQUIRE(_Cvinf > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); PORTABLE_ALWAYS_REQUIRE(_v0 > 0, "Reference specific volume, 'v0', must be positive"); PORTABLE_ALWAYS_REQUIRE(_T0 >= 0, "Reference temperature, 'T0', must be non-negative"); - PORTABLE_ALWAYS_REQUIRE(_Cvinf > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); PORTABLE_ALWAYS_REQUIRE(_Gc > 0 && _Gc < 1, "Gruneisen parameter, 'Gc', must be in the interval (0,1)"); _AZbar.CheckParams(); } @@ -249,9 +249,9 @@ class SimpleMACAW : public EosBase { static constexpr unsigned long PreferredInput() { return _preferred_input; } PORTABLE_INLINE_FUNCTION void PrintParams() const { - printf("Simple MACAW Parameters:\nA = %g\nB = %g\nCvinf = %g\nGc = " + printf("Simple MACAW Parameters:\nA = %g\nB = %g\nCvinf = %g\nv0 = %g\nT0 = %g\nGc = " "%g\n", - _A, _B, _Cvinf, _Gc); + _A, _B, _Cvinf, _v0, _T0, _Gc); _AZbar.PrintParams(); } template @@ -287,9 +287,7 @@ class SimpleMACAW : public EosBase { static std::string EosPyType() { return EosType(); } private: - Real _A, _B, _Cvinf, _Gc; - // reference values - Real _v0, _T0; + Real _A, _B, _Cvinf, _v0, _T0, _Gc; MeanAtomicProperties _AZbar; static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; From 7e7f43250feedfed04e72a0878d0a4a4422e49b3 Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Mon, 15 Sep 2025 22:41:59 -0600 Subject: [PATCH 08/25] Made suggested changes. A few more thing need to be finalized. --- doc/sphinx/src/models.rst | 3 +- singularity-eos/CMakeLists.txt | 1 + singularity-eos/eos/eos_simple_macaw.hpp | 59 ++++++------- test/test_eos_simple_macaw.cpp | 102 +++++++++++++---------- 4 files changed, 85 insertions(+), 80 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 3848b1f9340..f910e1b3fa3 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -23,6 +23,7 @@ .. _SimpleMACAW: https://www.osti.gov/biblio/2479474 +.. _MACAW EOS: https://pubs.aip.org/aip/jap/article/134/12/125102/2912256 EOS Models =========== @@ -1322,8 +1323,6 @@ This constructor also optionally accepts `MeanAtomicProperties` for the atomic mass and number as a final optional parameter. -.. _MACAW EOS: https://pubs.aip.org/aip/jap/article/134/12/125102/2912256 - Simple MACAW ```````````` The `Simple MACAW EOS <_SimpleMACAW>`_ is a simplified version of the `MACAW EOS `_ diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 81302eb703c..51b71e5a4af 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -53,6 +53,7 @@ register_headers( eos/eos_gruneisen.hpp eos/eos_vinet.hpp eos/eos_builder.hpp + eos/eos_simple_macaw.hpp eos/eos_jwl.hpp eos/eos_helmholtz.hpp eos/eos_sap_polynomial.hpp diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 29efe516133..3d088f308c4 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -64,7 +64,8 @@ class SimpleMACAW : public EosBase { PORTABLE_ALWAYS_REQUIRE(_Cvinf > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); PORTABLE_ALWAYS_REQUIRE(_v0 > 0, "Reference specific volume, 'v0', must be positive"); PORTABLE_ALWAYS_REQUIRE(_T0 >= 0, "Reference temperature, 'T0', must be non-negative"); - PORTABLE_ALWAYS_REQUIRE(_Gc > 0 && _Gc < 1, "Gruneisen parameter, 'Gc', must be in the interval (0,1)"); + PORTABLE_ALWAYS_REQUIRE(_Gc > 0, "Gruneisen parameter, 'Gc', must be positive"); + if (_Gc > 1) {PORTABLE_WARN("Warning: Gruneisen coefficient greater than 1. Thermodynamic stability may be violated.");} _AZbar.CheckParams(); } template @@ -224,9 +225,31 @@ class SimpleMACAW : public EosBase { } template PORTABLE_INLINE_FUNCTION void - FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, - const unsigned long output, - Indexer_t &&lambda = static_cast(nullptr)) const; + FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, + const unsigned long output, Indexer_t &&lambda) const { + if (output & thermalqs::density && output & thermalqs::specific_internal_energy) { + if (output & thermalqs::pressure || output & thermalqs::temperature) { + UNDEFINED_ERROR; + } + DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie); + } + if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) { + if (output & thermalqs::density || output & thermalqs::temperature) { + UNDEFINED_ERROR; + } + sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); + } + if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) { + sie = InternalEnergyFromDensityPressure(rho, press); + } + if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::temperature) + temp = TemperatureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::bulk_modulus) + bmod = BulkModulusFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::specific_heat) + cv = SpecificHeatFromDensityInternalEnergy(rho, sie); + } template PORTABLE_INLINE_FUNCTION void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, @@ -293,34 +316,6 @@ class SimpleMACAW : public EosBase { thermalqs::density | thermalqs::specific_internal_energy; }; -template -PORTABLE_INLINE_FUNCTION void -SimpleMACAW::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, - const unsigned long output, Indexer_t &&lambda) const { - if (output & thermalqs::density && output & thermalqs::specific_internal_energy) { - if (output & thermalqs::pressure || output & thermalqs::temperature) { - UNDEFINED_ERROR; - } - DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie); - } - if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) { - if (output & thermalqs::density || output & thermalqs::temperature) { - UNDEFINED_ERROR; - } - sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); - } - if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) { - sie = InternalEnergyFromDensityPressure(rho, press); - } - if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); - if (output & thermalqs::temperature) - temp = TemperatureFromDensityInternalEnergy(rho, sie); - if (output & thermalqs::bulk_modulus) - bmod = BulkModulusFromDensityInternalEnergy(rho, sie); - if (output & thermalqs::specific_heat) - cv = SpecificHeatFromDensityInternalEnergy(rho, sie); -} - } // namespace singularity #endif // _SINGULARITY_EOS_EOS_EOS_SIMPLE_MACAW_HPP_ diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index 4a6c54034b7..5dfb6b44828 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -31,61 +31,71 @@ using EOS = singularity::Variant; constexpr Real REAL_TOL = std::numeric_limits::epsilon() * 1e4; // 1e-12 for double -// Parameters for copper -constexpr Real A = 7.3; -constexpr Real B = 3.9; -constexpr Real Cvinf = 0.000389; -constexpr Real v0 = 1. / 8.952; -constexpr Real T0 = 150.; -constexpr Real Gc = 0.5; -// Create the EOS -auto host_eos = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); -auto eos = host_eos.GetOnDevice(); // Only run exception checking when we aren't offloading to the GPUs -SCENARIO("Zero temperature on cold curve", "[SimpleMACAWEOS][Temperature]") { +SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { GIVEN("Parameters for a Simple MACAW EOS") { - Real rho = 0.5; - for (int i = 0; i < 10; i++) { - rho += rho + i; // cylce through a variety of densities - THEN("The temperature evaluated on the cold curve should produce zero") { - Real v = 1.0 / rho; - Real e = eos.SieColdCurve(v); - REQUIRE(eos.TemperatureFromDensityInternalEnergy(rho, e) == 0.0); + // Parameters for copper (See Table 1 in the paper by Aslam & Lozano) + constexpr Real A = 7.3; + constexpr Real B = 3.9; + constexpr Real Cvinf = 0.000389; + constexpr Real v0 = 1. / 8.952; + constexpr Real T0 = 150.; + constexpr Real Gc = 0.5; + + // Create the EOS + auto host_eos = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); + auto eos = host_eos.GetOnDevice(); + + WHEN("A set of densities is provided") { + Real rho = 0.5; + for (int i = 0; i < 10; i++) { + rho += rho + i; // cylce through a variety of densities + DYNAMIC_SECTION("For a given density " << rho << " and energy from the cold curve") { + Real v = 1.0 / rho; + Real e = eos.SieColdCurve(v); + INFO("rho = " << rho << " v = " << v << " e = " << e); + THEN("The temperature at this density and energy should be zero") { + REQUIRE(eos.TemperatureFromDensityInternalEnergy(rho, e) == 0.0); + } // Then + } // Dynamic Section + } // for + } // When + + WHEN("A set of densities and energies are provided") { + const Real rho_0 = 0.3; + const Real sie_0 = 2.7; + const Real P = eos.PressureFromDensityInternalEnergy(rho_0, sie_0); + const Real T = eos.TemperatureFromDensityInternalEnergy(rho_0, sie_0); + Real rho, sie; + Real* lambda; + eos.DensityEnergyFromPressureTemperature(P, T, lambda, rho, sie); + sie = eos.InternalEnergyFromDensityPressure(rho, P); + INFO("rho_0 = " << rho_0 << " rho = " << rho << " e_0 = " << sie_0 << " e = " << sie); + INFO("P = " << P << " T = " << T); + THEN("The inversion back to rho and sie from P and T should be identital.") { + REQUIRE(isClose(rho, rho_0, REAL_TOL * rho_0)); + REQUIRE(isClose(sie, sie_0, REAL_TOL * sie_0)); } - } - } -} + } // When -SCENARIO("Inversion equivalence", "[SimpleMACAWEOS][PT-RhoSie Equivalence]") { - GIVEN("Initial density and specific internal energy") { - const Real rho_0 = 0.3; - const Real sie_0 = 2.7; - const Real P = eos.PressureFromDensityInternalEnergy(rho_0, sie_0); - const Real T = eos.TemperatureFromDensityInternalEnergy(rho_0, sie_0); - Real rho, sie; - Real* lambda; - eos.DensityEnergyFromPressureTemperature(P, T, lambda, rho, sie); - sie = eos.InternalEnergyFromDensityPressure(rho, P); - THEN("The inversion back to rho and sie from P and T should be identital.") { - REQUIRE(isClose(rho, rho_0, REAL_TOL * rho_0)); - REQUIRE(isClose(sie, sie_0, REAL_TOL * sie_0)); + WHEN("") { + const Real rho_0 = 3.; + const Real T_0 = 178.; + const Real Bs = eos.BulkModulusFromDensityTemperature(rho_0, T_0); + const Real BT = eos.IsothermalBulkModulusFromDensityTemperature(rho_0, T_0); + const Real cv = eos.SpecificHeatFromDensityTemperature(rho_0, T_0); + const Real cp = eos.ConstantPressureSpecificHeatFromDensityTemperature(rho_0, T_0); + THEN("Thermodynamic stability should satisfy: Bs >= BT and cp >= cv.") { + REQUIRE(Bs >= BT); + REQUIRE(cp >= cv); + } } - } -} + } // Given +} // Scenario SCENARIO("Thermodynamic consistency", "[SimpleMACAWEOS][Thermo Consistency]") { GIVEN("Initial density and temperature") { - const Real rho_0 = 3.; - const Real T_0 = 178.; - const Real Bs = eos.BulkModulusFromDensityTemperature(rho_0, T_0); - const Real BT = eos.IsothermalBulkModulusFromDensityTemperature(rho_0, T_0); - const Real cv = eos.SpecificHeatFromDensityTemperature(rho_0, T_0); - const Real cp = eos.ConstantPressureSpecificHeatFromDensityTemperature(rho_0, T_0); - THEN("Thermodynamic consistency is: Bs >= BT and cp >= cv.") { - REQUIRE(Bs >= BT); - REQUIRE(cp >= cv); - } } } From fe50dcb3d29f318f9992831baed48e628d1091a8 Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Tue, 16 Sep 2025 12:53:45 -0600 Subject: [PATCH 09/25] I think I got all the corrections done. --- doc/sphinx/src/models.rst | 8 +- singularity-eos/eos/eos_simple_macaw.hpp | 134 ++++++++++++----------- test/test_eos_simple_macaw.cpp | 88 +++++++++------ 3 files changed, 129 insertions(+), 101 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index f910e1b3fa3..31792f09f35 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1367,11 +1367,11 @@ The constructor for the Simple MACAW EOS is where ``A`` is :math:`A`, ``B`` is :math:`B`, ``Cvinf`` is :math:`C^\infty_v`, ``v0`` is :math:`v_0`, ``T0`` is :math:`T_0`, ``Gc`` is :math:`\Gamma_c`. -In order to maintain thermodynamic consistency, a sufficient set of constraints +In order to maintain thermodynamic stability, a sufficient set of constraints is given by :math:`A > 0`, :math:`B > 0`, :math:`C^\infty_v > 0`, :math:`v_0 > -0`, :math:`T_0 > 0`, and :math:`\Gamma_c > 0`. If one wants to maintain -thermodynamic consistency, we also require :math:`\Gamma_c \in (0,1)`. - +0`, :math:`T_0 > 0`, and :math:`\Gamma_c \in (0,1]`. One can still select +:math:`\Gamma_c > 1`, just note that the isothermal bulk modulus can be +negative (the isentropic bulk modulus will still be positive though). JWL EOS `````````` diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 3d088f308c4..562d3851d5d 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -45,7 +45,7 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION SimpleMACAW(Real A, Real B, Real Cvinf, Real v0, Real T0, Real Gc, const MeanAtomicProperties &AZbar = MeanAtomicProperties()) - : _A(A), _B(B), _Cvinf(Cvinf), _v0(v0), _T0(T0), _Gc(Gc), _AZbar(AZbar) { + : A_(A), B_(B), Cvinf_(Cvinf), v0_(v0), T0_(T0), Gc_(Gc), _AZbar(AZbar) { CheckParams(); } SimpleMACAW GetOnDevice() { return *this; } @@ -55,17 +55,17 @@ class SimpleMACAW : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (18) */ const Real Delta_e = sie - SieColdCurve(rho); - const Real discriminant = Delta_e * (Delta_e + 4.0 * _Cvinf * _T0 * std::pow(rho * _v0, _Gc)); - return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * _Cvinf); + const Real discriminant = Delta_e * (Delta_e + 4.0 * Cvinf_ * T0_ * std::pow(rho * v0_, Gc_)); + return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * Cvinf_); } PORTABLE_INLINE_FUNCTION void CheckParams() const { - PORTABLE_ALWAYS_REQUIRE(_A >= 0, "Parameter, 'A', must be non-negative"); - PORTABLE_ALWAYS_REQUIRE(_B >= 0, "Parameter, 'B', must be non-negative"); - PORTABLE_ALWAYS_REQUIRE(_Cvinf > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); - PORTABLE_ALWAYS_REQUIRE(_v0 > 0, "Reference specific volume, 'v0', must be positive"); - PORTABLE_ALWAYS_REQUIRE(_T0 >= 0, "Reference temperature, 'T0', must be non-negative"); - PORTABLE_ALWAYS_REQUIRE(_Gc > 0, "Gruneisen parameter, 'Gc', must be positive"); - if (_Gc > 1) {PORTABLE_WARN("Warning: Gruneisen coefficient greater than 1. Thermodynamic stability may be violated.");} + PORTABLE_ALWAYS_REQUIRE(A_ >= 0, "Parameter, 'A', must be non-negative"); + PORTABLE_ALWAYS_REQUIRE(B_ >= 0, "Parameter, 'B', must be non-negative"); + PORTABLE_ALWAYS_REQUIRE(Cvinf_ > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); + PORTABLE_ALWAYS_REQUIRE(v0_ > 0, "Reference specific volume, 'v0', must be positive"); + PORTABLE_ALWAYS_REQUIRE(T0_ >= 0, "Reference temperature, 'T0', must be non-negative"); + PORTABLE_ALWAYS_REQUIRE(Gc_ > 0, "Gruneisen parameter, 'Gc', must be positive"); + if (Gc_ > 1) {PORTABLE_WARN("Warning: Gruneisen coefficient greater than 1. Thermodynamic stability may be violated.");} _AZbar.CheckParams(); } template @@ -73,27 +73,27 @@ class SimpleMACAW : public EosBase { const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (16) */ - return rho * (SieColdCurve(rho) + SieThermalPortion(rho, temperature)); + return SieColdCurve(rho) + _SieThermalPortion(rho, temperature); } template PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (15) */ - return PressureColdCurve(rho) + PressureThermalPortion(rho, temperature); + return PressureColdCurve(rho) + _PressureThermalPortion(rho, temperature); } template PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (19) */ - return PressureColdCurve(rho) + _Gc * rho * (sie - SieColdCurve(rho)); + return PressureColdCurve(rho) + Gc_ * rho * (sie - SieColdCurve(rho)); } template PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityPressure( const Real rho, const Real P, Indexer_t &&lambda = static_cast(nullptr)) const { - return SieColdCurve(rho) + robust::ratio(P - PressureColdCurve(rho), _Gc * rho); + return SieColdCurve(rho) + robust::ratio(P - PressureColdCurve(rho), Gc_ * rho); } // Entropy member functions @@ -102,8 +102,8 @@ class SimpleMACAW : public EosBase { EntropyFromDensityTemperature(const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (8) */ - const Real tau = robust::ratio(temperature, TemperatureScale(rho)); - return _Cvinf * (robust::ratio(tau, 1.0 + tau) + std::log(1.0 + tau)); + const Real tau = robust::ratio(temperature, _TemperatureScale(rho)); + return Cvinf_ * (robust::ratio(tau, 1.0 + tau) + std::log(1.0 + tau)); } template PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( @@ -120,35 +120,14 @@ class SimpleMACAW : public EosBase { template PORTABLE_INLINE_FUNCTION Real SieColdCurve( const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real ratio = rho * _v0; - return _A * _v0 * (std::pow(ratio, _B) + robust::ratio(_B, ratio) - (_B + 1.)); - } - template - PORTABLE_INLINE_FUNCTION Real SieThermalPortion( - const Real rho, const Real T, - Indexer_t &&lambda = static_cast(nullptr)) const { - const Real ratio = rho * _v0; - return _Cvinf * T * (1.0 - robust::ratio(_T0, _T0 + T * std::pow(ratio, -_Gc))); + const Real ratio = rho * v0_; + return A_ * v0_ * (std::pow(ratio, B_) + robust::ratio(B_, ratio) - (B_ + 1.)); } template PORTABLE_INLINE_FUNCTION Real PressureColdCurve( const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real ratio = rho * _v0; - return _A * _B * (std::pow(ratio, _B + 1.0) - 1.0); - } - template - PORTABLE_INLINE_FUNCTION Real PressureThermalPortion( - const Real rho, const Real T, - Indexer_t &&lambda = static_cast(nullptr)) const { - const Real numerator = rho * _Cvinf * _Gc * math_utils::pow<2>(T); - const Real denominator = _T0 * std::pow(rho * _v0, _Gc) + T; - return robust::ratio(numerator, denominator); - } - template - PORTABLE_INLINE_FUNCTION Real TemperatureScale( - const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { - /* Equation (13) */ - return _T0 * std::pow(rho * _v0, _Gc); + const Real ratio = rho * v0_; + return A_ * B_ * (std::pow(ratio, B_ + 1.0) - 1.0); } // Specific heat capacity at constant volume @@ -157,10 +136,10 @@ class SimpleMACAW : public EosBase { const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (6) */ - const Real tau = robust::ratio(temperature, TemperatureScale(rho)); + const Real tau = robust::ratio(temperature, _TemperatureScale(rho)); const Real numerator = math_utils::pow<2>(tau) + 2.0 * tau; const Real denominator = math_utils::pow<2>(tau + 1.0); - return _Cvinf * robust::ratio(numerator, denominator); + return Cvinf_ * robust::ratio(numerator, denominator); } template PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( @@ -182,8 +161,8 @@ class SimpleMACAW : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (21) (multiplied by $\rho$ to get bulk mod) */ - const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0); - const Real term2 = _Gc * (_Gc + 1.0) * rho * (sie - SieColdCurve(rho)); + const Real term1 = A_ * B_ * (B_ + 1.0) * std::pow(rho * v0_, B_ + 1.0); + const Real term2 = Gc_ * (Gc_ + 1.0) * rho * (sie - SieColdCurve(rho)); return term1 + term2; } // Isothermal bulk modulus @@ -194,10 +173,10 @@ class SimpleMACAW : public EosBase { /* Equation (22) (multiplied by $\rho$ to get bulk mod) * As mentioned in the paper, from the constraints on the parameters, * the isothermal bulk modulus is guaranteed to be positive. */ - const Real term1 = _A * _B * (_B + 1.0) * std::pow(rho * _v0, _B + 1.0); - const Real numerator = rho * _T0 * (1.0 - _Gc) * std::pow(rho * _v0, 2.0 * _Gc) + temperature; - const Real denominator = _T0 * std::pow(rho * _v0, _Gc) + temperature; - return term1 + _Cvinf * _Gc * temperature * temperature * robust::ratio(numerator, denominator * denominator); + const Real term1 = A_ * B_ * (B_ + 1.0) * std::pow(rho * v0_, B_ + 1.0); + const Real numerator = rho * T0_ * (1.0 - Gc_) * std::pow(rho * v0_, 2.0 * Gc_) + temperature; + const Real denominator = T0_ * std::pow(rho * v0_, Gc_) + temperature; + return term1 + Cvinf_ * Gc_ * temperature * temperature * robust::ratio(numerator, denominator * denominator); } // Specific heat capacity at constant pressure template @@ -215,28 +194,26 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return _Gc; + return Gc_; } template PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - return _Gc; + return Gc_; } template PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda) const { if (output & thermalqs::density && output & thermalqs::specific_internal_energy) { - if (output & thermalqs::pressure || output & thermalqs::temperature) { - UNDEFINED_ERROR; - } + PORTABLE_ALWAYS_REQUIRE(!(output & thermalqs::pressure || output & thermalqs::temperature), + "Cannot request more than two thermodynamic variables (rho, T, e, P)"); DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie); } if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) { - if (output & thermalqs::density || output & thermalqs::temperature) { - UNDEFINED_ERROR; - } + PORTABLE_ALWAYS_REQUIRE(!(output & thermalqs::density || output & thermalqs::temperature), + "Cannot request more than two thermodynamic variables (rho, T, e, P)"); sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); } if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) { @@ -256,8 +233,8 @@ class SimpleMACAW : public EosBase { Real &bmod, Real &dpde, Real &dvdt, Indexer_t &&lambda = static_cast(nullptr)) const { // v_0 and T_0 are the only user selected reference state parameters - rho = 1.0 / _v0; - temp = _T0; + rho = 1.0 / v0_; + temp = T0_; press = PressureFromDensityTemperature(rho, temp); sie = InternalEnergyFromDensityTemperature(rho, temp); @@ -274,7 +251,7 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION void PrintParams() const { printf("Simple MACAW Parameters:\nA = %g\nB = %g\nCvinf = %g\nv0 = %g\nT0 = %g\nGc = " "%g\n", - _A, _B, _Cvinf, _v0, _T0, _Gc); + A_, B_, Cvinf_, v0_, T0_, Gc_); _AZbar.PrintParams(); } template @@ -286,19 +263,19 @@ class SimpleMACAW : public EosBase { // Setup lambda function for rho auto f = [=](const Real x /* density */){ - const Real term1 = _A * _B * (std::pow(_v0 * x, _B + 1.0) - 1.0) - press; - const Real term2 = _T0 * std::pow(_v0 * x, _Gc) + temp; - const Real term3 = _Cvinf * _Gc * temp * temp * x; + const Real term1 = A_ * B_ * (std::pow(v0_ * x, B_ + 1.0) - 1.0) - press; + const Real term2 = T0_ * std::pow(v0_ * x, Gc_) + temp; + const Real term3 = Cvinf_ * Gc_ * temp * temp * x; return term1 * term2 + term3; }; const RootFinding1D::RootCounts root_info; // Finding the density on pressure cold curve is a guaranteed upper bound - const Real rho_high = std::pow(press / (_A * _B) + 1.0, 1.0 / (_B + 1.0)) / _v0; + const Real rho_high = std::pow(press / (A_ * B_) + 1.0, 1.0 / (B_ + 1.0)) / v0_; const Real rho_low = 0.0; // zero density is valid for `f` defined above. - regula_falsi(f, 0.0 /*target*/, 1.0 / _v0 /*guess*/, + regula_falsi(f, 0.0 /*target*/, 1.0 / v0_ /*guess*/, rho_low /*left bracket*/, rho_high /*right bracket*/, 1.0e-9 /*x? tol*/, 1.0e-9 /*y? tol*/, rho, &root_info, true); @@ -310,10 +287,35 @@ class SimpleMACAW : public EosBase { static std::string EosPyType() { return EosType(); } private: - Real _A, _B, _Cvinf, _v0, _T0, _Gc; + Real A_, B_, Cvinf_, v0_, T0_, Gc_; MeanAtomicProperties _AZbar; static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; + + template + PORTABLE_INLINE_FUNCTION Real _TemperatureScale( + const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (13) */ + return T0_ * std::pow(rho * v0_, Gc_); + } + + template + PORTABLE_INLINE_FUNCTION Real _SieThermalPortion( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real ratio = rho * v0_; + return Cvinf_ * T * (1.0 - robust::ratio(T0_, T0_ + T * std::pow(ratio, -Gc_))); + } + + template + PORTABLE_INLINE_FUNCTION Real _PressureThermalPortion( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real numerator = rho * Cvinf_ * Gc_ * math_utils::pow<2>(T); + const Real denominator = T0_ * std::pow(rho * v0_, Gc_) + T; + return robust::ratio(numerator, denominator); + } + }; } // namespace singularity diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index 5dfb6b44828..39b2a033b0f 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -29,9 +29,7 @@ using singularity::SimpleMACAW; using EOS = singularity::Variant; -constexpr Real REAL_TOL = std::numeric_limits::epsilon() * 1e4; // 1e-12 for double - - +constexpr Real REAL_TOL = 1e-12; // Only run exception checking when we aren't offloading to the GPUs SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { @@ -53,49 +51,77 @@ SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { for (int i = 0; i < 10; i++) { rho += rho + i; // cylce through a variety of densities DYNAMIC_SECTION("For a given density " << rho << " and energy from the cold curve") { - Real v = 1.0 / rho; - Real e = eos.SieColdCurve(v); - INFO("rho = " << rho << " v = " << v << " e = " << e); + Real e = eos.SieColdCurve(rho); + INFO("rho = " << rho << " e = " << e); THEN("The temperature at this density and energy should be zero") { REQUIRE(eos.TemperatureFromDensityInternalEnergy(rho, e) == 0.0); } // Then } // Dynamic Section } // for } // When + host_eos.Finalize(); + eos.Finalize(); + } // Given +} // Scenario + +SCENARIO("Testing the Variant API Simple MACAW EOS", "[SimpleMACAWEOS]") { + GIVEN("Parameters for a Simple MACAW EOS") { + // Parameters for copper (See Table 1 in the paper by Aslam & Lozano) + constexpr Real A = 7.3; + constexpr Real B = 3.9; + constexpr Real Cvinf = 0.000389; + constexpr Real v0 = 1. / 8.952; + constexpr Real T0 = 150.; + constexpr Real Gc = 0.5; + + // Create the EOS + using EOS = singularity::Variant; + EOS eos_in_variant = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); + + WHEN("Densities are provided") { + // For large temperatures, the specific heat capacity should approach a constant value + // (The Dulong-Petit Law) + Real rho = 0.1; + Real T = 1e7; + Real Cv = eos_in_variant.SpecificHeatFromDensityTemperature(rho, T); + for (int i = 0; i < 10; i++) { + rho += rho + i; // cylce through a variety of densities + DYNAMIC_SECTION("For a given density " << rho << " and temperature " << T) { + INFO(std::fixed << std::setprecision(15) << + "Cv = " << Cv << ", Cvinf = " << Cvinf << ", T = " << T << ", rho = " << rho); + THEN("The Dulong-Petit law should be satisfied") { + REQUIRE(isClose(Cv, Cvinf, 1e-8 * Cvinf)); + } // Then + } // Dynamic section + } // for + + // Check for zero specific heat at zero temperature + Cv = eos_in_variant.SpecificHeatFromDensityTemperature(rho, 0.); + THEN("The specific heat capacity at constant volume should be zero at zero temperature") { + INFO(std::fixed << std::setprecision(15) << "Cv = " << Cv << ", rho = " << rho); + REQUIRE(Cv == 0.0); + } + } + // Check that the EOS is inverting correctly WHEN("A set of densities and energies are provided") { + using EOS = singularity::Variant; + EOS eos_in_variant = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); + const Real rho_0 = 0.3; - const Real sie_0 = 2.7; - const Real P = eos.PressureFromDensityInternalEnergy(rho_0, sie_0); - const Real T = eos.TemperatureFromDensityInternalEnergy(rho_0, sie_0); + const Real sie_0 = 190.0; + const Real P = eos_in_variant.PressureFromDensityInternalEnergy(rho_0, sie_0); + const Real T = eos_in_variant.TemperatureFromDensityInternalEnergy(rho_0, sie_0); Real rho, sie; Real* lambda; - eos.DensityEnergyFromPressureTemperature(P, T, lambda, rho, sie); - sie = eos.InternalEnergyFromDensityPressure(rho, P); - INFO("rho_0 = " << rho_0 << " rho = " << rho << " e_0 = " << sie_0 << " e = " << sie); + eos_in_variant.DensityEnergyFromPressureTemperature(P, T, lambda, rho, sie); + INFO("rho_0 = " << rho_0 << ", rho = " << rho << ", sie_0 = " << sie_0 << ", sie = " << sie); INFO("P = " << P << " T = " << T); THEN("The inversion back to rho and sie from P and T should be identital.") { REQUIRE(isClose(rho, rho_0, REAL_TOL * rho_0)); - REQUIRE(isClose(sie, sie_0, REAL_TOL * sie_0)); + REQUIRE(isClose(sie, sie_0, 100 * REAL_TOL * sie_0)); } } // When - - WHEN("") { - const Real rho_0 = 3.; - const Real T_0 = 178.; - const Real Bs = eos.BulkModulusFromDensityTemperature(rho_0, T_0); - const Real BT = eos.IsothermalBulkModulusFromDensityTemperature(rho_0, T_0); - const Real cv = eos.SpecificHeatFromDensityTemperature(rho_0, T_0); - const Real cp = eos.ConstantPressureSpecificHeatFromDensityTemperature(rho_0, T_0); - THEN("Thermodynamic stability should satisfy: Bs >= BT and cp >= cv.") { - REQUIRE(Bs >= BT); - REQUIRE(cp >= cv); - } - } + eos_in_variant.Finalize(); } // Given } // Scenario - -SCENARIO("Thermodynamic consistency", "[SimpleMACAWEOS][Thermo Consistency]") { - GIVEN("Initial density and temperature") { - } -} From afc172bca28eb97ad5072c8fc15c87d66be71d72 Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Wed, 17 Sep 2025 11:21:35 -0600 Subject: [PATCH 10/25] Updated copyright and changelog. --- CHANGELOG.md | 1 + singularity-eos/eos/eos_models.hpp | 2 +- singularity-eos/eos/singularity_eos.f90 | 2 +- singularity-eos/eos/singularity_eos.hpp | 2 +- test/test_eos_simple_macaw.cpp | 2 +- 5 files changed, 5 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5cd748f6984..d9252a816cf 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ ### Added (new features/APIs/variables/...) - [[PR556]](https://github.com/lanl/singularity-eos/pull/556) Add introspection into types available in the variant +- [[PR560]](https://github.com/lanl/singularity-eos/pull/560) Add Simple MACAW EOS - [[PR564]](https://github.com/lanl/singularity-eos/pull/564) Removed Get() function from IndexableTypes since it could have unexpected consequences when a type wasn't present ### Fixed (Repair bugs, etc) diff --git a/singularity-eos/eos/eos_models.hpp b/singularity-eos/eos/eos_models.hpp index f2142666cdc..1c6479b029c 100644 --- a/singularity-eos/eos/eos_models.hpp +++ b/singularity-eos/eos/eos_models.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index 71a556206b2..a63627072e8 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -1,5 +1,5 @@ !------------------------------------------------------------------------------ -! © 2021-2024. Triad National Security, LLC. All rights reserved. This +! © 2021-2025. Triad National Security, LLC. All rights reserved. This ! program was produced under U.S. Government contract 89233218CNA000001 ! for Los Alamos National Laboratory (LANL), which is operated by Triad ! National Security, LLC for the U.S. Department of Energy/National diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index d196f01564e..07c75880ad1 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index 39b2a033b0f..c391304fd49 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National From 136ed5e802eb1164751e4cb157e5be44250be64b Mon Sep 17 00:00:00 2001 From: Bennett Clayton Date: Wed, 17 Sep 2025 21:13:26 -0600 Subject: [PATCH 11/25] Fixed a bug in the isothermal bulk modulus which is also a mistake in the paper --- singularity-eos/eos/eos_simple_macaw.hpp | 13 ++++++- test/test_eos_simple_macaw.cpp | 47 ++++++++++++++++++++++++ 2 files changed, 59 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 562d3851d5d..030c7f1f178 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -174,7 +174,8 @@ class SimpleMACAW : public EosBase { * As mentioned in the paper, from the constraints on the parameters, * the isothermal bulk modulus is guaranteed to be positive. */ const Real term1 = A_ * B_ * (B_ + 1.0) * std::pow(rho * v0_, B_ + 1.0); - const Real numerator = rho * T0_ * (1.0 - Gc_) * std::pow(rho * v0_, 2.0 * Gc_) + temperature; + // There is a mistake in the paper and the numerator term should have a -Gc_ instead of -2*Gc_ + const Real numerator = rho * (T0_ * (1.0 - Gc_) * std::pow(rho * v0_, Gc_) + temperature); const Real denominator = T0_ * std::pow(rho * v0_, Gc_) + temperature; return term1 + Cvinf_ * Gc_ * temperature * temperature * robust::ratio(numerator, denominator * denominator); } @@ -188,6 +189,16 @@ class SimpleMACAW : public EosBase { const Real cv = SpecificHeatFromDensityTemperature(rho, temperature); return robust::ratio(Bs * cv, BT); /* General thermodynamic identity */ } + // Coefficient of thermal expansivity + template + PORTABLE_INLINE_FUNCTION Real CoefficientThermalExpansionFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real BT = IsothermalBulkModulusFromDensityTemperature(rho, temperature); + const Real cv = SpecificHeatFromDensityTemperature(rho, temperature); + return robust::ratio(Gc_ * rho * cv, BT); /* General thermodynamic identity */ + } + // Gruneisen parameter template diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index c391304fd49..410f498ad6f 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -59,6 +59,53 @@ SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { } // Dynamic Section } // for } // When + + WHEN("A density and temperature are provided") { + Real rho = 0.56; + Real T = 123.0; + for (int i = 0; i < 10; i++) { + rho += i; // cylce through a variety of densities + T += 15.0 * i; + DYNAMIC_SECTION("For a given density, " << rho << ", and temperature, " << T) { + // Setup thermodynamic derivatives + Real cp = eos.ConstantPressureSpecificHeatFromDensityTemperature(rho, T); + Real cv = eos.SpecificHeatFromDensityTemperature(rho, T); + Real Bs = eos.BulkModulusFromDensityTemperature(rho, T); + Real BT = eos.IsothermalBulkModulusFromDensityTemperature(rho, T); + Real beta = eos.CoefficientThermalExpansionFromDensityTemperature(rho, T); + Real G = eos.GruneisenParamFromDensityTemperature(rho, T); + + // Test a battery of different thermodynamic identity tests + INFO("rho = " << rho << ", T = " << T); + INFO("cp = " << cp << ", cv = " << cv << ", Bs = " << Bs << ", BT = " << BT); + INFO("beta = " << beta << ", G = " << G); + THEN("The thermodynamic stability ensures: cp >= cv and Bs >= BT") { + REQUIRE(Bs >= BT); + REQUIRE(cp >= cv); + } // Then + + Real term1 = cp * BT / (cv * Bs); + INFO("cp * BT / (cv * Bs) = " << term1); + THEN("The thermodynamic equality satisfies: cp * BT / (cv * Bs) = 1") { + REQUIRE(isClose(term1, 1.0, REAL_TOL)); + } // Then + + Real term2 = BT / Bs + beta * beta * T * BT / (rho * cp); + INFO("cv / cp + beta * beta * T * BT / (rho * cp) = " << term2); + THEN("The thermodynamic equality satisfies: cv / cp + beta^2 * T * BT / (rho * cp)") { + REQUIRE(isClose(term2, 1.0, REAL_TOL)); + } // Then + + Real gruneisen = beta * BT / (rho * cv); + INFO("beta * BT / (rho * cv) = " << gruneisen); + THEN("The thermodynamic equality for the Gruneisen parameter satisfies: Gamma = beta * BT / (rho * cv)") { + REQUIRE(isClose(G, gruneisen, REAL_TOL)); + } // Then + + } // Dynamic Section + } // for + } // When + host_eos.Finalize(); eos.Finalize(); } // Given From aba7fc1d7e1f193f2955da4d47d2ae385dcb0934 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 19 Sep 2025 17:31:52 -0400 Subject: [PATCH 12/25] formatting --- singularity-eos/eos/eos_models.hpp | 2 +- singularity-eos/eos/eos_simple_macaw.hpp | 138 +++++++++++--------- singularity-eos/eos/eos_type_lists.hpp | 3 +- singularity-eos/eos/get_sg_eos_functors.hpp | 4 +- singularity-eos/eos/singularity_eos.cpp | 8 +- singularity-eos/eos/singularity_eos.hpp | 7 +- test/test_eos_simple_macaw.cpp | 39 +++--- 7 files changed, 110 insertions(+), 91 deletions(-) diff --git a/singularity-eos/eos/eos_models.hpp b/singularity-eos/eos/eos_models.hpp index 1c6479b029c..497baf1f5bf 100644 --- a/singularity-eos/eos/eos_models.hpp +++ b/singularity-eos/eos/eos_models.hpp @@ -24,11 +24,11 @@ #include #include #include -#include #include #include #include #include +#include #include #include #include diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 030c7f1f178..a0022319ba0 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -34,17 +34,18 @@ namespace singularity { using namespace eos_base; /* Most information regarding this equation of state can be found here: - * https://www.osti.gov/biblio/2479474 - * - * Note that all equations are given as functions of the density ($\rho$) rather than functions - * of the specific volume ($v$). The reason for this is it simplifies the computational complexity. + * https://www.osti.gov/biblio/2479474 + * + * Note that all equations are given as functions of the density ($\rho$) rather than + * functions of the specific volume ($v$). The reason for this is it simplifies the + * computational complexity. */ class SimpleMACAW : public EosBase { public: SimpleMACAW() = default; PORTABLE_INLINE_FUNCTION SimpleMACAW(Real A, Real B, Real Cvinf, Real v0, Real T0, Real Gc, - const MeanAtomicProperties &AZbar = MeanAtomicProperties()) + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : A_(A), B_(B), Cvinf_(Cvinf), v0_(v0), T0_(T0), Gc_(Gc), _AZbar(AZbar) { CheckParams(); } @@ -55,38 +56,45 @@ class SimpleMACAW : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (18) */ const Real Delta_e = sie - SieColdCurve(rho); - const Real discriminant = Delta_e * (Delta_e + 4.0 * Cvinf_ * T0_ * std::pow(rho * v0_, Gc_)); + const Real discriminant = + Delta_e * (Delta_e + 4.0 * Cvinf_ * T0_ * std::pow(rho * v0_, Gc_)); return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * Cvinf_); } PORTABLE_INLINE_FUNCTION void CheckParams() const { PORTABLE_ALWAYS_REQUIRE(A_ >= 0, "Parameter, 'A', must be non-negative"); PORTABLE_ALWAYS_REQUIRE(B_ >= 0, "Parameter, 'B', must be non-negative"); - PORTABLE_ALWAYS_REQUIRE(Cvinf_ > 0, "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); + PORTABLE_ALWAYS_REQUIRE( + Cvinf_ > 0, + "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); PORTABLE_ALWAYS_REQUIRE(v0_ > 0, "Reference specific volume, 'v0', must be positive"); - PORTABLE_ALWAYS_REQUIRE(T0_ >= 0, "Reference temperature, 'T0', must be non-negative"); + PORTABLE_ALWAYS_REQUIRE(T0_ >= 0, + "Reference temperature, 'T0', must be non-negative"); PORTABLE_ALWAYS_REQUIRE(Gc_ > 0, "Gruneisen parameter, 'Gc', must be positive"); - if (Gc_ > 1) {PORTABLE_WARN("Warning: Gruneisen coefficient greater than 1. Thermodynamic stability may be violated.");} + if (Gc_ > 1) { + PORTABLE_WARN("Warning: Gruneisen coefficient greater than 1. Thermodynamic " + "stability may be violated."); + } _AZbar.CheckParams(); } template PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - /* Equation (16) */ + /* Equation (16) */ return SieColdCurve(rho) + _SieThermalPortion(rho, temperature); } template PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - /* Equation (15) */ + /* Equation (15) */ return PressureColdCurve(rho) + _PressureThermalPortion(rho, temperature); } template PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - /* Equation (19) */ + /* Equation (19) */ return PressureColdCurve(rho) + Gc_ * rho * (sie - SieColdCurve(rho)); } template @@ -115,11 +123,11 @@ class SimpleMACAW : public EosBase { // Cold curve, thermal portion, and other helper functions - /* Note: The cold curves for the simple MACAW are presented as functions of `v` rather than `rho`. - * We keep them as functions of `rho` for computational efficiency. */ + /* Note: The cold curves for the simple MACAW are presented as functions of `v` rather + * than `rho`. We keep them as functions of `rho` for computational efficiency. */ template - PORTABLE_INLINE_FUNCTION Real SieColdCurve( - const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { + PORTABLE_INLINE_FUNCTION Real + SieColdCurve(const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { const Real ratio = rho * v0_; return A_ * v0_ * (std::pow(ratio, B_) + robust::ratio(B_, ratio) - (B_ + 1.)); } @@ -148,7 +156,7 @@ class SimpleMACAW : public EosBase { const Real T = TemperatureFromDensityInternalEnergy(rho, sie); return SpecificHeatFromDensityTemperature(rho, T); } - // Isentropic bulk modulus + // Isentropic bulk modulus template PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( const Real rho, const Real temperature, @@ -174,13 +182,16 @@ class SimpleMACAW : public EosBase { * As mentioned in the paper, from the constraints on the parameters, * the isothermal bulk modulus is guaranteed to be positive. */ const Real term1 = A_ * B_ * (B_ + 1.0) * std::pow(rho * v0_, B_ + 1.0); - // There is a mistake in the paper and the numerator term should have a -Gc_ instead of -2*Gc_ - const Real numerator = rho * (T0_ * (1.0 - Gc_) * std::pow(rho * v0_, Gc_) + temperature); + // There is a mistake in the paper and the numerator term should have a -Gc_ instead + // of -2*Gc_ + const Real numerator = + rho * (T0_ * (1.0 - Gc_) * std::pow(rho * v0_, Gc_) + temperature); const Real denominator = T0_ * std::pow(rho * v0_, Gc_) + temperature; - return term1 + Cvinf_ * Gc_ * temperature * temperature * robust::ratio(numerator, denominator * denominator); + return term1 + Cvinf_ * Gc_ * temperature * temperature * + robust::ratio(numerator, denominator * denominator); } // Specific heat capacity at constant pressure - template + template PORTABLE_INLINE_FUNCTION Real ConstantPressureSpecificHeatFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { @@ -190,7 +201,7 @@ class SimpleMACAW : public EosBase { return robust::ratio(Bs * cv, BT); /* General thermodynamic identity */ } // Coefficient of thermal expansivity - template + template PORTABLE_INLINE_FUNCTION Real CoefficientThermalExpansionFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { @@ -199,7 +210,6 @@ class SimpleMACAW : public EosBase { return robust::ratio(Gc_ * rho * cv, BT); /* General thermodynamic identity */ } - // Gruneisen parameter template PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( @@ -214,29 +224,31 @@ class SimpleMACAW : public EosBase { return Gc_; } template - PORTABLE_INLINE_FUNCTION void - FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, - const unsigned long output, Indexer_t &&lambda) const { - if (output & thermalqs::density && output & thermalqs::specific_internal_energy) { - PORTABLE_ALWAYS_REQUIRE(!(output & thermalqs::pressure || output & thermalqs::temperature), - "Cannot request more than two thermodynamic variables (rho, T, e, P)"); - DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie); - } - if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) { - PORTABLE_ALWAYS_REQUIRE(!(output & thermalqs::density || output & thermalqs::temperature), - "Cannot request more than two thermodynamic variables (rho, T, e, P)"); - sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); - } - if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) { - sie = InternalEnergyFromDensityPressure(rho, press); - } - if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); - if (output & thermalqs::temperature) - temp = TemperatureFromDensityInternalEnergy(rho, sie); - if (output & thermalqs::bulk_modulus) - bmod = BulkModulusFromDensityInternalEnergy(rho, sie); - if (output & thermalqs::specific_heat) - cv = SpecificHeatFromDensityInternalEnergy(rho, sie); + PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &sie, Real &press, + Real &cv, Real &bmod, const unsigned long output, + Indexer_t &&lambda) const { + if (output & thermalqs::density && output & thermalqs::specific_internal_energy) { + PORTABLE_ALWAYS_REQUIRE( + !(output & thermalqs::pressure || output & thermalqs::temperature), + "Cannot request more than two thermodynamic variables (rho, T, e, P)"); + DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie); + } + if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) { + PORTABLE_ALWAYS_REQUIRE( + !(output & thermalqs::density || output & thermalqs::temperature), + "Cannot request more than two thermodynamic variables (rho, T, e, P)"); + sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); + } + if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) { + sie = InternalEnergyFromDensityPressure(rho, press); + } + if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::temperature) + temp = TemperatureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::bulk_modulus) + bmod = BulkModulusFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::specific_heat) + cv = SpecificHeatFromDensityInternalEnergy(rho, sie); } template PORTABLE_INLINE_FUNCTION void @@ -260,7 +272,8 @@ class SimpleMACAW : public EosBase { static constexpr unsigned long PreferredInput() { return _preferred_input; } PORTABLE_INLINE_FUNCTION void PrintParams() const { - printf("Simple MACAW Parameters:\nA = %g\nB = %g\nCvinf = %g\nv0 = %g\nT0 = %g\nGc = " + printf("Simple MACAW Parameters:\nA = %g\nB = %g\nCvinf = %g\nv0 = %g\nT0 = " + "%g\nGc = " "%g\n", A_, B_, Cvinf_, v0_, T0_, Gc_); _AZbar.PrintParams(); @@ -269,11 +282,12 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - /* Since the isothermal bulk modulus is always positive (assuming parameter constraints), - * this guarantees the function to solve is monotonic hence we always have a unique solution. */ + /* Since the isothermal bulk modulus is always positive (assuming parameter + * constraints), this guarantees the function to solve is monotonic hence we always + * have a unique solution. */ // Setup lambda function for rho - auto f = [=](const Real x /* density */){ + auto f = [=](const Real x /* density */) { const Real term1 = A_ * B_ * (std::pow(v0_ * x, B_ + 1.0) - 1.0) - press; const Real term2 = T0_ * std::pow(v0_ * x, Gc_) + temp; const Real term3 = Cvinf_ * Gc_ * temp * temp * x; @@ -286,10 +300,9 @@ class SimpleMACAW : public EosBase { const Real rho_high = std::pow(press / (A_ * B_) + 1.0, 1.0 / (B_ + 1.0)) / v0_; const Real rho_low = 0.0; // zero density is valid for `f` defined above. - regula_falsi(f, 0.0 /*target*/, 1.0 / v0_ /*guess*/, - rho_low /*left bracket*/, rho_high /*right bracket*/, - 1.0e-9 /*x? tol*/, 1.0e-9 /*y? tol*/, - rho, &root_info, true); + regula_falsi(f, 0.0 /*target*/, 1.0 / v0_ /*guess*/, rho_low /*left bracket*/, + rho_high /*right bracket*/, 1.0e-9 /*x? tol*/, 1.0e-9 /*y? tol*/, rho, + &root_info, true); sie = InternalEnergyFromDensityTemperature(rho, temp); } @@ -311,22 +324,21 @@ class SimpleMACAW : public EosBase { } template - PORTABLE_INLINE_FUNCTION Real _SieThermalPortion( - const Real rho, const Real T, - Indexer_t &&lambda = static_cast(nullptr)) const { + PORTABLE_INLINE_FUNCTION Real + _SieThermalPortion(const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { const Real ratio = rho * v0_; return Cvinf_ * T * (1.0 - robust::ratio(T0_, T0_ + T * std::pow(ratio, -Gc_))); } template - PORTABLE_INLINE_FUNCTION Real _PressureThermalPortion( - const Real rho, const Real T, - Indexer_t &&lambda = static_cast(nullptr)) const { - const Real numerator = rho * Cvinf_ * Gc_ * math_utils::pow<2>(T); - const Real denominator = T0_ * std::pow(rho * v0_, Gc_) + T; + PORTABLE_INLINE_FUNCTION Real + _PressureThermalPortion(const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real numerator = rho * Cvinf_ * Gc_ * math_utils::pow<2>(T); + const Real denominator = T0_ * std::pow(rho * v0_, Gc_) + T; return robust::ratio(numerator, denominator); } - }; } // namespace singularity diff --git a/singularity-eos/eos/eos_type_lists.hpp b/singularity-eos/eos/eos_type_lists.hpp index 848f80edda0..8895691bd78 100644 --- a/singularity-eos/eos/eos_type_lists.hpp +++ b/singularity-eos/eos/eos_type_lists.hpp @@ -49,7 +49,8 @@ using singularity::variadic_utils::transform_variadic_list; // all eos's static constexpr const auto full_eos_list = - tl= 0); EOS eosi = SGAPPLYMODSIMPLE(SimpleMACAW(A, B, Cvinf, v0, T0, Gc)); if (enabled[3] == 1) { @@ -139,11 +139,11 @@ int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const doub } int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, - const double Cvinf, const double v0, const double T0, const double Gc) { + const double Cvinf, const double v0, const double T0, + const double Gc) { return init_sg_SimpleMACAW(matindex, eos, A, B, Cvinf, v0, T0, Gc, def_en, def_v); } - 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, int const *const enabled, diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index 07c75880ad1..1bfef69fd69 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -35,8 +35,8 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, const double Cv, int const *const enabled, double *const vals); int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, - const double Cvinf, const double v0, const double T0, const double Gc, - int const *const enabled, double *const vals); + const double Cvinf, const double v0, const double T0, + const double Gc, int const *const enabled, double *const vals); int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const double s1, const double s2, const double s3, const double G0, const double b, @@ -149,7 +149,8 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, const double Cv); int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, - const double Cvinf, const double v0, const double T0, const double Gc); + const double Cvinf, const double v0, const double T0, + const double Gc); int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const double s1, const double s2, const double s3, const double G0, const double b, diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index 410f498ad6f..9ef1074a851 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -50,15 +50,16 @@ SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { Real rho = 0.5; for (int i = 0; i < 10; i++) { rho += rho + i; // cylce through a variety of densities - DYNAMIC_SECTION("For a given density " << rho << " and energy from the cold curve") { + DYNAMIC_SECTION("For a given density " << rho + << " and energy from the cold curve") { Real e = eos.SieColdCurve(rho); INFO("rho = " << rho << " e = " << e); THEN("The temperature at this density and energy should be zero") { REQUIRE(eos.TemperatureFromDensityInternalEnergy(rho, e) == 0.0); } // Then - } // Dynamic Section - } // for - } // When + } // Dynamic Section + } // for + } // When WHEN("A density and temperature are provided") { Real rho = 0.56; @@ -92,19 +93,21 @@ SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { Real term2 = BT / Bs + beta * beta * T * BT / (rho * cp); INFO("cv / cp + beta * beta * T * BT / (rho * cp) = " << term2); - THEN("The thermodynamic equality satisfies: cv / cp + beta^2 * T * BT / (rho * cp)") { + THEN("The thermodynamic equality satisfies: cv / cp + beta^2 * T * BT / (rho * " + "cp)") { REQUIRE(isClose(term2, 1.0, REAL_TOL)); } // Then Real gruneisen = beta * BT / (rho * cv); INFO("beta * BT / (rho * cv) = " << gruneisen); - THEN("The thermodynamic equality for the Gruneisen parameter satisfies: Gamma = beta * BT / (rho * cv)") { + THEN("The thermodynamic equality for the Gruneisen parameter satisfies: Gamma " + "= beta * BT / (rho * cv)") { REQUIRE(isClose(G, gruneisen, REAL_TOL)); } // Then } // Dynamic Section - } // for - } // When + } // for + } // When host_eos.Finalize(); eos.Finalize(); @@ -126,25 +129,26 @@ SCENARIO("Testing the Variant API Simple MACAW EOS", "[SimpleMACAWEOS]") { EOS eos_in_variant = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); WHEN("Densities are provided") { - // For large temperatures, the specific heat capacity should approach a constant value - // (The Dulong-Petit Law) + // For large temperatures, the specific heat capacity should approach a constant + // value (The Dulong-Petit Law) Real rho = 0.1; Real T = 1e7; Real Cv = eos_in_variant.SpecificHeatFromDensityTemperature(rho, T); for (int i = 0; i < 10; i++) { rho += rho + i; // cylce through a variety of densities DYNAMIC_SECTION("For a given density " << rho << " and temperature " << T) { - INFO(std::fixed << std::setprecision(15) << - "Cv = " << Cv << ", Cvinf = " << Cvinf << ", T = " << T << ", rho = " << rho); + INFO(std::fixed << std::setprecision(15) << "Cv = " << Cv << ", Cvinf = " + << Cvinf << ", T = " << T << ", rho = " << rho); THEN("The Dulong-Petit law should be satisfied") { REQUIRE(isClose(Cv, Cvinf, 1e-8 * Cvinf)); } // Then - } // Dynamic section - } // for + } // Dynamic section + } // for // Check for zero specific heat at zero temperature Cv = eos_in_variant.SpecificHeatFromDensityTemperature(rho, 0.); - THEN("The specific heat capacity at constant volume should be zero at zero temperature") { + THEN("The specific heat capacity at constant volume should be zero at zero " + "temperature") { INFO(std::fixed << std::setprecision(15) << "Cv = " << Cv << ", rho = " << rho); REQUIRE(Cv == 0.0); } @@ -160,9 +164,10 @@ SCENARIO("Testing the Variant API Simple MACAW EOS", "[SimpleMACAWEOS]") { const Real P = eos_in_variant.PressureFromDensityInternalEnergy(rho_0, sie_0); const Real T = eos_in_variant.TemperatureFromDensityInternalEnergy(rho_0, sie_0); Real rho, sie; - Real* lambda; + Real *lambda; eos_in_variant.DensityEnergyFromPressureTemperature(P, T, lambda, rho, sie); - INFO("rho_0 = " << rho_0 << ", rho = " << rho << ", sie_0 = " << sie_0 << ", sie = " << sie); + INFO("rho_0 = " << rho_0 << ", rho = " << rho << ", sie_0 = " << sie_0 + << ", sie = " << sie); INFO("P = " << P << " T = " << T); THEN("The inversion back to rho and sie from P and T should be identital.") { REQUIRE(isClose(rho, rho_0, REAL_TOL * rho_0)); From 1e5bd4aafe424a0af49e545f70cd0de3bf80c2b1 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 19 Sep 2025 16:47:16 -0600 Subject: [PATCH 13/25] Clang format again? --- singularity-eos/eos/get_sg_eos_functors.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/get_sg_eos_functors.hpp b/singularity-eos/eos/get_sg_eos_functors.hpp index aee80f19338..1cd3b7c0f40 100644 --- a/singularity-eos/eos/get_sg_eos_functors.hpp +++ b/singularity-eos/eos/get_sg_eos_functors.hpp @@ -58,8 +58,8 @@ struct init_functor { frac_vol_v{frac_vol_v_}, frac_ie_v{frac_ie_v_}, pte_mats{pte_mats_}, vfrac_pte{vfrac_pte_}, sie_pte{sie_pte_}, temp_pte{temp_pte_}, press_pte{press_pte_}, rho_pte{rho_pte_}, spvol_v{spvol_v_}, temp_v{temp_v_}, - press_v{press_v_}, sie_v{sie_v_}, nmat{nmat}, mass_frac_cutoff{ - mass_frac_cutoff_} {} + press_v{press_v_}, sie_v{sie_v_}, nmat{nmat}, + mass_frac_cutoff{mass_frac_cutoff_} {} PORTABLE_INLINE_FUNCTION void operator()(const int i, const int tid, double &mass_sum, int &npte, From 6d4776b7dc348dc21c1708bb7ef1d996d32dbf88 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 23 Sep 2025 16:39:10 -0600 Subject: [PATCH 14/25] Floor values at cold curve to maintain EOS surface --- singularity-eos/eos/eos_simple_macaw.hpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index a0022319ba0..cf4a7c9fb59 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -55,7 +55,7 @@ class SimpleMACAW : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (18) */ - const Real Delta_e = sie - SieColdCurve(rho); + const Real Delta_e = std::max(sie - SieColdCurve(rho), 0.); const Real discriminant = Delta_e * (Delta_e + 4.0 * Cvinf_ * T0_ * std::pow(rho * v0_, Gc_)); return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * Cvinf_); @@ -95,13 +95,15 @@ class SimpleMACAW : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (19) */ - return PressureColdCurve(rho) + Gc_ * rho * (sie - SieColdCurve(rho)); + const Real Delta_e = std::max(sie - SieColdCurve(rho), 0.); + return PressureColdCurve(rho) + Gc_ * rho * Delta_e; } template PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityPressure( const Real rho, const Real P, Indexer_t &&lambda = static_cast(nullptr)) const { - return SieColdCurve(rho) + robust::ratio(P - PressureColdCurve(rho), Gc_ * rho); + const Real delta_p = std::max(P - PressureColdCurve(rho), 0.); + return SieColdCurve(rho) + robust::ratio(delta_p, Gc_ * rho); } // Entropy member functions @@ -168,9 +170,10 @@ class SimpleMACAW : public EosBase { PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { + const Real Delta_e = std::max(sie - SieColdCurve(rho), 0.); /* Equation (21) (multiplied by $\rho$ to get bulk mod) */ const Real term1 = A_ * B_ * (B_ + 1.0) * std::pow(rho * v0_, B_ + 1.0); - const Real term2 = Gc_ * (Gc_ + 1.0) * rho * (sie - SieColdCurve(rho)); + const Real term2 = Gc_ * (Gc_ + 1.0) * rho * Delta_e; return term1 + term2; } // Isothermal bulk modulus From e152611250b91ee6bafc49a35ff4077fbc0cbda4 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 23 Sep 2025 19:36:59 -0600 Subject: [PATCH 15/25] Remove THEN from inside DYNAMIC_SELECTION and add WithinRel comparison --- test/test_eos_simple_macaw.cpp | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index 9ef1074a851..49dbf3d317c 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -21,6 +21,7 @@ #define CATCH_CONFIG_FAST_COMPILE #include #endif +#include #include #include @@ -48,18 +49,20 @@ SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { WHEN("A set of densities is provided") { Real rho = 0.5; - for (int i = 0; i < 10; i++) { - rho += rho + i; // cylce through a variety of densities - DYNAMIC_SECTION("For a given density " << rho - << " and energy from the cold curve") { - Real e = eos.SieColdCurve(rho); - INFO("rho = " << rho << " e = " << e); - THEN("The temperature at this density and energy should be zero") { - REQUIRE(eos.TemperatureFromDensityInternalEnergy(rho, e) == 0.0); - } // Then - } // Dynamic Section - } // for - } // When + THEN("The temperatue at this density and an energy on the cold curve should be " + "zero") { + for (int i = 0; i < 10; i++) { + rho += rho + i; // cylce through a variety of densities + DYNAMIC_SECTION("i: " << i << "For a given density " << rho << + " and energy from the cold curve") { + const Real e = eos.SieColdCurve(rho); + INFO("rho = " << rho << " e = " << e); + REQUIRE_THAT(eos.TemperatureFromDensityInternalEnergy(rho, e), + Catch::Matchers::WithinRel(0.0, 1.0e-12)); + } // Dynamic Section + } // for + } // Then + } // When WHEN("A density and temperature are provided") { Real rho = 0.56; From 5909470cd41b8f97ccc0f33067046c5cba208080 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 23 Sep 2025 19:37:37 -0600 Subject: [PATCH 16/25] Clang format --- test/test_eos_simple_macaw.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index 49dbf3d317c..0f34417526d 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -53,16 +53,16 @@ SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { "zero") { for (int i = 0; i < 10; i++) { rho += rho + i; // cylce through a variety of densities - DYNAMIC_SECTION("i: " << i << "For a given density " << rho << - " and energy from the cold curve") { + DYNAMIC_SECTION("i: " << i << "For a given density " << rho + << " and energy from the cold curve") { const Real e = eos.SieColdCurve(rho); INFO("rho = " << rho << " e = " << e); REQUIRE_THAT(eos.TemperatureFromDensityInternalEnergy(rho, e), Catch::Matchers::WithinRel(0.0, 1.0e-12)); - } // Dynamic Section - } // for - } // Then - } // When + } // Dynamic Section + } // for + } // Then + } // When WHEN("A density and temperature are provided") { Real rho = 0.56; From 8a217b50b720c973252b5d1577879649d6d19203 Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Thu, 25 Sep 2025 10:57:36 -0600 Subject: [PATCH 17/25] Remove DYNAMIC_SELECTION macro in case it's the problem --- test/test_eos_simple_macaw.cpp | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index 0f34417526d..d252ced6275 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -53,13 +53,10 @@ SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { "zero") { for (int i = 0; i < 10; i++) { rho += rho + i; // cylce through a variety of densities - DYNAMIC_SECTION("i: " << i << "For a given density " << rho - << " and energy from the cold curve") { - const Real e = eos.SieColdCurve(rho); - INFO("rho = " << rho << " e = " << e); - REQUIRE_THAT(eos.TemperatureFromDensityInternalEnergy(rho, e), - Catch::Matchers::WithinRel(0.0, 1.0e-12)); - } // Dynamic Section + const Real e = eos.SieColdCurve(rho); + INFO("i: " << i << " rho = " << rho << " e = " << e); + REQUIRE_THAT(eos.TemperatureFromDensityInternalEnergy(rho, e), + Catch::Matchers::WithinRel(0.0, 1.0e-12)); } // for } // Then } // When From ca76fbdcd00a2782f2dd455cae1c78c314097160 Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Thu, 25 Sep 2025 11:08:27 -0600 Subject: [PATCH 18/25] Clang format --- test/test_eos_simple_macaw.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index d252ced6275..084eeaf51e0 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -57,9 +57,9 @@ SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { INFO("i: " << i << " rho = " << rho << " e = " << e); REQUIRE_THAT(eos.TemperatureFromDensityInternalEnergy(rho, e), Catch::Matchers::WithinRel(0.0, 1.0e-12)); - } // for - } // Then - } // When + } // for + } // Then + } // When WHEN("A density and temperature are provided") { Real rho = 0.56; From d7656b05b5bff94da2dce37cb3731e6da526ba89 Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Thu, 25 Sep 2025 11:13:09 -0600 Subject: [PATCH 19/25] static cast integer to Real --- test/test_eos_simple_macaw.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index 084eeaf51e0..e3e0af43bca 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -52,7 +52,7 @@ SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { THEN("The temperatue at this density and an energy on the cold curve should be " "zero") { for (int i = 0; i < 10; i++) { - rho += rho + i; // cylce through a variety of densities + rho += rho + static_cast(i); // cylce through a variety of densities const Real e = eos.SieColdCurve(rho); INFO("i: " << i << " rho = " << rho << " e = " << e); REQUIRE_THAT(eos.TemperatureFromDensityInternalEnergy(rho, e), From 1cf972b4e9db16186e0cc5f7f40a243f406c2261 Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Thu, 25 Sep 2025 11:15:13 -0600 Subject: [PATCH 20/25] REQUIRE -> CHECK in loop --- test/test_eos_simple_macaw.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp index e3e0af43bca..f63110a18cb 100644 --- a/test/test_eos_simple_macaw.cpp +++ b/test/test_eos_simple_macaw.cpp @@ -55,8 +55,8 @@ SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { rho += rho + static_cast(i); // cylce through a variety of densities const Real e = eos.SieColdCurve(rho); INFO("i: " << i << " rho = " << rho << " e = " << e); - REQUIRE_THAT(eos.TemperatureFromDensityInternalEnergy(rho, e), - Catch::Matchers::WithinRel(0.0, 1.0e-12)); + CHECK_THAT(eos.TemperatureFromDensityInternalEnergy(rho, e), + Catch::Matchers::WithinRel(0.0, 1.0e-12)); } // for } // Then } // When From 213557b2815307e69d92213dc9c50b7426195067 Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Thu, 25 Sep 2025 13:27:07 -0600 Subject: [PATCH 21/25] Return cold curve quantities when near cold curve --- singularity-eos/eos/eos_simple_macaw.hpp | 29 ++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index cf4a7c9fb59..d63542706f4 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -56,8 +56,15 @@ class SimpleMACAW : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (18) */ const Real Delta_e = std::max(sie - SieColdCurve(rho), 0.); + + // Early return ensure we're on the cold curve + if (_IsNearOrBelowZero(Delta_e)) { + return 0.; + } + const Real discriminant = Delta_e * (Delta_e + 4.0 * Cvinf_ * T0_ * std::pow(rho * v0_, Gc_)); + return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * Cvinf_); } PORTABLE_INLINE_FUNCTION void CheckParams() const { @@ -96,6 +103,12 @@ class SimpleMACAW : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const { /* Equation (19) */ const Real Delta_e = std::max(sie - SieColdCurve(rho), 0.); + + // Early return ensure we're on the cold curve + if (_IsNearOrBelowZero(Delta_e)) { + return PressureColdCurve(rho); + } + return PressureColdCurve(rho) + Gc_ * rho * Delta_e; } template @@ -103,6 +116,12 @@ class SimpleMACAW : public EosBase { const Real rho, const Real P, Indexer_t &&lambda = static_cast(nullptr)) const { const Real delta_p = std::max(P - PressureColdCurve(rho), 0.); + + // Early return ensure we're on the cold curve + if (_IsNearOrBelowZero(delta_p)) { + return SieColdCurve(rho); + } + return SieColdCurve(rho) + robust::ratio(delta_p, Gc_ * rho); } @@ -173,6 +192,12 @@ class SimpleMACAW : public EosBase { const Real Delta_e = std::max(sie - SieColdCurve(rho), 0.); /* Equation (21) (multiplied by $\rho$ to get bulk mod) */ const Real term1 = A_ * B_ * (B_ + 1.0) * std::pow(rho * v0_, B_ + 1.0); + + // Early return ensure we're on the cold curve + if (_IsNearOrBelowZero(Delta_e)) { + return term1; + } + const Real term2 = Gc_ * (Gc_ + 1.0) * rho * Delta_e; return term1 + term2; } @@ -319,6 +344,10 @@ class SimpleMACAW : public EosBase { static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; + PORTABLE_FORCEINLINE_FUNCTION bool _IsNearOrBelowZero(Real value) { + return value < std::numeric_limitsmin() * 5; + } + template PORTABLE_INLINE_FUNCTION Real _TemperatureScale( const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { From 5a7ae76a1913ad0389be73155be65ecf2dd1c4b6 Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Thu, 25 Sep 2025 13:27:30 -0600 Subject: [PATCH 22/25] Clang format --- singularity-eos/eos/eos_simple_macaw.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index d63542706f4..2cf5a2e9a9c 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -345,7 +345,7 @@ class SimpleMACAW : public EosBase { thermalqs::density | thermalqs::specific_internal_energy; PORTABLE_FORCEINLINE_FUNCTION bool _IsNearOrBelowZero(Real value) { - return value < std::numeric_limitsmin() * 5; + return value < std::numeric_limits min() * 5; } template From 342b5fe2a2b28eb9feae3f54cc1a4ccd32f9c2ca Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Thu, 25 Sep 2025 13:46:09 -0600 Subject: [PATCH 23/25] Add const qualifiers to _IsNearOrBelowZero and input --- singularity-eos/eos/eos_simple_macaw.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 2cf5a2e9a9c..05c0fcba67c 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -344,7 +344,7 @@ class SimpleMACAW : public EosBase { static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; - PORTABLE_FORCEINLINE_FUNCTION bool _IsNearOrBelowZero(Real value) { + PORTABLE_FORCEINLINE_FUNCTION bool _IsNearOrBelowZero(const Real value) const { return value < std::numeric_limits min() * 5; } From 0631d2e28f92c8fa61aa13fa17fd46f127de0ae0 Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Thu, 25 Sep 2025 13:59:38 -0600 Subject: [PATCH 24/25] Silly mistakes --- singularity-eos/eos/eos_simple_macaw.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 05c0fcba67c..54eeddcb52f 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -345,7 +345,7 @@ class SimpleMACAW : public EosBase { thermalqs::density | thermalqs::specific_internal_energy; PORTABLE_FORCEINLINE_FUNCTION bool _IsNearOrBelowZero(const Real value) const { - return value < std::numeric_limits min() * 5; + return (value <= std::numeric_limits::min() * 5); } template From 98c83ce746236d667e89c2cdcf05032fedd82b70 Mon Sep 17 00:00:00 2001 From: Jeff Peterson Date: Thu, 25 Sep 2025 17:00:29 -0600 Subject: [PATCH 25/25] Try epsilon rather than min --- singularity-eos/eos/eos_simple_macaw.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp index 54eeddcb52f..bffb692ebeb 100644 --- a/singularity-eos/eos/eos_simple_macaw.hpp +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -345,7 +345,7 @@ class SimpleMACAW : public EosBase { thermalqs::density | thermalqs::specific_internal_energy; PORTABLE_FORCEINLINE_FUNCTION bool _IsNearOrBelowZero(const Real value) const { - return (value <= std::numeric_limits::min() * 5); + return (value <= std::numeric_limits::epsilon() * 2); } template