From bb4b7972fae127374284579376540f13f32a3d89 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 29 Nov 2024 14:39:40 -0700 Subject: [PATCH 01/18] start adding mean atomic mass and number --- singularity-eos/eos/eos_base.hpp | 33 ++++++++++++++++++- singularity-eos/eos/eos_variant.hpp | 30 +++++++++++++++++ .../eos/modifiers/eos_unitsystem.hpp | 19 +++++++++++ singularity-eos/eos/modifiers/ramps_eos.hpp | 18 ++++++++++ .../eos/modifiers/relativistic_eos.hpp | 17 ++++++++++ singularity-eos/eos/modifiers/scaled_eos.hpp | 17 ++++++++++ singularity-eos/eos/modifiers/shifted_eos.hpp | 17 ++++++++++ 7 files changed, 150 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 83d0bf03a9..e9df3dbe7b 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -82,7 +82,9 @@ char *StrCat(char *destination, const char *source) { using EosBase::EntropyFromDensityTemperature; \ using EosBase::EntropyFromDensityInternalEnergy; \ using EosBase::GibbsFreeEnergyFromDensityTemperature; \ - using EosBase::GibbsFreeEnergyFromDensityInternalEnergy; + using EosBase::GibbsFreeEnergyFromDensityInternalEnergy; \ + using EosBase::MeanAtomicMass; \ + using EosBase::MeanAtomicNumber; // This macro adds several methods that most modifiers will // want. Not ALL modifiers will want these methods as written here, @@ -739,6 +741,35 @@ class EosBase { PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { return 0.0; } + // Defaults. Likely not accurate. You should almost always + // overwrite. For a given atom, always an integer, but this is a + // mean. + // JMM: EOS's which encapsulate a mix or reactions may wish to vary + // this. For example, Helmholtz and StellarCollapse. This isn't the + // default, so by default the base class provides a specialization. + // for models where density and temperature are required, the EOS + // developer is in charge of either throwing an error or choosing + // reasonable defaults. + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { return 1.0; } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { return 1.0; } + // TODO(JMM): Should we provide vector implementations if we depend on rho, T, etc? + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + CRTP copy = *(static_cast(this)); + return copy.MeanAtomicMass(rho, T, lambda); + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( + const Real rho, const Rela T, + Indexer_t &&lambda = static_cast(nullptr)) const { + CRTP copy = *(static_cast(this)); + return copy.MeanAtomicNumber(rho, T, lambda); + } + // Default entropy behavior is to cause an error PORTABLE_FORCEINLINE_FUNCTION void EntropyIsNotEnabled(const char *eosname) const { diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index 36711df58f..b89415dc2a 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -333,6 +333,36 @@ class Variant { return mpark::visit([](const auto &eos) { return eos.MinimumTemperature(); }, eos_); } + // Atomic mass/atomic number functions + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { + return mpark::visit([](const auto &eos) { return eos.MeanAtomicMass(); }, eos_); + } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { + return mpark::visit([](const auto &eos) { return eos.MeanAtomicNumber(); }, eos_); + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return mpark::visit( + [&rho, &T, &lambda](const auto &eos) { + return eos.MeanAtomicMassFromDensityTemperature(rho, T, lambda); + }, + eos_); + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( + const Real rho, const Rela T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return mpark::visit( + [&rho, &T, &lambda](const auto &eos) { + return eos.MeanAtomicNumberFromDensityTemperature(rho, T, lambda); + }, + eos_); + } + /* Vector versions of the member functions run on the host but the scalar lookups will run on the device diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 53d9435aab..06ac894e0f 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -249,6 +249,25 @@ class UnitSystem : public EosBase> { return inv_temp_unit_ * t_.MinimumTemperature(); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return t_.MeanAtomicMassFromDensityTemperature(rho * rho_unit_, T * temp_unit_, + lambda); + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return t_.MeanAtomicNumberFromDensityTemperature(rho * rho_unit_, T * temp_unit_, + lambda); + } + // vector implementations template inline void TemperatureFromDensityInternalEnergy( diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index a0ce08db2e..5db6f340ab 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -210,6 +210,7 @@ class BilinearRampEOS : public EosBase> { template PORTABLE_FUNCTION Real GruneisenParamFromDensityTemperature( const Real rho, const Real temperature, Indexer_t &&lambda = nullptr) const { + // TODO(JMM): This doesn't seem right. dpdrho relevant. return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda); } template @@ -242,6 +243,23 @@ class BilinearRampEOS : public EosBase> { return; } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return t_.MeanAtomicMassFromDensityTemperature(rho, T, lambda); + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return t_.MeanAtomicNumberFromDensityTemperature(rho, T, lambda); + } + // vector implementations template inline void TemperatureFromDensityInternalEnergy( diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index 0fd296d78d..72a34865f0 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -182,6 +182,23 @@ class RelativisticEOS : public EosBase> { t_.ValuesAtReferenceState(rho, temp, sie, press, cv, bmod, dpde, dvdt, lambda); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return t_.MeanAtomicMassFromDensityTemperature(rho, T, lambda); + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return t_.MeanAtomicNumberFromDensityTemperature(rho, T, lambda); + } + SG_ADD_MODIFIER_METHODS(T, t_); private: diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index bf64e90a04..76db90a01a 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -342,6 +342,23 @@ class ScaledEOS : public EosBase> { return t_.MinimumTemperature(); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return t_.MeanAtomicMassFromDensityTemperature(scale_ * rho, T, lambda); + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return t_.MeanAtomicNumberFromDensityTemperature(scale_ * rho, T, lambda); + } + SG_ADD_MODIFIER_METHODS(T, t_); private: diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index 79eaad6e4e..7a3e0762ff 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -352,6 +352,23 @@ class ShiftedEOS : public EosBase> { return t_.MinimumTemperature(); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return t_.MeanAtomicMassFromDensityTemperature(rho, T, lambda); + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( + const Real rho, const Rela T, + Indexer_t &&lambda = static_cast(nullptr)) const { + return t_.MeanAtomicNumberFromDensityTemperature(rho, T, lambda); + } + SG_ADD_MODIFIER_METHODS(T, t_); private: From 3e9a1f2920a84362fb9a939b9baf521e86dceffb Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sat, 30 Nov 2024 16:14:29 -0700 Subject: [PATCH 02/18] Rela -> Real --- singularity-eos/eos/eos_base.hpp | 2 +- singularity-eos/eos/eos_variant.hpp | 2 +- singularity-eos/eos/modifiers/shifted_eos.hpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index e9df3dbe7b..e107f2ee82 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -764,7 +764,7 @@ class EosBase { } template PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( - const Real rho, const Rela T, + const Real rho, const Real T, Indexer_t &&lambda = static_cast(nullptr)) const { CRTP copy = *(static_cast(this)); return copy.MeanAtomicNumber(rho, T, lambda); diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index b89415dc2a..31dcb7867a 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -354,7 +354,7 @@ class Variant { } template PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( - const Real rho, const Rela T, + const Real rho, const Real T, Indexer_t &&lambda = static_cast(nullptr)) const { return mpark::visit( [&rho, &T, &lambda](const auto &eos) { diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index 7a3e0762ff..465b6e585f 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -364,7 +364,7 @@ class ShiftedEOS : public EosBase> { } template PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( - const Real rho, const Rela T, + const Real rho, const Real T, Indexer_t &&lambda = static_cast(nullptr)) const { return t_.MeanAtomicNumberFromDensityTemperature(rho, T, lambda); } From f37996cb8bdb489e6929c8aad5c8262848036fa7 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Sat, 30 Nov 2024 16:23:08 -0700 Subject: [PATCH 03/18] dont shadow T parameter --- singularity-eos/eos/modifiers/eos_unitsystem.hpp | 12 ++++++------ singularity-eos/eos/modifiers/ramps_eos.hpp | 8 ++++---- singularity-eos/eos/modifiers/relativistic_eos.hpp | 8 ++++---- singularity-eos/eos/modifiers/scaled_eos.hpp | 8 ++++---- singularity-eos/eos/modifiers/shifted_eos.hpp | 8 ++++---- 5 files changed, 22 insertions(+), 22 deletions(-) diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 06ac894e0f..7fec533f65 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -255,17 +255,17 @@ class UnitSystem : public EosBase> { Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( - const Real rho, const Real T, + const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return t_.MeanAtomicMassFromDensityTemperature(rho * rho_unit_, T * temp_unit_, - lambda); + return t_.MeanAtomicMassFromDensityTemperature(rho * rho_unit_, + temperature * temp_unit_, lambda); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( - const Real rho, const Real T, + const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return t_.MeanAtomicNumberFromDensityTemperature(rho * rho_unit_, T * temp_unit_, - lambda); + return t_.MeanAtomicNumberFromDensityTemperature(rho * rho_unit_, + temperature * temp_unit_, lambda); } // vector implementations diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 5db6f340ab..202936186e 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -249,15 +249,15 @@ class BilinearRampEOS : public EosBase> { Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( - const Real rho, const Real T, + const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return t_.MeanAtomicMassFromDensityTemperature(rho, T, lambda); + return t_.MeanAtomicMassFromDensityTemperature(rho, temperature, lambda); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( - const Real rho, const Real T, + const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return t_.MeanAtomicNumberFromDensityTemperature(rho, T, lambda); + return t_.MeanAtomicNumberFromDensityTemperature(rho, temperature, lambda); } // vector implementations diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index 72a34865f0..ad2ba15e5a 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -188,15 +188,15 @@ class RelativisticEOS : public EosBase> { Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( - const Real rho, const Real T, + const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return t_.MeanAtomicMassFromDensityTemperature(rho, T, lambda); + return t_.MeanAtomicMassFromDensityTemperature(rho, temperature, lambda); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( - const Real rho, const Real T, + const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return t_.MeanAtomicNumberFromDensityTemperature(rho, T, lambda); + return t_.MeanAtomicNumberFromDensityTemperature(rho, temperature, lambda); } SG_ADD_MODIFIER_METHODS(T, t_); diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index 76db90a01a..d1ff50a2c0 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -348,15 +348,15 @@ class ScaledEOS : public EosBase> { Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( - const Real rho, const Real T, + const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return t_.MeanAtomicMassFromDensityTemperature(scale_ * rho, T, lambda); + return t_.MeanAtomicMassFromDensityTemperature(scale_ * rho, temperature, lambda); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( - const Real rho, const Real T, + const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return t_.MeanAtomicNumberFromDensityTemperature(scale_ * rho, T, lambda); + return t_.MeanAtomicNumberFromDensityTemperature(scale_ * rho, temperature, lambda); } SG_ADD_MODIFIER_METHODS(T, t_); diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index 465b6e585f..e006a0ce68 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -358,15 +358,15 @@ class ShiftedEOS : public EosBase> { Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( - const Real rho, const Real T, + const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return t_.MeanAtomicMassFromDensityTemperature(rho, T, lambda); + return t_.MeanAtomicMassFromDensityTemperature(rho, temperature, lambda); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( - const Real rho, const Real T, + const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const { - return t_.MeanAtomicNumberFromDensityTemperature(rho, T, lambda); + return t_.MeanAtomicNumberFromDensityTemperature(rho, temperature, lambda); } SG_ADD_MODIFIER_METHODS(T, t_); From 6b2c05c5a824aa67a2dd9f533ee58f176e5d832e Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 2 Dec 2024 10:39:32 -0700 Subject: [PATCH 04/18] Add MeanAtomicProperties for analytic EOS's to use and implement for ideal gas --- singularity-eos/eos/eos_base.hpp | 31 +++++++++++++++++++++++++++---- singularity-eos/eos/eos_ideal.hpp | 20 ++++++++++++++++---- 2 files changed, 43 insertions(+), 8 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index e107f2ee82..3f97babe09 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -134,6 +134,23 @@ struct Transform { Factor x, y, f; }; +/* + This is a utility struct used for analytic equations of state, + allowing them to store/report mean atomic mass/number easily. + */ +struct MeanAtomicProperties { + Real Abar, Zbar; + + // default is hydrogen + static constexpr Real DEFAULT_ABAR = 1.0; + static constexpr Real DEFAULT_ZBAR = 1.0; + + PORTABLE_INLINE_FUNCTION + MeanAtomicProperties(Real Abar_, Real Zbar_) : Abar(Abar_), Zbar(Zbar_) {} + PORTABLE_INLINE_FUNCTION + MeanAtomicProperties() : Abar(DEFAULT_ABAR), Zbar(DEFAULT_ZBAR) {} +}; + /* This is a CRTP that allows for static inheritance so that default behavior for various member functions can be defined. @@ -751,23 +768,29 @@ class EosBase { // developer is in charge of either throwing an error or choosing // reasonable defaults. PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { return 1.0; } + Real MeanAtomicMass() const { + PORTABLE_THROW_OR_ABORT("Mean atomic mass not implemented!"); + return 1.0; + } PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { return 1.0; } + Real MeanAtomicNumber() const { + PORTABLE_THROW_OR_ABORT("Mean atomic number not implemented!"); + return 1.0; + } // TODO(JMM): Should we provide vector implementations if we depend on rho, T, etc? template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( const Real rho, const Real T, Indexer_t &&lambda = static_cast(nullptr)) const { CRTP copy = *(static_cast(this)); - return copy.MeanAtomicMass(rho, T, lambda); + return copy.MeanAtomicMass(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( const Real rho, const Real T, Indexer_t &&lambda = static_cast(nullptr)) const { CRTP copy = *(static_cast(this)); - return copy.MeanAtomicNumber(rho, T, lambda); + return copy.MeanAtomicNumber(); } // Default entropy behavior is to cause an error diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp index d395928404..ddb17c243f 100644 --- a/singularity-eos/eos/eos_ideal.hpp +++ b/singularity-eos/eos/eos_ideal.hpp @@ -38,16 +38,20 @@ using namespace eos_base; class IdealGas : public EosBase { public: IdealGas() = default; - PORTABLE_INLINE_FUNCTION IdealGas(Real gm1, Real Cv) + PORTABLE_INLINE_FUNCTION + IdealGas(Real gm1, Real Cv, const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _Cv(Cv), _gm1(gm1), _rho0(_P0 / (_gm1 * _Cv * _T0)), _sie0(_Cv * _T0), _bmod0((_gm1 + 1) * _gm1 * _rho0 * _Cv * _T0), _dpde0(_gm1 * _rho0), - _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(_T0), _EntropyRho0(_rho0) { + _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(_T0), _EntropyRho0(_rho0), _AZbar(AZbar) { CheckParams(); } - PORTABLE_INLINE_FUNCTION IdealGas(Real gm1, Real Cv, Real EntropyT0, Real EntropyRho0) + PORTABLE_INLINE_FUNCTION + IdealGas(Real gm1, Real Cv, Real EntropyT0, Real EntropyRho0, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _Cv(Cv), _gm1(gm1), _rho0(_P0 / (_gm1 * _Cv * _T0)), _sie0(_Cv * _T0), _bmod0((_gm1 + 1) * _gm1 * _rho0 * _Cv * _T0), _dpde0(_gm1 * _rho0), - _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(EntropyT0), _EntropyRho0(EntropyRho0) { + _dvdt0(1. / (_rho0 * _T0)), _EntropyT0(EntropyT0), _EntropyRho0(EntropyRho0), + _AZbar(AZbar) { CheckParams(); } @@ -148,6 +152,12 @@ class IdealGas : 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; + + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { return _AZbar.Abar; } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { return _AZbar.Zbar; } + template PORTABLE_INLINE_FUNCTION void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, @@ -198,6 +208,8 @@ class IdealGas : public EosBase { thermalqs::density | thermalqs::specific_internal_energy; // optional entropy reference state variables Real _EntropyT0, _EntropyRho0; + // optional mean atomic mass and number + MeanAtomicProperties _AZbar; }; template From 2b7d0b5f316639c4f0889ca5e296f04b29bf90f1 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 2 Dec 2024 10:55:11 -0700 Subject: [PATCH 05/18] test --- test/test_eos_ideal.cpp | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp index 11b8b6447a..0ef1d9bc93 100644 --- a/test/test_eos_ideal.cpp +++ b/test/test_eos_ideal.cpp @@ -30,6 +30,7 @@ #include using singularity::IdealGas; +using singularity::MeanAtomicProperties; using EOS = singularity::Variant; SCENARIO("Ideal gas entropy", "[IdealGas][Entropy][GibbsFreeEnergy]") { @@ -80,6 +81,40 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy][GibbsFreeEnergy]") { } } +SCENARIO("Ideal gas mean atomic properties", "[IdealGas][MeanAtomicMass][MeanAtomicNumber]") { + constexpr Real Cv = 5.0; + constexpr Real gm1 = 0.4; + constexpr Real Abar = 4.0; // Helium + constexpr Real Zbar = 2.0; + const MeanAtomicProperties azbar(Abar, Zbar); + GIVEN("An ideal gas initialized with mean atomic poroperties") { + EOS host_eos = IdealGas(gm1, Cv, azbar); + WHEN("We evaluate it on host") { + Real Ab_eval = host_eos.MeanAtomicMass(); + Real Zb_eval = host_eos.MeanAtomicNumber(); + THEN("We get the right answer") { + REQUIRE( isClose(Ab_eval, Abar, 1e-12) ); + REQUIRE( isClose(Zb_eval, Zbar, 1e-12) ); + } + } + WHEN("We evaluate it on device, using a loop") { + constexpr int N = 100; + auto device_eos = host_eos.GetOnDevice(); + int nwrong = 0; + portableReduce("Check mean atomic number", 0, N, PORTABLE_LAMBDA(const int i, int &nw) { + double rho = i; + double T = 100.0*i; + double Ab_eval = device_eos.MeanAtomicMassFromDensityTemperature(rho, T); + double Zb_eval = device_eos.MeanAtomicNumberFromDensityTemperature(rho, T); + nw += !(isClose(Ab_eval, Abar, 1e-12)) + !(isClose(Zb_eval, Zbar, 1e-12)); + }, nwrong); + REQUIRE(nwrong == 0); + device_eos.Finalize(); + } + host_eos.Finalize(); + } +} + // A non-standard pattern where we do a reduction class CheckPofRE { public: From ade7cc33314015e4bed49330fad9cfea354aaa31 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 2 Dec 2024 13:58:26 -0700 Subject: [PATCH 06/18] formatting. Add carnahan starling AZbar --- singularity-eos/eos/eos_carnahan_starling.hpp | 19 +++++++++++---- test/test_eos_ideal.cpp | 24 +++++++++++-------- 2 files changed, 28 insertions(+), 15 deletions(-) diff --git a/singularity-eos/eos/eos_carnahan_starling.hpp b/singularity-eos/eos/eos_carnahan_starling.hpp index e77e844332..65ad4d79cf 100644 --- a/singularity-eos/eos/eos_carnahan_starling.hpp +++ b/singularity-eos/eos/eos_carnahan_starling.hpp @@ -38,7 +38,9 @@ using namespace eos_base; class CarnahanStarling : public EosBase { public: CarnahanStarling() = default; - PORTABLE_INLINE_FUNCTION CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq) + PORTABLE_INLINE_FUNCTION + CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(ROOM_TEMPERATURE), _P0(ATMOSPHERIC_PRESSURE), _qp(0.0), _rho0(DensityFromPressureTemperature(_P0, _T0)), _vol0(robust::ratio(1.0, _rho0)), @@ -46,18 +48,19 @@ class CarnahanStarling : public EosBase { _bmod0(_rho0 * Cv * _T0 * (PartialRhoZedFromDensity(_rho0) + ZedFromDensity(_rho0) * ZedFromDensity(_rho0) * gm1)), - _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1) { + _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1), _AZbar(AZbar) { CheckParams(); } - PORTABLE_INLINE_FUNCTION CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq, Real qp, - Real T0, Real P0) + PORTABLE_INLINE_FUNCTION + CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq, Real qp, Real T0, Real P0, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(T0), _P0(P0), _qp(qp), _rho0(DensityFromPressureTemperature(P0, T0)), _vol0(robust::ratio(1.0, _rho0)), _sie0(Cv * T0 + qq), _bmod0(_rho0 * Cv * T0 * (PartialRhoZedFromDensity(_rho0) + ZedFromDensity(_rho0) * ZedFromDensity(_rho0) * gm1)), - _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1) { + _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1), _AZbar(AZbar) { CheckParams(); } CarnahanStarling GetOnDevice() { return *this; } @@ -247,12 +250,18 @@ class CarnahanStarling : public EosBase { static std::string EosType() { return std::string("CarnahanStarling"); } static std::string EosPyType() { return EosType(); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { return _AZbar.Abar; } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { return _AZbar.Zbar; } + private: Real _Cv, _gm1, _bb, _qq; // optional reference state variables Real _T0, _P0, _qp; // reference values Real _rho0, _vol0, _sie0, _bmod0, _dpde0; + MeanAtomicProperties _AZbar; // static constexpr const Real _T0 = ROOM_TEMPERATURE; // static constexpr const Real _P0 = ATMOSPHERIC_PRESSURE; static constexpr const unsigned long _preferred_input = diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp index 0ef1d9bc93..dafc287b73 100644 --- a/test/test_eos_ideal.cpp +++ b/test/test_eos_ideal.cpp @@ -81,7 +81,8 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy][GibbsFreeEnergy]") { } } -SCENARIO("Ideal gas mean atomic properties", "[IdealGas][MeanAtomicMass][MeanAtomicNumber]") { +SCENARIO("Ideal gas mean atomic properties", + "[IdealGas][MeanAtomicMass][MeanAtomicNumber]") { constexpr Real Cv = 5.0; constexpr Real gm1 = 0.4; constexpr Real Abar = 4.0; // Helium @@ -93,21 +94,24 @@ SCENARIO("Ideal gas mean atomic properties", "[IdealGas][MeanAtomicMass][MeanAto Real Ab_eval = host_eos.MeanAtomicMass(); Real Zb_eval = host_eos.MeanAtomicNumber(); THEN("We get the right answer") { - REQUIRE( isClose(Ab_eval, Abar, 1e-12) ); - REQUIRE( isClose(Zb_eval, Zbar, 1e-12) ); + REQUIRE(isClose(Ab_eval, Abar, 1e-12)); + REQUIRE(isClose(Zb_eval, Zbar, 1e-12)); } } WHEN("We evaluate it on device, using a loop") { constexpr int N = 100; auto device_eos = host_eos.GetOnDevice(); int nwrong = 0; - portableReduce("Check mean atomic number", 0, N, PORTABLE_LAMBDA(const int i, int &nw) { - double rho = i; - double T = 100.0*i; - double Ab_eval = device_eos.MeanAtomicMassFromDensityTemperature(rho, T); - double Zb_eval = device_eos.MeanAtomicNumberFromDensityTemperature(rho, T); - nw += !(isClose(Ab_eval, Abar, 1e-12)) + !(isClose(Zb_eval, Zbar, 1e-12)); - }, nwrong); + portableReduce( + "Check mean atomic number", 0, N, + PORTABLE_LAMBDA(const int i, int &nw) { + double rho = i; + double T = 100.0 * i; + double Ab_eval = device_eos.MeanAtomicMassFromDensityTemperature(rho, T); + double Zb_eval = device_eos.MeanAtomicNumberFromDensityTemperature(rho, T); + nw += !(isClose(Ab_eval, Abar, 1e-12)) + !(isClose(Zb_eval, Zbar, 1e-12)); + }, + nwrong); REQUIRE(nwrong == 0); device_eos.Finalize(); } From 5edeed40aef04d2e6701ddb0e472448e34a49dd7 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 2 Dec 2024 14:06:48 -0700 Subject: [PATCH 07/18] eos davis --- singularity-eos/eos/eos_davis.hpp | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 2bfd96ad76..a9a6ac228a 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -35,9 +35,10 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION DavisReactants(const Real rho0, const Real e0, const Real P0, const Real T0, const Real A, const Real B, const Real C, const Real G0, const Real Z, - const Real alpha, const Real Cv0) + const Real alpha, const Real Cv0, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _rho0(rho0), _e0(e0), _P0(P0), _T0(T0), _A(A), _B(B), _C(C), _G0(G0), _Z(Z), - _alpha(alpha), _Cv0(Cv0) { + _alpha(alpha), _Cv0(Cv0), _AZbar(AZbar) { CheckParams(); } DavisReactants GetOnDevice() { return *this; } @@ -146,6 +147,10 @@ class DavisReactants : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const { return Gamma(rho); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { return _AZbar.Abar; } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { return _AZbar.Zbar; } template PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, @@ -183,6 +188,7 @@ class DavisReactants : public EosBase { private: static constexpr Real onethird = 1.0 / 3.0; Real _rho0, _e0, _P0, _T0, _A, _B, _C, _G0, _Z, _alpha, _Cv0; + MeanAtomicProperties _AZbar; // static constexpr const char _eos_type[] = "DavisReactants"; static constexpr unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; @@ -198,8 +204,9 @@ class DavisProducts : public EosBase { DavisProducts() = default; PORTABLE_INLINE_FUNCTION DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc, - const Real pc, const Real Cv) - : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv) {} + const Real pc, const Real Cv, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) + : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv), _AZbar(AZbar) {} PORTABLE_INLINE_FUNCTION void CheckParams() const { // TODO(JMM): Stub. @@ -285,6 +292,10 @@ class DavisProducts : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const { return Gamma(rho); } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { return _AZbar.Abar; } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { return _AZbar.Zbar; } template PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, @@ -321,6 +332,7 @@ class DavisProducts : public EosBase { private: static constexpr Real onethird = 1.0 / 3.0; Real _a, _b, _k, _n, _vc, _pc, _Cv; + MeanAtomicProperties _AZbar; // static constexpr const char _eos_type[] = "DavisProducts"; static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; From 388a596114b353bccd525fb4e3884d62605d6770 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 2 Dec 2024 18:59:20 -0700 Subject: [PATCH 08/18] Move default implementations of MeanAtomicMass/MeanAtomicNumber to macro --- singularity-eos/eos/eos_base.hpp | 50 +++++++++++-------- singularity-eos/eos/eos_carnahan_starling.hpp | 9 ++-- singularity-eos/eos/eos_davis.hpp | 16 +++--- singularity-eos/eos/eos_eospac.hpp | 8 +++ singularity-eos/eos/eos_gruneisen.hpp | 15 ++++-- singularity-eos/eos/eos_helmholtz.hpp | 25 ++++++++++ singularity-eos/eos/eos_ideal.hpp | 8 ++- singularity-eos/eos/eos_jwl.hpp | 14 ++++-- singularity-eos/eos/eos_mgusup.hpp | 10 +++- singularity-eos/eos/eos_noble_abel.hpp | 17 +++++-- singularity-eos/eos/eos_powermg.hpp | 9 +++- singularity-eos/eos/eos_sap_polynomial.hpp | 10 ++-- singularity-eos/eos/eos_spiner.hpp | 1 + singularity-eos/eos/eos_stellar_collapse.hpp | 26 ++++++++++ singularity-eos/eos/eos_stiff.hpp | 17 +++++-- singularity-eos/eos/eos_vinet.hpp | 12 +++-- .../eos/modifiers/eos_unitsystem.hpp | 4 -- singularity-eos/eos/modifiers/ramps_eos.hpp | 4 -- .../eos/modifiers/relativistic_eos.hpp | 4 -- singularity-eos/eos/modifiers/scaled_eos.hpp | 4 -- singularity-eos/eos/modifiers/shifted_eos.hpp | 4 -- 21 files changed, 181 insertions(+), 86 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 3f97babe09..f9422141ab 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -82,9 +82,17 @@ char *StrCat(char *destination, const char *source) { using EosBase::EntropyFromDensityTemperature; \ using EosBase::EntropyFromDensityInternalEnergy; \ using EosBase::GibbsFreeEnergyFromDensityTemperature; \ - using EosBase::GibbsFreeEnergyFromDensityInternalEnergy; \ - using EosBase::MeanAtomicMass; \ - using EosBase::MeanAtomicNumber; + using EosBase::GibbsFreeEnergyFromDensityInternalEnergy; + +// This macro adds these methods to a derived class. Due to scope, +// these can't be implemented in the base class, unless we make +// _AZbar public. Not all EOS's may want these default functions +// TODO(JMM): Should we go the alternate route and make _AZbar public? +#define SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) \ + PORTABLE_INLINE_FUNCTION \ + Real MeanAtomicMass() const { return _AZbar.Abar; } \ + PORTABLE_INLINE_FUNCTION \ + Real MeanAtomicNumber() const { return _AZbar.Zbar; } // This macro adds several methods that most modifiers will // want. Not ALL modifiers will want these methods as written here, @@ -103,7 +111,11 @@ char *StrCat(char *destination, const char *source) { } \ constexpr bool AllDynamicMemoryIsShareable() const { \ return t_.AllDynamicMemoryIsShareable(); \ - } + } \ + PORTABLE_INLINE_FUNCTION \ + Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } \ + PORTABLE_INLINE_FUNCTION \ + Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } class Factor { Real value_ = 1.0; @@ -135,8 +147,9 @@ struct Transform { }; /* - This is a utility struct used for analytic equations of state, - allowing them to store/report mean atomic mass/number easily. + This is a utility struct used to bundle mean atomic + mass/number. Used in the default implementations of MeanAtomicMass + and MeanAtomicNumber provided by the base class. */ struct MeanAtomicProperties { Real Abar, Zbar; @@ -149,6 +162,15 @@ struct MeanAtomicProperties { MeanAtomicProperties(Real Abar_, Real Zbar_) : Abar(Abar_), Zbar(Zbar_) {} PORTABLE_INLINE_FUNCTION MeanAtomicProperties() : Abar(DEFAULT_ABAR), Zbar(DEFAULT_ZBAR) {} + PORTABLE_INLINE_FUNCTION + void CheckParams() const { + PORTABLE_ALWAYS_REQUIRE(Abar > 0, "Positive mean atomic mass"); + PORTABLE_ALWAYS_REQUIRE(Zbar > 0, "Positive mean atomic number"); + } + void PrintParams() const { + printf(" Abar = %g\n", Abar); + printf(" Zbar = %g\n", Zbar); + } }; /* @@ -758,26 +780,14 @@ class EosBase { PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { return 0.0; } - // Defaults. Likely not accurate. You should almost always - // overwrite. For a given atom, always an integer, but this is a - // mean. // JMM: EOS's which encapsulate a mix or reactions may wish to vary // this. For example, Helmholtz and StellarCollapse. This isn't the // default, so by default the base class provides a specialization. // for models where density and temperature are required, the EOS // developer is in charge of either throwing an error or choosing // reasonable defaults. - PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { - PORTABLE_THROW_OR_ABORT("Mean atomic mass not implemented!"); - return 1.0; - } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { - PORTABLE_THROW_OR_ABORT("Mean atomic number not implemented!"); - return 1.0; - } - // TODO(JMM): Should we provide vector implementations if we depend on rho, T, etc? + // TODO(JMM): Should we provide vector implementations if we depend + // on rho, T, etc? template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( const Real rho, const Real T, diff --git a/singularity-eos/eos/eos_carnahan_starling.hpp b/singularity-eos/eos/eos_carnahan_starling.hpp index 65ad4d79cf..ba07356441 100644 --- a/singularity-eos/eos/eos_carnahan_starling.hpp +++ b/singularity-eos/eos/eos_carnahan_starling.hpp @@ -74,6 +74,7 @@ class CarnahanStarling : public EosBase { PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive"); PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive"); PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive"); + _AZbar.CheckParams(); } template PORTABLE_INLINE_FUNCTION Real ZedFromDensity( @@ -224,6 +225,8 @@ class CarnahanStarling : public EosBase { bmod = _bmod0; dpde = _dpde0; } + // Default accessors for mean atomic mass/number + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) // 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(CarnahanStarling) @@ -238,6 +241,7 @@ class CarnahanStarling : public EosBase { printf("Carnahan-Starling Parameters:\nGamma = %g\nCv = %g\nb = %g\nq = " "%g\n", _gm1 + 1.0, _Cv, _bb, _qq); + _AZbar.PrintParams(); } template PORTABLE_INLINE_FUNCTION void @@ -250,11 +254,6 @@ class CarnahanStarling : public EosBase { static std::string EosType() { return std::string("CarnahanStarling"); } static std::string EosPyType() { return EosType(); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { return _AZbar.Abar; } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { return _AZbar.Zbar; } - private: Real _Cv, _gm1, _bb, _qq; // optional reference state variables diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index a9a6ac228a..ff7091c969 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -46,6 +46,7 @@ class DavisReactants : public EosBase { void CheckParams() const { PORTABLE_REQUIRE(_rho0 >= 0, "Density must be positive"); PORTABLE_REQUIRE(_T0 >= 0, "Temperature must be positive"); + _AZbar.CheckParams(); } template PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( @@ -147,10 +148,6 @@ class DavisReactants : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const { return Gamma(rho); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { return _AZbar.Abar; } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { return _AZbar.Zbar; } template PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, @@ -165,6 +162,8 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const; + // Default accessors for mean atomic mass/number + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) // 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(DavisReactants) @@ -180,6 +179,7 @@ class DavisReactants : public EosBase { printf("%srho0:%e e0:%e P0:%e\nT0:%e A:%e B:%e\nC:%e G0:%e Z:%e\nalpha:%e " "Cv0:%e\n", s1, _rho0, _e0, _P0, _T0, _A, _B, _C, _G0, _Z, _alpha, _Cv0); + _AZbar.PrintParams(); } void inline Finalize() {} static std::string EosType() { return std::string("DavisReactants"); } @@ -209,6 +209,7 @@ class DavisProducts : public EosBase { : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv), _AZbar(AZbar) {} PORTABLE_INLINE_FUNCTION void CheckParams() const { + _AZbar.CheckParams(); // TODO(JMM): Stub. } DavisProducts GetOnDevice() { return *this; } @@ -292,10 +293,6 @@ class DavisProducts : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const { return Gamma(rho); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { return _AZbar.Abar; } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { return _AZbar.Zbar; } template PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, @@ -310,6 +307,8 @@ class DavisProducts : 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; + // Default accessors for mean atomic mass/number + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) // 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(DavisProducts) @@ -324,6 +323,7 @@ class DavisProducts : public EosBase { static constexpr char s1[]{"DavisProducts Params: "}; printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e \n", s1, _a, _b, _k, _n, _vc, _pc, _Cv); + _AZbar.PrintParams(); } inline void Finalize() {} static std::string EosType() { return std::string("DavisProducts"); } diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index 1f08a65873..3817ec06ee 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -141,6 +141,7 @@ class EOSPAC : public EosBase { // TODO(JMM): More validation checks? PORTABLE_ALWAYS_REQUIRE(rho_min_ >= 0, "Non-negative minimum density"); PORTABLE_ALWAYS_REQUIRE(temp_min_ >= 0, "Non-negative minimum temperature"); + AZbar_.CheckParams(); } inline EOSPAC GetOnDevice() { return *this; } @@ -457,6 +458,8 @@ class EOSPAC : public EosBase { using EosBase::FillEos; using EosBase::EntropyIsNotEnabled; + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_) + // EOSPAC vector implementations template inline void @@ -1097,6 +1100,7 @@ class EOSPAC : public EosBase { static std::string EosPyType() { return EosType(); } PORTABLE_INLINE_FUNCTION void PrintParams() const { printf("EOSPAC parameters:\nmatid = %i\n", matid_); + _AZbar.PrintParams(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rho_min_; } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return temp_min_; } @@ -1127,6 +1131,7 @@ class EOSPAC : public EosBase { // TODO(JMM): Is the fact that EOS_INTEGER isn't a size_t a // problem? Could it ever realistically overflow? EOS_INTEGER shared_size_, packed_size_; + MeanAtomicProperties AZbar_; static inline std::map &scratch_nbuffers() { static std::map nbuffers = { @@ -1194,6 +1199,9 @@ inline EOSPAC::EOSPAC(const int matid, bool invert_at_setup, Real insert_data, rho_ref_ = m.normalDensity; rho_min_ = m.rhoMin; temp_min_ = m.TMin; + // use std::max to hydrogen, in case of bad table + AZbar_.Abar = std::max(1.0, m.meanAtomicMass); + AZbar_.Zbar = std::max(1.0, m.meanAtomicNumber); EOS_REAL R[1] = {rho_ref_}; EOS_REAL T[1] = {temperatureToSesame(temp_ref_)}; diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index 93570a37e7..a2892e9a2a 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -39,9 +39,10 @@ class Gruneisen : public EosBase { PORTABLE_INLINE_FUNCTION Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0, const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv, - const Real rho_max) + const Real rho_max, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _C0(C0), _s1(s1), _s2(s2), _s3(s3), _G0(G0), _b(b), _rho0(rho0), _T0(T0), _P0(P0), - _Cv(Cv), _rho_max(rho_max) { + _Cv(Cv), _rho_max(rho_max), _AZbar(AZbar) { // Warn user when provided rho_max is greater than the computed rho_max #ifndef NDEBUG const Real computed_rho_max = ComputeRhoMax(s1, s2, s3, rho0); @@ -58,9 +59,11 @@ class Gruneisen : public EosBase { // Constructor when rho_max isn't specified automatically determines _rho_max PORTABLE_INLINE_FUNCTION Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0, - const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv) + const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _C0(C0), _s1(s1), _s2(s2), _s3(s3), _G0(G0), _b(b), _rho0(rho0), _T0(T0), _P0(P0), - _Cv(Cv), _rho_max(RHOMAX_SAFETY * ComputeRhoMax(s1, s2, s3, rho0)) {} + _Cv(Cv), _rho_max(RHOMAX_SAFETY * ComputeRhoMax(s1, s2, s3, rho0)), + _AZbar(AZbar) {} static PORTABLE_INLINE_FUNCTION Real ComputeRhoMax(const Real s1, const Real s2, const Real s3, const Real rho0); PORTABLE_INLINE_FUNCTION @@ -71,6 +74,7 @@ class Gruneisen : public EosBase { PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Non-negative heat capacity required"); PORTABLE_ALWAYS_REQUIRE(_rho_max > _rho0, "Maximum density must be greater than reference"); + _AZbar.CheckParams(); } PORTABLE_INLINE_FUNCTION Real MaxStableDensityAtTemperature(const Real temperature) const; @@ -150,6 +154,7 @@ class Gruneisen : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const; // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) SG_ADD_BASE_CLASS_USINGS(Gruneisen) PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } @@ -163,6 +168,7 @@ class Gruneisen : public EosBase { printf("%s C0:%e s1:%e s2:%e s3:%e\n G0:%e b:%e rho0:%e T0:%e\n P0:%eCv:%e " "rho_max:%e\n", s1, _C0, _s1, _s2, _s3, _G0, _b, _rho0, _T0, _P0, _Cv, _rho_max); + _AZbar.PrintParams(); } template PORTABLE_INLINE_FUNCTION void @@ -174,6 +180,7 @@ class Gruneisen : public EosBase { private: Real _C0, _s1, _s2, _s3, _G0, _b, _rho0, _T0, _P0, _Cv, _rho_max; + MeanAtomicProperties _AZbar; // static constexpr const char _eos_type[] = {"Gruneisen"}; PORTABLE_INLINE_FUNCTION Real Gamma(const Real rho_in) const { diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp index e5fbb0be78..b6ade5ac19 100644 --- a/singularity-eos/eos/eos_helmholtz.hpp +++ b/singularity-eos/eos/eos_helmholtz.hpp @@ -640,6 +640,31 @@ class Helmholtz : public EosBase { return gamma3 - 1.0; } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { + PORTABLE_THROW_OR_ABORT("For Helmholtz EOS, mean atomic mass is an input!\n"); + return 1.0; + } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { + PORTABLE_THROW_OR_ABORT("For Helmholtz EOS, mean atomic number is an input!\n"); + return 1.0; + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + using namespace HelmUtils; + return lambda[Lambda::Abar]; + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + using namespace HelmUtils; + return lambda[Lambda::Zbar]; + } + template PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp index ddb17c243f..02a0514674 100644 --- a/singularity-eos/eos/eos_ideal.hpp +++ b/singularity-eos/eos/eos_ideal.hpp @@ -71,6 +71,7 @@ class IdealGas : public EosBase { "Entropy reference temperature must be positive"); PORTABLE_ALWAYS_REQUIRE(_EntropyRho0 >= 0, "Entropy reference density must be positive"); + _AZbar.CheckParams(); return; } template @@ -153,11 +154,6 @@ class IdealGas : public EosBase { const unsigned long output, Indexer_t &&lambda = static_cast(nullptr)) const; - PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { return _AZbar.Abar; } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { return _AZbar.Zbar; } - template PORTABLE_INLINE_FUNCTION void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, @@ -175,6 +171,7 @@ class IdealGas : public EosBase { } // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) SG_ADD_BASE_CLASS_USINGS(IdealGas) PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } @@ -185,6 +182,7 @@ class IdealGas : public EosBase { static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { printf("Ideal Gas Parameters:\nGamma = %g\nCv = %g\n", _gm1 + 1.0, _Cv); + _AZbar.PrintParams(); } template PORTABLE_INLINE_FUNCTION void diff --git a/singularity-eos/eos/eos_jwl.hpp b/singularity-eos/eos/eos_jwl.hpp index 1931541630..31cf5477b6 100644 --- a/singularity-eos/eos/eos_jwl.hpp +++ b/singularity-eos/eos/eos_jwl.hpp @@ -40,17 +40,20 @@ class JWL : public EosBase { public: JWL() = default; PORTABLE_INLINE_FUNCTION JWL(const Real A, const Real B, const Real R1, const Real R2, - const Real w, const Real rho0, const Real Cv) + const Real w, const Real rho0, const Real Cv, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _A(A), _B(B), _R1(R1), _R2(R2), _w(w), _rho0(rho0), _Cv(Cv), - _c1(robust::ratio(A, rho0 * R1)), _c2(robust::ratio(B, rho0 * R2)) { - assert(R1 > 0.0); - assert(R2 > 0.0); + _c1(robust::ratio(A, rho0 * R1)), _c2(robust::ratio(B, rho0 * R2)), + _AZbar(AZbar) { + PORTABLE_REQUIRE(R1 > 0.0, "Ratio of rho0 to A > 0"); + PORTABLE_REQUIRE(R2 > 0.0, "Ratio of rho0 to B > 0"); CheckParams(); } PORTABLE_INLINE_FUNCTION void CheckParams() const { PORTABLE_ALWAYS_REQUIRE(_w > 0.0, "w > 0"); PORTABLE_ALWAYS_REQUIRE(_rho0 > 0.0, "Positive reference density"); PORTABLE_ALWAYS_REQUIRE(_Cv > 0.0, "Positive specific heat"); + _AZbar.CheckParams(); } JWL GetOnDevice() { return *this; } template @@ -111,6 +114,7 @@ class JWL : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const; // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) SG_ADD_BASE_CLASS_USINGS(JWL) PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } @@ -128,6 +132,7 @@ class JWL : public EosBase { static constexpr char s1[]{"JWL Params: "}; printf("%sA:%e B:%e R1: %e\nR2:%e w:%e rho0:%e\nCv:%e\n", s1, _A, _B, _R1, _R2, _w, _rho0, _Cv); + _AZbar.PrintParams(); } template PORTABLE_INLINE_FUNCTION void @@ -139,6 +144,7 @@ class JWL : public EosBase { private: Real _A, _B, _R1, _R2, _w, _rho0, _Cv, _c1, _c2; + MeanAtomicProperties _AZbar; PORTABLE_INLINE_FUNCTION Real ReferenceEnergy(const Real rho) const; PORTABLE_INLINE_FUNCTION Real ReferencePressure(const Real rho) const; // static constexpr const char _eos_type[] = "JWL"; diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index 409e596fe1..c315c1619a 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -38,8 +38,10 @@ class MGUsup : public EosBase { MGUsup() = default; // Constructor MGUsup(const Real rho0, const Real T0, const Real Cs, const Real s, const Real G0, - const Real Cv0, const Real E0, const Real S0) - : _rho0(rho0), _T0(T0), _Cs(Cs), _s(s), _G0(G0), _Cv0(Cv0), _E0(E0), _S0(S0) { + const Real Cv0, const Real E0, const Real S0, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) + : _rho0(rho0), _T0(T0), _Cs(Cs), _s(s), _G0(G0), _Cv0(Cv0), _E0(E0), _S0(S0), + _AZbar(AZbar) { CheckParams(); } PORTABLE_INLINE_FUNCTION void CheckParams() const; @@ -130,6 +132,7 @@ class MGUsup : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const; // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) SG_ADD_BASE_CLASS_USINGS(MGUsup) PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } @@ -142,6 +145,7 @@ class MGUsup : public EosBase { static constexpr char st[]{"MGUsup Params: "}; printf("%s rho0:%e T0:%e Cs:%e s:%e\n G0:%e Cv0:%e E0:%e S0:%e\n", st, _rho0, _T0, _Cs, _s, _G0, _Cv0, _E0, _S0); + _AZbar.PrintParams(); printf("\n\n"); } // Density/Energy from P/T not unique, if used will give error @@ -157,6 +161,7 @@ class MGUsup : public EosBase { static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; Real _rho0, _T0, _Cs, _s, _G0, _Cv0, _E0, _S0; + MeanAtomicProperties _AZbar; }; PORTABLE_INLINE_FUNCTION void MGUsup::CheckParams() const { @@ -178,6 +183,7 @@ PORTABLE_INLINE_FUNCTION void MGUsup::CheckParams() const { if (_Cv0 < 0.0) { PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter Cv0 < 0"); } + _AZbar.CheckParams(); } PORTABLE_INLINE_FUNCTION Real MGUsup::HugPressureFromDensity(Real rho) const { diff --git a/singularity-eos/eos/eos_noble_abel.hpp b/singularity-eos/eos/eos_noble_abel.hpp index a39816996b..e937d153b9 100644 --- a/singularity-eos/eos/eos_noble_abel.hpp +++ b/singularity-eos/eos/eos_noble_abel.hpp @@ -36,24 +36,27 @@ using namespace eos_base; class NobleAbel : public EosBase { public: NobleAbel() = default; - PORTABLE_INLINE_FUNCTION NobleAbel(Real gm1, Real Cv, Real bb, Real qq) + PORTABLE_INLINE_FUNCTION + NobleAbel(Real gm1, Real Cv, Real bb, Real qq, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(ROOM_TEMPERATURE), _P0(ATMOSPHERIC_PRESSURE), _qp(0.0), _rho0(robust::ratio(_P0, gm1 * Cv * _T0 + bb * _P0)), _vol0(robust::ratio(gm1 * Cv * _T0 + bb * _P0, _P0)), _sie0(Cv * _T0 + qq), _bmod0(robust::ratio(_rho0 * Cv * _T0 * gm1 * (gm1 + 1.0), (1.0 - bb * _rho0) * (1.0 - bb * _rho0))), - _dpde0(robust::ratio(_rho0 * gm1, 1.0 - bb * _rho0)) { + _dpde0(robust::ratio(_rho0 * gm1, 1.0 - bb * _rho0)), _AZbar(AZbar) { CheckParams(); } - PORTABLE_INLINE_FUNCTION NobleAbel(Real gm1, Real Cv, Real bb, Real qq, Real qp, - Real T0, Real P0) + PORTABLE_INLINE_FUNCTION + NobleAbel(Real gm1, Real Cv, Real bb, Real qq, Real qp, Real T0, Real P0, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(T0), _P0(P0), _qp(qp), _rho0(robust::ratio(P0, gm1 * Cv * T0 + bb * P0)), _vol0(robust::ratio(gm1 * Cv * T0 + bb * P0, P0)), _sie0(Cv * T0 + qq), _bmod0(robust::ratio(_rho0 * Cv * T0 * gm1 * (gm1 + 1.0), (1.0 - bb * _rho0) * (1.0 - bb * _rho0))), - _dpde0(robust::ratio(_rho0 * gm1, 1.0 - bb * _rho0)) { + _dpde0(robust::ratio(_rho0 * gm1, 1.0 - bb * _rho0)), _AZbar(AZbar) { CheckParams(); } NobleAbel GetOnDevice() { return *this; } @@ -67,6 +70,7 @@ class NobleAbel : public EosBase { PORTABLE_ALWAYS_REQUIRE(_Cv > 0, "Heat capacity must be positive"); PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive"); PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive"); + _AZbar.CheckParams(); } template PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( @@ -175,6 +179,7 @@ class NobleAbel : public EosBase { } // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) SG_ADD_BASE_CLASS_USINGS(NobleAbel) PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } @@ -187,6 +192,7 @@ class NobleAbel : public EosBase { printf("Noble-Abel Parameters:\nGamma = %g\nCv = %g\nb = %g\nq = " "%g\n", _gm1 + 1.0, _Cv, _bb, _qq); + _AZbar.PrintParams(); } template PORTABLE_INLINE_FUNCTION void @@ -206,6 +212,7 @@ class NobleAbel : public EosBase { Real _T0, _P0, _qp; // reference values Real _rho0, _vol0, _sie0, _bmod0, _dpde0; + MeanAtomicProperties _AZbar; // static constexpr const Real _T0 = ROOM_TEMPERATURE; // static constexpr const Real _P0 = ATMOSPHERIC_PRESSURE; static constexpr const unsigned long _preferred_input = diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp index ea8904285a..7e400b89a7 100644 --- a/singularity-eos/eos/eos_powermg.hpp +++ b/singularity-eos/eos/eos_powermg.hpp @@ -38,8 +38,10 @@ class PowerMG : public EosBase { PowerMG() = default; // Constructor PowerMG(const Real rho0, const Real T0, const Real G0, const Real Cv0, const Real E0, - const Real S0, const Real Pmin, const Real *expconsts) - : _rho0(rho0), _T0(T0), _G0(G0), _Cv0(Cv0), _E0(E0), _S0(S0), _Pmin(Pmin) { + const Real S0, const Real Pmin, const Real *expconsts, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) + : _rho0(rho0), _T0(T0), _G0(G0), _Cv0(Cv0), _E0(E0), _S0(S0), _Pmin(Pmin), + _AZbar(AZbar) { _InitializePowerMG(expconsts); CheckParams(); if (_Pmin >= 0.0) { @@ -136,6 +138,7 @@ class PowerMG : public EosBase { Indexer_t &&lambda = static_cast(nullptr)) const; // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) SG_ADD_BASE_CLASS_USINGS(PowerMG) PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } @@ -167,6 +170,7 @@ class PowerMG : public EosBase { static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; Real _rho0, _T0, _G0, _Cv0, _E0, _S0, _Pmin; + MeanAtomicProperties _AZbar; static constexpr const int PressureCoeffsK0toK40Size = 41; int _M; Real _K0toK40[PressureCoeffsK0toK40Size]; @@ -194,6 +198,7 @@ PORTABLE_INLINE_FUNCTION void PowerMG::CheckParams() const { if (_K0toK40[0] < 0.0) { PORTABLE_ALWAYS_THROW_OR_ABORT("Required PowerMG model parameter K0 < 0"); } + _AZbar.CheckParams(); } inline void PowerMG::_InitializePowerMG(const Real *K0toK40input) { diff --git a/singularity-eos/eos/eos_sap_polynomial.hpp b/singularity-eos/eos/eos_sap_polynomial.hpp index c26e77d429..7d361190c9 100644 --- a/singularity-eos/eos/eos_sap_polynomial.hpp +++ b/singularity-eos/eos/eos_sap_polynomial.hpp @@ -39,15 +39,17 @@ class SAP_Polynomial : public EosBase { PORTABLE_INLINE_FUNCTION SAP_Polynomial(const Real rho0, const Real a0, const Real a1, const Real a2c, const Real a2e, const Real a3, const Real b0, const Real b1, - const Real b2c, const Real b2e, const Real b3) + const Real b2c, const Real b2e, const Real b3, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _a0(a0), _a1(a1), _a2c(a2c), _a2e(a2e), _a3(a3), _b0(b0), _b1(b1), _b2c(b2c), - _b2e(b2e), _b3(b3), _rho0(rho0) { + _b2e(b2e), _b3(b3), _rho0(rho0), _AZbar(AZbar) { CheckParams(); } SAP_Polynomial GetOnDevice() { return *this; } PORTABLE_INLINE_FUNCTION void CheckParams() const { PORTABLE_ALWAYS_REQUIRE(_rho0 >= 0, "Reference density must be non-negative"); + _AZbar.CheckParams(); } template PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( @@ -183,6 +185,7 @@ class SAP_Polynomial : public EosBase { } // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) SG_ADD_BASE_CLASS_USINGS(SAP_Polynomial) PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } @@ -204,6 +207,7 @@ class SAP_Polynomial : public EosBase { printf(" b2c = %g\n", _b2c); printf(" b2e = %g\n", _b2e); printf(" b3 = %g\n", _b3); + _AZbar.PrintParams(); } template PORTABLE_INLINE_FUNCTION void @@ -221,7 +225,7 @@ class SAP_Polynomial : public EosBase { Real _a0, _a1, _a2c, _a2e, _a3, _b0, _b1, _b2c, _b2e, _b3; // reference values Real _rho0; - + MeanAtomicProperties _AZbar; static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; }; diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 6da8e47c4a..1a48967e23 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -307,6 +307,7 @@ class SpinerEOSDependsRhoT : public EosBase { Real rhoNormal_, TNormal_, sieNormal_, PNormal_; Real CvNormal_, bModNormal_, dPdENormal_, dVdTNormal_; Real lRhoOffset_, lTOffset_; // offsets must be non-negative + MeanAtomicProperties _AZbar; // TODO(JMM): Load from table. int matid_; bool reproducible_; // whereAmI_ and status_ used only for reporting. They are not thread-safe. diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 9901fffa8c..31c3821bf8 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -162,6 +162,32 @@ class StellarCollapse : public EosBase { Indexer_t &&lambda, Real &rho, Real &sie) const; // Properties of an NSE EOS + PORTABLE_INLINE_FUNCTION + Real MeanAtomicMass() const { + PORTABLE_THROW_OR_ABORT("For Helmholtz EOS, mean atomic mass is an input!\n"); + return 1.0; + } + PORTABLE_INLINE_FUNCTION + Real MeanAtomicNumber() const { + PORTABLE_THROW_OR_ABORT("For Helmholtz EOS, mean atomic number is an input!\n"); + return 1.0; + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + Real lRho, lT, Ye; + getLogsFromRhoT_(rho, temperature, lambda, lRho, lT, Ye); + return Abar_.interpToReal(Ye, lT, lRho); + } + template + PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature( + const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + Real lRho, lT, Ye; + getLogsFromRhoT_(rho, temperature, lambda, lRho, lT, Ye); + return Zbar_.interpToReal(Ye, lT, lRho); + } template PORTABLE_INLINE_FUNCTION void MassFractionsFromDensityTemperature( const Real rho, const Real temperature, Real &Xa, Real &Xh, Real &Xn, Real &Xp, diff --git a/singularity-eos/eos/eos_stiff.hpp b/singularity-eos/eos/eos_stiff.hpp index 4a64aeffc5..4c33ca77bc 100644 --- a/singularity-eos/eos/eos_stiff.hpp +++ b/singularity-eos/eos/eos_stiff.hpp @@ -36,24 +36,27 @@ using namespace eos_base; class StiffGas : public EosBase { public: StiffGas() = default; - PORTABLE_INLINE_FUNCTION StiffGas(Real gm1, Real Cv, Real Pinf, Real qq) + PORTABLE_INLINE_FUNCTION + StiffGas(Real gm1, Real Cv, Real Pinf, Real qq, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _Cv(Cv), _gm1(gm1), _Pinf(Pinf), _qq(qq), _T0(ROOM_TEMPERATURE), _P0(ATMOSPHERIC_PRESSURE), _qp(0.0), _rho0(robust::ratio(_P0 + Pinf, gm1 * Cv * _T0)), _vol0(robust::ratio(gm1 * Cv * _T0, _P0 + Pinf)), _sie0(robust::ratio(_P0 + (gm1 + 1.0) * Pinf, _P0 + Pinf) * Cv * _T0 + qq), _bmod0(gm1 * (gm1 + 1.0) * robust::ratio(_P0 + Pinf, gm1 * Cv * _T0) * Cv * _T0), - _dpde0(_rho0 * _gm1) { + _dpde0(_rho0 * _gm1), _AZbar(AZbar) { CheckParams(); } - PORTABLE_INLINE_FUNCTION StiffGas(Real gm1, Real Cv, Real Pinf, Real qq, Real qp, - Real T0, Real P0) + PORTABLE_INLINE_FUNCTION + StiffGas(Real gm1, Real Cv, Real Pinf, Real qq, Real qp, Real T0, Real P0, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) : _Cv(Cv), _gm1(gm1), _Pinf(Pinf), _qq(qq), _T0(T0), _P0(P0), _qp(qp), _rho0(robust::ratio(P0 + Pinf, gm1 * Cv * T0)), _vol0(robust::ratio(gm1 * Cv * _T0, _P0 + Pinf)), _sie0(robust::ratio(P0 + (gm1 + 1.0) * Pinf, P0 + Pinf) * Cv * T0 + qq), _bmod0(gm1 * (gm1 + 1.0) * robust::ratio(P0 + Pinf, gm1 * Cv * T0) * Cv * T0), - _dpde0(_rho0 * _gm1) { + _dpde0(_rho0 * _gm1), _AZbar(AZbar) { CheckParams(); } StiffGas GetOnDevice() { return *this; } @@ -66,6 +69,7 @@ class StiffGas : public EosBase { PORTABLE_INLINE_FUNCTION void CheckParams() const { PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive"); PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive"); + _AZbar.CheckParams(); } template PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( @@ -168,6 +172,7 @@ class StiffGas : public EosBase { // 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(StiffGas) + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } static constexpr unsigned long PreferredInput() { return _preferred_input; } @@ -179,6 +184,7 @@ class StiffGas : public EosBase { printf("Stiff Gas Parameters:\nGamma = %g\nCv = %g\nPinf = %g\nq = " "%g\n", _gm1 + 1.0, _Cv, _Pinf, _qq); + _AZbar.PrintParams(); } template PORTABLE_INLINE_FUNCTION void @@ -199,6 +205,7 @@ class StiffGas : public EosBase { Real _T0, _P0, _qp; // reference values Real _rho0, _vol0, _sie0, _bmod0, _dpde0; + MeanAtomicProperties _AZbar; // static constexpr const Real _T0 = ROOM_TEMPERATURE; // static constexpr const Real _P0 = ATMOSPHERIC_PRESSURE; static constexpr const unsigned long _preferred_input = diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp index 570d1bb319..8ea913e404 100644 --- a/singularity-eos/eos/eos_vinet.hpp +++ b/singularity-eos/eos/eos_vinet.hpp @@ -38,10 +38,13 @@ using namespace robust; class Vinet : public EosBase { public: Vinet() = default; - // Constructor + // TODO(JMM): Should Vinet and MGPower take AZbar before + // the series coefficients? Would require an overload... Vinet(const Real rho0, const Real T0, const Real B0, const Real BP0, const Real A0, - const Real Cv0, const Real E0, const Real S0, const Real *expconsts) - : _rho0(rho0), _T0(T0), _B0(B0), _BP0(BP0), _A0(A0), _Cv0(Cv0), _E0(E0), _S0(S0) { + const Real Cv0, const Real E0, const Real S0, const Real *expconsts, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) + : _rho0(rho0), _T0(T0), _B0(B0), _BP0(BP0), _A0(A0), _Cv0(Cv0), _E0(E0), _S0(S0), + _AZbar(AZbar) { CheckParams(); InitializeVinet(expconsts); } @@ -131,6 +134,7 @@ class Vinet : public EosBase { // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here SG_ADD_BASE_CLASS_USINGS(Vinet) + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } static constexpr unsigned long PreferredInput() { return _preferred_input; } @@ -161,6 +165,7 @@ class Vinet : public EosBase { static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::temperature; Real _rho0, _T0, _B0, _BP0, _A0, _Cv0, _E0, _S0; + MeanAtomicProperties _AZbar; static constexpr const int PressureCoeffsd2tod40Size = 39; static constexpr const int VinetInternalParametersSize = PressureCoeffsd2tod40Size + 4; Real _VIP[VinetInternalParametersSize], _d2tod40[PressureCoeffsd2tod40Size]; @@ -189,6 +194,7 @@ PORTABLE_INLINE_FUNCTION void Vinet::CheckParams() const { if (_Cv0 < 0.0) { PORTABLE_ALWAYS_THROW_OR_ABORT("Required Vinet model parameter Cv0 < 0"); } + _AZbar.CheckParams(); } inline void Vinet::InitializeVinet(const Real *d2tod40input) { diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 7fec533f65..957aa5eed6 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -249,10 +249,6 @@ class UnitSystem : public EosBase> { return inv_temp_unit_ * t_.MinimumTemperature(); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( const Real rho, const Real temperature, diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 202936186e..f316c57dec 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -243,10 +243,6 @@ class BilinearRampEOS : public EosBase> { return; } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( const Real rho, const Real temperature, diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index ad2ba15e5a..8729c30a83 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -182,10 +182,6 @@ class RelativisticEOS : public EosBase> { t_.ValuesAtReferenceState(rho, temp, sie, press, cv, bmod, dpde, dvdt, lambda); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( const Real rho, const Real temperature, diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index d1ff50a2c0..9576c306c4 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -342,10 +342,6 @@ class ScaledEOS : public EosBase> { return t_.MinimumTemperature(); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( const Real rho, const Real temperature, diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index e006a0ce68..c943bec517 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -352,10 +352,6 @@ class ShiftedEOS : public EosBase> { return t_.MinimumTemperature(); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicMass() const { return t_.MeanAtomicMass(); } - PORTABLE_INLINE_FUNCTION - Real MeanAtomicNumber() const { return t_.MeanAtomicNumber(); } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( const Real rho, const Real temperature, From 44b9fc8f1e6ff68d0c6b9241a12b83c73f36547f Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 2 Dec 2024 19:06:02 -0700 Subject: [PATCH 09/18] spiner already had Abar and Zbar in metadata. Just add it to class. --- singularity-eos/eos/eos_spiner.hpp | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 1a48967e23..b8f8ce0e19 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -80,6 +80,7 @@ class SpinerEOSDependsRhoT : public EosBase { struct Lambda { enum Index { lRho = 0, lT = 1 }; }; + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_) SG_ADD_BASE_CLASS_USINGS(SpinerEOSDependsRhoT); inline SpinerEOSDependsRhoT(const std::string &filename, int matid, bool reproduciblity_mode = false); @@ -307,7 +308,7 @@ class SpinerEOSDependsRhoT : public EosBase { Real rhoNormal_, TNormal_, sieNormal_, PNormal_; Real CvNormal_, bModNormal_, dPdENormal_, dVdTNormal_; Real lRhoOffset_, lTOffset_; // offsets must be non-negative - MeanAtomicProperties _AZbar; // TODO(JMM): Load from table. + MeanAtomicProperties AZbar_; // TODO(JMM): Load from table. int matid_; bool reproducible_; // whereAmI_ and status_ used only for reporting. They are not thread-safe. @@ -350,6 +351,7 @@ class SpinerEOSDependsRhoSie : public EosBase { }; using STricks = table_utils::SpinerTricks; + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_) SG_ADD_BASE_CLASS_USINGS(SpinerEOSDependsRhoSie); PORTABLE_INLINE_FUNCTION SpinerEOSDependsRhoSie() : memoryStatus_(DataStatus::Deallocated) {} @@ -539,6 +541,7 @@ class SpinerEOSDependsRhoSie : public EosBase { thermalqs::density | thermalqs::temperature; // static constexpr const char _eos_type[] = "SpinerEOSDependsRhoSie"; int matid_; + MeanAtomicProperties AZbar_; // TODO(JMM): Load from table. bool reproducible_; mutable RootFinding1D::Status status_; static constexpr const int _n_lambda = 1; @@ -717,6 +720,11 @@ inline herr_t SpinerEOSDependsRhoT::loadDataboxes_(const std::string &matid_str, status += H5LTget_attribute_double(file, matid_str.c_str(), SP5::Material::normalDensity, &rhoNormal_); rhoNormal_ = std::abs(rhoNormal_); + // Mean atomic mass and mean atomic number + status += H5LTget_attribute_double(file, matid_str.c_str(), + SP5::Material::meanAtomicMass, &(AZbar_.Abar)); + status += H5LTget_attribute_double(file, matid_str.c_str(), + SP5::Material::meanAtomicNumber, &(AZbar_.Zbar)); // tables status += P_.loadHDF(lTGroup, SP5::Fields::P); @@ -1553,6 +1561,11 @@ herr_t SpinerEOSDependsRhoSie::loadDataboxes_(const std::string &matid_str, hid_ status += H5LTget_attribute_double(file, matid_str.c_str(), SP5::Material::normalDensity, &rhoNormal_); rhoNormal_ = std::abs(rhoNormal_); + // Mean atomic mass and mean atomic number + status += H5LTget_attribute_double(file, matid_str.c_str(), + SP5::Material::meanAtomicMass, &(AZbar_.Abar)); + status += H5LTget_attribute_double(file, matid_str.c_str(), + SP5::Material::meanAtomicNumber, &(AZbar_.Zbar)); // sometimes independent variables status += sie_.loadHDF(lTGroup, SP5::Fields::sie); From 738b6f57d1f3c2ca8aff16d7b1fd9bfab6d06052 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 2 Dec 2024 19:16:21 -0700 Subject: [PATCH 10/18] temperature -> T --- singularity-eos/eos/eos_stellar_collapse.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 31c3821bf8..1a1e1ff0af 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -177,7 +177,7 @@ class StellarCollapse : public EosBase { const Real rho, const Real T, Indexer_t &&lambda = static_cast(nullptr)) const { Real lRho, lT, Ye; - getLogsFromRhoT_(rho, temperature, lambda, lRho, lT, Ye); + getLogsFromRhoT_(rho, T, lambda, lRho, lT, Ye); return Abar_.interpToReal(Ye, lT, lRho); } template @@ -185,7 +185,7 @@ class StellarCollapse : public EosBase { const Real rho, const Real T, Indexer_t &&lambda = static_cast(nullptr)) const { Real lRho, lT, Ye; - getLogsFromRhoT_(rho, temperature, lambda, lRho, lT, Ye); + getLogsFromRhoT_(rho, T, lambda, lRho, lT, Ye); return Zbar_.interpToReal(Ye, lT, lRho); } template From df746cc7e5013ead86fb4a262f170f657c6e6dd0 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 2 Dec 2024 19:32:50 -0700 Subject: [PATCH 11/18] CHANGELOG --- CHANGELOG.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e2e855636a..6f04083d0b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,7 +5,8 @@ ### Added (new features/APIs/variables/...) ### Fixed (Repair bugs, etc) -- [[OR437]](https://github.com/lanl/singularity-eos/pull/437) Fix segfault on HIP, clean up warnings, add strict sanitizer test +- [[PR439]](https://github.com/lanl/singularity-eos/pull/439) Add mean atomic mass and number to EOS API +- [[PR437]](https://github.com/lanl/singularity-eos/pull/437) Fix segfault on HIP, clean up warnings, add strict sanitizer test ### Changed (changing behavior/API/variables/...) From 25b959edc462ef77897569409e86e820b020ef6f Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Dec 2024 15:58:25 -0700 Subject: [PATCH 12/18] docs --- doc/sphinx/src/models.rst | 92 ++++++++++++++++++++++++++++++++++++ doc/sphinx/src/using-eos.rst | 50 ++++++++++++++++++++ 2 files changed, 142 insertions(+) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 1087807b2b..b8c8271d8b 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -338,6 +338,38 @@ Note: sometimes temperatures are measured in eV for which the conversion is Sesame units are equivalent to the mm-mg-µs unit system. +The ``MeanAtomicProperties`` struct +------------------------------------ + +Several analytic equations of state optionally accept mean atomic mass +and number as physics properties. These are the average number of +nucleons and protons in a constituent nucleus respectively. They are +not necessarily integers, as a given material may be made up of +multiple kinds of atom. For example, dry air contains both nitrogen +and oxygen. + +The mean atomic mass and number are frequently carried in the +container struct + +.. code-block:: cpp + + struct MeanAtomicProperties { + Real Abar, Zbar; + + // default is hydrogen + static constexpr Real DEFAULT_ABAR = 1.0; + static constexpr Real DEFAULT_ZBAR = 1.0; + + PORTABLE_INLINE_FUNCTION + MeanAtomicProperties(Real Abar_, Real Zbar_) : Abar(Abar_), Zbar(Zbar_) {} + PORTABLE_INLINE_FUNCTION + MeanAtomicProperties() : Abar(DEFAULT_ABAR), Zbar(DEFAULT_ZBAR) {} + }; + +which owns the atomic mass ``Abar`` and atomic number ``Zbar``. You +may set these by constructing the struct or by setting the fields in a +pre-constructed struct. The defaults are for hydrogen. + Implemented EOS models ---------------------- @@ -429,6 +461,14 @@ these values are not set, they will be the same as those returned by the conditions are given, the return values of the :code:`ValuesAtReferenceState()` function will not be the same. +Both constructors also optionally accept `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter, e.g., + +.. code-block:: cpp + + IdealGas(Real gm1, Real Cv, MeanAtomicProperties(Abar, Zbar)); + IdealGas(Real gm1, Real Cv, Real EntropyT0, Real EntropyRho0, MeanAtomicProperties(Abar, Zbar)); + Stiffened Gas ````````````` @@ -486,6 +526,15 @@ these values are not set, they will be the same as those returned by the conditions are given, the return values of the :code:`ValuesAtReferenceState()` function will not be the same. +Both constructors also optionally accept `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter, e.g., + +.. code-block:: cpp + + StiffGas(Real gm1, Real Cv, Real Pinf, Real q, MeanAtomicProperties(Abar, Zbar)); + StiffGas(Real gm1, Real Cv, Real Pinf, Real q, Real qp, Real T0, Real P0, + EntropyRho0, MeanAtomicProperties(Abar, Zbar)); + Noble-Abel `````````` @@ -544,6 +593,17 @@ these values are not set, they will be the same as those returned by the conditions are given, the return values of the :code:`ValuesAtReferenceState()` function will not be the same. +Both constructors also optionally accept `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter, e.g., + +.. code-block:: cpp + + NobleAbel(Real gm1, Real Cv, Real b, Real q, + MeanAtomicProperties(Abar, Zbar)); + NobleAbel(Real gm1, Real Cv, Real b, Real q, Real qp, Real T0, Real P0, + MeanAtomicProperties(Abar, Zbar)); + + Carnahan-Starling ````````````````` @@ -630,6 +690,16 @@ these values are not set, they will be the same as those returned by the conditions are given, the return values of the :code:`ValuesAtReferenceState()` function will not be the same. +Both constructors also optionally accept `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter, e.g., + +.. code-block:: cpp + + CarnahanStarling(Real gm1, Real Cv, Real b, Real q, + MeanAtomicProperties(Abar, Zbar)) + CarnahanStarling(Real gm1, Real Cv, Real b, Real q, Real qp, Real T0, Real P0, + MeanAtomicProperties(Abar, Zbar)) + Gruneisen EOS ````````````` @@ -767,6 +837,9 @@ There is an overload of the ``Gruneisen`` class which computes Gruneisen(const Real C0, const Real s1, const Real s2, const Real s3, const Real G0, const Real b, const Real rho0, const Real T0, const Real P0, const Real Cv) +Both constructors also optionally accept `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter. + Extendended Vinet EOS ````````````````````` @@ -849,6 +922,9 @@ is :math:`S_0`, and ``expconsts`` is a pointer to the constant array of length :math:`d_2` to :math:`d_{40}`. Expansion coefficients not used should be set to 0.0. +This constructor also optionally accepts `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter. + Mie-Gruneisen linear :math:`U_s`- :math:`u_p` EOS ````````````````````````````````````````````````` @@ -995,6 +1071,9 @@ where ``G0`` is :math:`\Gamma(\rho_0)`, ``Cv0`` is :math:`C_V`, ``E0`` is :math:`E_0`, and ``S0`` is :math:`S_0`. +This constructor also optionally accepts `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter. + Mie-Gruneisen power expansion EOS ````````````````````````````````` As we noted above, the assumption of a linear :math:`U_s`- :math:`u_p` relation is simply not valid at large compressions. At @@ -1075,6 +1154,8 @@ where :math:`K_0` to :math:`K_{40}`. Expansion coefficients not used should be set to :math:`0.0`. +This constructor also optionally accepts `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter. JWL EOS @@ -1127,6 +1208,9 @@ where ``A`` is :math:`A`, ``B`` is :math:`B`, ``R1`` is :math:`R_1`, ``R2`` is :math:`R_2`, ``w`` is :math:`w`, ``rho0`` is :math:`\rho_0`, and ``Cv`` is :math:`C_V`. +This constructor also optionally accepts `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter. + Davis EOS ``````````` @@ -1274,6 +1358,9 @@ where ``rho0`` is :math:`\rho_0`, ``e0`` is :math:`e_0`, ``P0`` is :math:`Z`, ``alpha`` is :math:`\alpha`, and ``Cv0`` is the specific heat capacity at the reference state. +This constructor also optionally accepts `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter. + Davis Products EOS ''''''''''''''''''' @@ -1381,6 +1468,9 @@ where ``a`` is :math:`a`, ``b`` is :math:`b`, ``k`` is :math:`k`, ``n`` is :math:`n`, ``vc`` is :math:`V_\mathrm{C}`, ``pc`` is :math:`P_\mathrm{C}`, ``Cv`` is :math:`C_{V,0}`. +This constructor also optionally accepts `MeanAtomicProperties` for +the atomic mass and number as a final optional parameter. + Spiner EOS ```````````` @@ -1630,6 +1720,8 @@ values in expansion and compression. and similar expressions for :math:`b_2^*`. +The SAP polynomial EOS also optionally accepts a +``MeanAtomicProperties`` struct. Stellar Collapse EOS diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index 8810914adb..e9a2ecf848 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1037,6 +1037,56 @@ specific internal energy in :math:`erg/g`. The function +.. code-block:: cpp + + template + void MeanAtomicMassFromDensityTemperature(const Real rho, const Real T, + Indexer_t &&lambda = nullpter) const; + +returns the mean atomic mass (i.e., the number of nucleons) of a +material given density in :math:`g/cm^3` and temperature in +Kelvin. The reason this is allowed to vary with density and +temperature is that some equations of state, such as the Stellar +Collapse and Helmholtz equations of state encapsulate reactive flows +where the average nucleus may depend on thermodynamic variables. For +most materials, however, this is not the case and a convenience +function that drops the dependence is available: + +.. code-block:: cpp + + Real MeanAtomicMass() const; + +The function + +.. code-block:: cpp + + template + void MeanAtomicNumberFromDensityTemperature(const Real rho, const Real T, + Indexer_t &&lambda = nullpter) const; + +returns the mean atomic number (i.e., the number of protons in the +nucleus) of a material given density in :math:`g/cm^3` and temperature +in Kelvin. The reason this is allowed to vary with density and +temperature is that some equations of state, such as the Stellar +Collapse and Helmholtz equations of state encapsulate reactive flows +where the average nucleus may depend on thermodynamic variables. For +most materials, however, this is not the case and a convenience +function that drops the dependence is available: + +.. code-block:: cpp + + Real MeanAtomicNumber() const; + + + +.. warning:: + + For materials where the mean atomic mass and number **do** vary with + density and temperature, the convenience call without this + dependence will produce an error. + +The function + .. code-block:: cpp template From 7829f318dfbadb1180b4d39999595dfeaec6eed3 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Dec 2024 16:05:41 -0700 Subject: [PATCH 13/18] docs for tables --- doc/sphinx/src/models.rst | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index b8c8271d8b..5f3b914856 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1533,6 +1533,12 @@ of the material in the file, and ``reproducability_mode`` is a boolean which slightly changes how initial guesses for root finds are computed. The constructor for ``SpinerEOSDependsRhoSie`` is identical. +.. note:: + + Mean atomic mass and number are loaded from input tables. The + ``SpinerEOS`` model does **not** support the + ``MeanAtomicProperties`` struct. + ``sp5`` files and ``sesame2spiner`` ````````````````````````````````````` @@ -1803,6 +1809,17 @@ where ``filename`` is the file containing the tabulated model, original `Stellar Collapse`_ format, and ``filter_bmod`` specifies whether or not to apply the above-described median filter. +.. note:: + + The ``StellarCollapse`` EOS assumes nuclear statistical equilibrium + and as such mean atomic mass and number are state variables. As such + class does not accept the ``MeanAtomicProperties`` struct. The + ``MeanAtomicMassFromDensityTemperature`` and + ``MeanAtomicNumberFromDensityTemperature`` functions return the + relevant quantities for some thermodynamic state. The + ``MeanAtomicMass()`` and ``MeanAtomicNumber()`` functions raise an + error. + ``StellarCollapse`` also provides .. cpp:function:: void Save(const std::string &filename) @@ -1888,6 +1905,14 @@ of the Helmholtz free energy (hence the name Helmholtz EOS). The free energy is pre-computed via integrals over the Fermi sphere and tabulated in a file provided from `Frank Timmes's website`_. +.. note:: + + Since mean atomic mass and number are required inputs, the + ``MeanAtomicMassFromDensityTemperature`` and + ``MeanAtomicNumberFromDensityAndTemperature`` functions simply + return the input values. The ``MeanAtomicMass()`` and + ``MeanAtomicNumber`` functions produce an error. + The table is a simple small ascii file. To ensure thermodyanic consistency, the table is interpolated using either biquintic or bicubic Hermite polynomials, which are sufficiently high order that @@ -1989,6 +2014,13 @@ and :math:`splitCowan` uses the cold curve plus Cowan-nuclear model for ions and the final option ``linear_interp`` uses linear instead of bilinear interpolation. +.. note:: + + Mean atomic mass and number are loaded from input tables. The + ``EOSPAC`` model does **not** support the ``MeanAtomicProperties`` + struct. + + Note for performance reasons this EOS uses a slightly different vector API. See :ref:`EOSPAC Vector Functions ` for more details. From 8e0c2fb9c3b85ac624e4b3fcfbf1f16a71907438 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Dec 2024 16:10:31 -0700 Subject: [PATCH 14/18] minor cleanup of comments/error messages --- singularity-eos/eos/eos_spiner.hpp | 4 ++-- singularity-eos/eos/eos_stellar_collapse.hpp | 6 ++++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index b8f8ce0e19..7a9bb2e285 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -308,7 +308,7 @@ class SpinerEOSDependsRhoT : public EosBase { Real rhoNormal_, TNormal_, sieNormal_, PNormal_; Real CvNormal_, bModNormal_, dPdENormal_, dVdTNormal_; Real lRhoOffset_, lTOffset_; // offsets must be non-negative - MeanAtomicProperties AZbar_; // TODO(JMM): Load from table. + MeanAtomicProperties AZbar_; int matid_; bool reproducible_; // whereAmI_ and status_ used only for reporting. They are not thread-safe. @@ -541,7 +541,7 @@ class SpinerEOSDependsRhoSie : public EosBase { thermalqs::density | thermalqs::temperature; // static constexpr const char _eos_type[] = "SpinerEOSDependsRhoSie"; int matid_; - MeanAtomicProperties AZbar_; // TODO(JMM): Load from table. + MeanAtomicProperties AZbar_; bool reproducible_; mutable RootFinding1D::Status status_; static constexpr const int _n_lambda = 1; diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 1a1e1ff0af..4ff7788dda 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -164,12 +164,14 @@ class StellarCollapse : public EosBase { // Properties of an NSE EOS PORTABLE_INLINE_FUNCTION Real MeanAtomicMass() const { - PORTABLE_THROW_OR_ABORT("For Helmholtz EOS, mean atomic mass is an input!\n"); + PORTABLE_THROW_OR_ABORT( + "For Stellar Collapse EOS, mean atomic mass is a state variable!\n"); return 1.0; } PORTABLE_INLINE_FUNCTION Real MeanAtomicNumber() const { - PORTABLE_THROW_OR_ABORT("For Helmholtz EOS, mean atomic number is an input!\n"); + PORTABLE_THROW_OR_ABORT( + "For Stellar Collapse EOS, mean atomic number is a state variable!\n"); return 1.0; } template From 2e64b7bbdd35196956205b4fa37b4ad3a277db04 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Dec 2024 16:12:26 -0700 Subject: [PATCH 15/18] add mean atomic mass/number test for tabulated --- test/test_eos_tabulated.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp index 9d0dcb72d4..ee01f2f04a 100644 --- a/test/test_eos_tabulated.cpp +++ b/test/test_eos_tabulated.cpp @@ -86,6 +86,11 @@ SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS][DependsRhoT][EOSPAC]") { REQUIRE(isClose(T, T_pac)); } + AND_THEN("We get the correct mean atomic mass and number") { + REQUIRE(steelEOS_host.MeanAtomicMass() == eospac.MeanAtomicMass()); + REQUIRE(steelEOS_host.MeanAtomicNumber() == eospac.MeanAtomicNumber()); + } + // TODO: this needs to be a much more rigorous test AND_THEN("Quantities can be read from density and temperature") { const Real sie_pac = eospac.InternalEnergyFromDensityTemperature(1e0, 1e6); @@ -233,6 +238,11 @@ SCENARIO("SpinerEOS depends on rho and sie", "[SpinerEOS][DependsRhoSie]") { THEN("The correct metadata is read in") { REQUIRE(steelEOS_host.matid() == steelID); + AND_THEN("We get the correct mean atomic mass and number") { + REQUIRE(steelEOS_host.MeanAtomicMass() == eospac.MeanAtomicMass()); + REQUIRE(steelEOS_host.MeanAtomicNumber() == eospac.MeanAtomicNumber()); + } + int nw_ie2{0}, nw_te2{0}; #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::fence(); From c397e6c763d719186af4bf28de931e517eb4c19d Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Dec 2024 16:19:25 -0700 Subject: [PATCH 16/18] AZbar_ in EOSPAC, not _AZbar --- singularity-eos/eos/eos_eospac.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index 3817ec06ee..763cb53695 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -1100,7 +1100,7 @@ class EOSPAC : public EosBase { static std::string EosPyType() { return EosType(); } PORTABLE_INLINE_FUNCTION void PrintParams() const { printf("EOSPAC parameters:\nmatid = %i\n", matid_); - _AZbar.PrintParams(); + AZbar_.PrintParams(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rho_min_; } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return temp_min_; } From 1e0f5631ce53f149447ac11b6482d83bb402464a Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Dec 2024 16:37:34 -0700 Subject: [PATCH 17/18] printparams should be host device --- singularity-eos/eos/eos_base.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index f9422141ab..8122925f2a 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -167,6 +167,7 @@ struct MeanAtomicProperties { PORTABLE_ALWAYS_REQUIRE(Abar > 0, "Positive mean atomic mass"); PORTABLE_ALWAYS_REQUIRE(Zbar > 0, "Positive mean atomic number"); } + PORTABLE_INLINE_FUNCTION void PrintParams() const { printf(" Abar = %g\n", Abar); printf(" Zbar = %g\n", Zbar); From 0084c806afd3edfb03e27cc8476007da388b62d8 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Dec 2024 16:47:02 -0700 Subject: [PATCH 18/18] Atomic mass and number --- test/test_eos_tabulated.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp index ee01f2f04a..6a31ed9451 100644 --- a/test/test_eos_tabulated.cpp +++ b/test/test_eos_tabulated.cpp @@ -239,8 +239,8 @@ SCENARIO("SpinerEOS depends on rho and sie", "[SpinerEOS][DependsRhoSie]") { REQUIRE(steelEOS_host.matid() == steelID); AND_THEN("We get the correct mean atomic mass and number") { - REQUIRE(steelEOS_host.MeanAtomicMass() == eospac.MeanAtomicMass()); - REQUIRE(steelEOS_host.MeanAtomicNumber() == eospac.MeanAtomicNumber()); + REQUIRE(isClose(steelEOS_host.MeanAtomicMass(), 55.37)); + REQUIRE(isClose(steelEOS_host.MeanAtomicNumber(), 25.80)); } int nw_ie2{0}, nw_te2{0};