diff --git a/CHANGELOG.md b/CHANGELOG.md index 215f57bf6c4..1ba7aa04258 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ ### Added (new features/APIs/variables/...) - [[PR556]](https://github.com/lanl/singularity-eos/pull/556) Add introspection into types available in the variant +- [[PR560]](https://github.com/lanl/singularity-eos/pull/560) Add Simple MACAW EOS - [[PR564]](https://github.com/lanl/singularity-eos/pull/564) Removed Get() function from IndexableTypes since it could have unexpected consequences when a type wasn't present ### Fixed (Repair bugs, etc) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 0bc2cf1350f..31792f09f35 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -21,6 +21,9 @@ .. _PowerMG: https://www.osti.gov/biblio/1762624 +.. _SimpleMACAW: https://www.osti.gov/biblio/2479474 + +.. _MACAW EOS: https://pubs.aip.org/aip/jap/article/134/12/125102/2912256 EOS Models =========== @@ -1320,6 +1323,56 @@ This constructor also optionally accepts `MeanAtomicProperties` for the atomic mass and number as a final optional parameter. +Simple MACAW +```````````` +The `Simple MACAW EOS <_SimpleMACAW>`_ is a simplified version of the `MACAW EOS `_ +and is thermodynamically consistent. It is constructed from a the Helmholtz +free energy using a Murnaghan cold curve and a simplified thermal component +from the MACAW EOS. + +Fundamentally, the equation of state can be written in Mie-Gruneisen form (with constant Gruneisen parameter) as: + +.. math:: + + P(v, e) = P_{\text{cold}}(v) + \Gamma_c \rho (e - e_{\text{cold}}(v)) + +where the cold curves are given by: + +.. math:: + + e_{\text{cold}}(v) = A v_0 \Big[ \Big( \frac{v}{v_0} \Big)^{-B} + \Big( \frac{v}{v_0} \Big) B - (B+1) \Big] + +and + +.. math:: + + p_{\text{cold}}(v) := -e'_{\text{cold}}(v) = AB \Big( \Big( \frac{v}{v_0} \Big)^{-(B+1)} - 1 \Big) + +The specific heat capacity at constant volume for this model is given by, + +.. math:: + + C_v(v, \tau) = C^\infty_v \frac{\tau^2 + 2\tau}{(\tau + 1)^2} \quad \text{where } \tau = \frac{T}{\theta(v)} \quad \text{ and } \quad \theta(v) := T_0 \Big( \frac{v}{v_c} \Big)^{-\Gamma_c} + +Note it obeys the expected physical behavior; that, :math:`\lim_{T\to 0^+} C_v(v,\tau(v,T)) = 0` and +:math:`\lim_{T\to\infty} C_v(v,\tau(v,T) = C^\infty_v < \infty` (Dulong-Petit law). + +The constructor for the Simple MACAW EOS is + +.. code-block:: cpp + + SimpleMACAW(const Real A, const Real B, const Real Cvinf, const Real v0, + const Real T0, const Real Gc) + +where ``A`` is :math:`A`, ``B`` is :math:`B`, ``Cvinf`` is :math:`C^\infty_v`, +``v0`` is :math:`v_0`, ``T0`` is :math:`T_0`, ``Gc`` is :math:`\Gamma_c`. + +In order to maintain thermodynamic stability, a sufficient set of constraints +is given by :math:`A > 0`, :math:`B > 0`, :math:`C^\infty_v > 0`, :math:`v_0 > +0`, :math:`T_0 > 0`, and :math:`\Gamma_c \in (0,1]`. One can still select +:math:`\Gamma_c > 1`, just note that the isothermal bulk modulus can be +negative (the isentropic bulk modulus will still be positive though). + JWL EOS `````````` diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 81302eb703c..51b71e5a4af 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -53,6 +53,7 @@ register_headers( eos/eos_gruneisen.hpp eos/eos_vinet.hpp eos/eos_builder.hpp + eos/eos_simple_macaw.hpp eos/eos_jwl.hpp eos/eos_helmholtz.hpp eos/eos_sap_polynomial.hpp diff --git a/singularity-eos/eos/eos_models.hpp b/singularity-eos/eos/eos_models.hpp index f0ce27612ca..497baf1f5bf 100644 --- a/singularity-eos/eos/eos_models.hpp +++ b/singularity-eos/eos/eos_models.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include diff --git a/singularity-eos/eos/eos_simple_macaw.hpp b/singularity-eos/eos/eos_simple_macaw.hpp new file mode 100644 index 00000000000..bffb692ebeb --- /dev/null +++ b/singularity-eos/eos/eos_simple_macaw.hpp @@ -0,0 +1,378 @@ +//------------------------------------------------------------------------------ +// © 2021-2025. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#ifndef _SINGULARITY_EOS_EOS_EOS_SIMPLE_MACAW_HPP_ +#define _SINGULARITY_EOS_EOS_EOS_SIMPLE_MACAW_HPP_ + +// stdlib +#include +#include +#include + +// Ports-of-call +#include + +// Base stuff +#include +#include +#include +#include + +namespace singularity { + +using namespace eos_base; + +/* Most information regarding this equation of state can be found here: + * https://www.osti.gov/biblio/2479474 + * + * Note that all equations are given as functions of the density ($\rho$) rather than + * functions of the specific volume ($v$). The reason for this is it simplifies the + * computational complexity. + */ +class SimpleMACAW : public EosBase { + public: + SimpleMACAW() = default; + PORTABLE_INLINE_FUNCTION + SimpleMACAW(Real A, Real B, Real Cvinf, Real v0, Real T0, Real Gc, + const MeanAtomicProperties &AZbar = MeanAtomicProperties()) + : A_(A), B_(B), Cvinf_(Cvinf), v0_(v0), T0_(T0), Gc_(Gc), _AZbar(AZbar) { + CheckParams(); + } + SimpleMACAW GetOnDevice() { return *this; } + template + PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (18) */ + const Real Delta_e = std::max(sie - SieColdCurve(rho), 0.); + + // Early return ensure we're on the cold curve + if (_IsNearOrBelowZero(Delta_e)) { + return 0.; + } + + const Real discriminant = + Delta_e * (Delta_e + 4.0 * Cvinf_ * T0_ * std::pow(rho * v0_, Gc_)); + + return robust::ratio((Delta_e + std::sqrt(discriminant)), 2.0 * Cvinf_); + } + PORTABLE_INLINE_FUNCTION void CheckParams() const { + PORTABLE_ALWAYS_REQUIRE(A_ >= 0, "Parameter, 'A', must be non-negative"); + PORTABLE_ALWAYS_REQUIRE(B_ >= 0, "Parameter, 'B', must be non-negative"); + PORTABLE_ALWAYS_REQUIRE( + Cvinf_ > 0, + "Specific heat capacity (at constant volume), `Cvinf`, must be positive"); + PORTABLE_ALWAYS_REQUIRE(v0_ > 0, "Reference specific volume, 'v0', must be positive"); + PORTABLE_ALWAYS_REQUIRE(T0_ >= 0, + "Reference temperature, 'T0', must be non-negative"); + PORTABLE_ALWAYS_REQUIRE(Gc_ > 0, "Gruneisen parameter, 'Gc', must be positive"); + if (Gc_ > 1) { + PORTABLE_WARN("Warning: Gruneisen coefficient greater than 1. Thermodynamic " + "stability may be violated."); + } + _AZbar.CheckParams(); + } + template + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (16) */ + return SieColdCurve(rho) + _SieThermalPortion(rho, temperature); + } + template + PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (15) */ + return PressureColdCurve(rho) + _PressureThermalPortion(rho, temperature); + } + template + PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (19) */ + const Real Delta_e = std::max(sie - SieColdCurve(rho), 0.); + + // Early return ensure we're on the cold curve + if (_IsNearOrBelowZero(Delta_e)) { + return PressureColdCurve(rho); + } + + return PressureColdCurve(rho) + Gc_ * rho * Delta_e; + } + template + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityPressure( + const Real rho, const Real P, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real delta_p = std::max(P - PressureColdCurve(rho), 0.); + + // Early return ensure we're on the cold curve + if (_IsNearOrBelowZero(delta_p)) { + return SieColdCurve(rho); + } + + return SieColdCurve(rho) + robust::ratio(delta_p, Gc_ * rho); + } + + // Entropy member functions + template + PORTABLE_INLINE_FUNCTION Real + EntropyFromDensityTemperature(const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (8) */ + const Real tau = robust::ratio(temperature, _TemperatureScale(rho)); + return Cvinf_ * (robust::ratio(tau, 1.0 + tau) + std::log(1.0 + tau)); + } + template + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real T = TemperatureFromDensityInternalEnergy(rho, sie); + return EntropyFromDensityTemperature(rho, T); + } + + // Cold curve, thermal portion, and other helper functions + + /* Note: The cold curves for the simple MACAW are presented as functions of `v` rather + * than `rho`. We keep them as functions of `rho` for computational efficiency. */ + template + PORTABLE_INLINE_FUNCTION Real + SieColdCurve(const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { + const Real ratio = rho * v0_; + return A_ * v0_ * (std::pow(ratio, B_) + robust::ratio(B_, ratio) - (B_ + 1.)); + } + template + PORTABLE_INLINE_FUNCTION Real PressureColdCurve( + const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { + const Real ratio = rho * v0_; + return A_ * B_ * (std::pow(ratio, B_ + 1.0) - 1.0); + } + + // Specific heat capacity at constant volume + template + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (6) */ + const Real tau = robust::ratio(temperature, _TemperatureScale(rho)); + const Real numerator = math_utils::pow<2>(tau) + 2.0 * tau; + const Real denominator = math_utils::pow<2>(tau + 1.0); + return Cvinf_ * robust::ratio(numerator, denominator); + } + template + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real T = TemperatureFromDensityInternalEnergy(rho, sie); + return SpecificHeatFromDensityTemperature(rho, T); + } + // Isentropic bulk modulus + template + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real sie = InternalEnergyFromDensityTemperature(rho, temperature); + return BulkModulusFromDensityInternalEnergy(rho, sie); + } + template + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real Delta_e = std::max(sie - SieColdCurve(rho), 0.); + /* Equation (21) (multiplied by $\rho$ to get bulk mod) */ + const Real term1 = A_ * B_ * (B_ + 1.0) * std::pow(rho * v0_, B_ + 1.0); + + // Early return ensure we're on the cold curve + if (_IsNearOrBelowZero(Delta_e)) { + return term1; + } + + const Real term2 = Gc_ * (Gc_ + 1.0) * rho * Delta_e; + return term1 + term2; + } + // Isothermal bulk modulus + template + PORTABLE_INLINE_FUNCTION Real IsothermalBulkModulusFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (22) (multiplied by $\rho$ to get bulk mod) + * As mentioned in the paper, from the constraints on the parameters, + * the isothermal bulk modulus is guaranteed to be positive. */ + const Real term1 = A_ * B_ * (B_ + 1.0) * std::pow(rho * v0_, B_ + 1.0); + // There is a mistake in the paper and the numerator term should have a -Gc_ instead + // of -2*Gc_ + const Real numerator = + rho * (T0_ * (1.0 - Gc_) * std::pow(rho * v0_, Gc_) + temperature); + const Real denominator = T0_ * std::pow(rho * v0_, Gc_) + temperature; + return term1 + Cvinf_ * Gc_ * temperature * temperature * + robust::ratio(numerator, denominator * denominator); + } + // Specific heat capacity at constant pressure + template + PORTABLE_INLINE_FUNCTION Real ConstantPressureSpecificHeatFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real BT = IsothermalBulkModulusFromDensityTemperature(rho, temperature); + const Real Bs = BulkModulusFromDensityTemperature(rho, temperature); + const Real cv = SpecificHeatFromDensityTemperature(rho, temperature); + return robust::ratio(Bs * cv, BT); /* General thermodynamic identity */ + } + // Coefficient of thermal expansivity + template + PORTABLE_INLINE_FUNCTION Real CoefficientThermalExpansionFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real BT = IsothermalBulkModulusFromDensityTemperature(rho, temperature); + const Real cv = SpecificHeatFromDensityTemperature(rho, temperature); + return robust::ratio(Gc_ * rho * cv, BT); /* General thermodynamic identity */ + } + + // Gruneisen parameter + template + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( + const Real rho, const Real temperature, + Indexer_t &&lambda = static_cast(nullptr)) const { + return Gc_; + } + template + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy( + const Real rho, const Real sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + return Gc_; + } + template + PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &sie, Real &press, + Real &cv, Real &bmod, const unsigned long output, + Indexer_t &&lambda) const { + if (output & thermalqs::density && output & thermalqs::specific_internal_energy) { + PORTABLE_ALWAYS_REQUIRE( + !(output & thermalqs::pressure || output & thermalqs::temperature), + "Cannot request more than two thermodynamic variables (rho, T, e, P)"); + DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie); + } + if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) { + PORTABLE_ALWAYS_REQUIRE( + !(output & thermalqs::density || output & thermalqs::temperature), + "Cannot request more than two thermodynamic variables (rho, T, e, P)"); + sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); + } + if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) { + sie = InternalEnergyFromDensityPressure(rho, press); + } + if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::temperature) + temp = TemperatureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::bulk_modulus) + bmod = BulkModulusFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::specific_heat) + cv = SpecificHeatFromDensityInternalEnergy(rho, sie); + } + template + PORTABLE_INLINE_FUNCTION void + ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, + Real &bmod, Real &dpde, Real &dvdt, + Indexer_t &&lambda = static_cast(nullptr)) const { + // v_0 and T_0 are the only user selected reference state parameters + rho = 1.0 / v0_; + temp = T0_; + press = PressureFromDensityTemperature(rho, temp); + sie = InternalEnergyFromDensityTemperature(rho, temp); + + cv = SpecificHeatFromDensityInternalEnergy(rho, sie); + bmod = BulkModulusFromDensityInternalEnergy(rho, sie); + dpde = rho * GruneisenParamFromDensityTemperature(rho, temp); + } + // Generic functions provided by the base class. These contain e.g. the vector + // overloads that use the scalar versions declared here + SG_ADD_BASE_CLASS_USINGS(SimpleMACAW) + SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(_AZbar) + + static constexpr unsigned long PreferredInput() { return _preferred_input; } + PORTABLE_INLINE_FUNCTION void PrintParams() const { + printf("Simple MACAW Parameters:\nA = %g\nB = %g\nCvinf = %g\nv0 = %g\nT0 = " + "%g\nGc = " + "%g\n", + A_, B_, Cvinf_, v0_, T0_, Gc_); + _AZbar.PrintParams(); + } + template + PORTABLE_INLINE_FUNCTION void + DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Indexer_t &&lambda, Real &rho, Real &sie) const { + /* Since the isothermal bulk modulus is always positive (assuming parameter + * constraints), this guarantees the function to solve is monotonic hence we always + * have a unique solution. */ + + // Setup lambda function for rho + auto f = [=](const Real x /* density */) { + const Real term1 = A_ * B_ * (std::pow(v0_ * x, B_ + 1.0) - 1.0) - press; + const Real term2 = T0_ * std::pow(v0_ * x, Gc_) + temp; + const Real term3 = Cvinf_ * Gc_ * temp * temp * x; + return term1 * term2 + term3; + }; + + const RootFinding1D::RootCounts root_info; + + // Finding the density on pressure cold curve is a guaranteed upper bound + const Real rho_high = std::pow(press / (A_ * B_) + 1.0, 1.0 / (B_ + 1.0)) / v0_; + const Real rho_low = 0.0; // zero density is valid for `f` defined above. + + regula_falsi(f, 0.0 /*target*/, 1.0 / v0_ /*guess*/, rho_low /*left bracket*/, + rho_high /*right bracket*/, 1.0e-9 /*x? tol*/, 1.0e-9 /*y? tol*/, rho, + &root_info, true); + + sie = InternalEnergyFromDensityTemperature(rho, temp); + } + inline void Finalize() {} + static std::string EosType() { return std::string("SimpleMACAW"); } + static std::string EosPyType() { return EosType(); } + + private: + Real A_, B_, Cvinf_, v0_, T0_, Gc_; + MeanAtomicProperties _AZbar; + static constexpr const unsigned long _preferred_input = + thermalqs::density | thermalqs::specific_internal_energy; + + PORTABLE_FORCEINLINE_FUNCTION bool _IsNearOrBelowZero(const Real value) const { + return (value <= std::numeric_limits::epsilon() * 2); + } + + template + PORTABLE_INLINE_FUNCTION Real _TemperatureScale( + const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { + /* Equation (13) */ + return T0_ * std::pow(rho * v0_, Gc_); + } + + template + PORTABLE_INLINE_FUNCTION Real + _SieThermalPortion(const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real ratio = rho * v0_; + return Cvinf_ * T * (1.0 - robust::ratio(T0_, T0_ + T * std::pow(ratio, -Gc_))); + } + + template + PORTABLE_INLINE_FUNCTION Real + _PressureThermalPortion(const Real rho, const Real T, + Indexer_t &&lambda = static_cast(nullptr)) const { + const Real numerator = rho * Cvinf_ * Gc_ * math_utils::pow<2>(T); + const Real denominator = T0_ * std::pow(rho * v0_, Gc_) + T; + return robust::ratio(numerator, denominator); + } +}; + +} // namespace singularity + +#endif // _SINGULARITY_EOS_EOS_EOS_SIMPLE_MACAW_HPP_ diff --git a/singularity-eos/eos/eos_type_lists.hpp b/singularity-eos/eos/eos_type_lists.hpp index 30aaced18fa..8895691bd78 100644 --- a/singularity-eos/eos/eos_type_lists.hpp +++ b/singularity-eos/eos/eos_type_lists.hpp @@ -49,7 +49,8 @@ using singularity::variadic_utils::transform_variadic_list; // all eos's static constexpr const auto full_eos_list = - tl= 0); + EOS eosi = SGAPPLYMODSIMPLE(SimpleMACAW(A, B, Cvinf, v0, T0, Gc)); + if (enabled[3] == 1) { + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); + } + EOS eos_ = SGAPPLYMOD(SimpleMACAW(A, B, Cvinf, v0, T0, Gc)); + eos[matindex] = eos_.GetOnDevice(); + return 0; +} + +int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, + const double Cvinf, const double v0, const double T0, + const double Gc) { + return init_sg_SimpleMACAW(matindex, eos, A, B, Cvinf, v0, T0, Gc, def_en, def_v); +} + int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, const double pc, const double Cv, int const *const enabled, diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index e13e3d93eb0..a63627072e8 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -1,5 +1,5 @@ !------------------------------------------------------------------------------ -! © 2021-2024. Triad National Security, LLC. All rights reserved. This +! © 2021-2025. Triad National Security, LLC. All rights reserved. This ! program was produced under U.S. Government contract 89233218CNA000001 ! for Los Alamos National Laboratory (LANL), which is operated by Triad ! National Security, LLC for the U.S. Department of Energy/National @@ -35,6 +35,7 @@ module singularity_eos init_sg_IdealGas_f,& init_sg_Gruneisen_f,& init_sg_JWL_f,& + init_sg_SimpleMACAW_f,& init_sg_DavisProducts_f,& init_sg_DavisReactants_f,& init_sg_NobleAbel_f,& @@ -112,6 +113,19 @@ end function init_sg_Gruneisen end function init_sg_JWL end interface + interface + integer(kind=c_int) function & + init_sg_SimpleMACAW(matindex, eos, A, B, Cvinf, v0, T0, Gc, sg_mods_enabled, & + sg_mods_values) & + bind(C, name='init_sg_SimpleMACAW') + import + integer(c_int), value, intent(in) :: matindex + type(c_ptr), value, intent(in) :: eos + real(kind=c_double), value, intent(in) :: A, B, Cvinf, v0, T0, Gc + type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values + end function init_sg_SimpleMACAW + end interface + interface integer(kind=c_int) function & init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, & @@ -503,6 +517,28 @@ integer function init_sg_JWL_f(matindex, eos, A, B, R1, R2, w, rho0, Cv, & c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) end function init_sg_JWL_f + integer function init_sg_SimpleMACAW_f(matindex, eos, A, B, Cvinf, v0, T0, Gc, & + sg_mods_enabled, sg_mods_values) & + result(err) + integer(c_int), value, intent(in) :: matindex + type(sg_eos_ary_t), intent(in) :: eos + real(kind=8), value, intent(in) :: A, B, Cvinf, v0, T0, Gc + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + + err = init_sg_SimpleMACAW(matindex-1, eos%ptr, A, B, Cvinf, v0, T0, Gc, & + c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) + end function init_sg_SimpleMACAW_f + + integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & Cv, sg_mods_enabled, & sg_mods_values) & diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index d374eb40dd1..1bfef69fd69 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -34,6 +34,10 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, const double R1, const double R2, const double w, const double rho0, const double Cv, int const *const enabled, double *const vals); +int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, + const double Cvinf, const double v0, const double T0, + const double Gc, int const *const enabled, double *const vals); + int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const double s1, const double s2, const double s3, const double G0, const double b, const double rho0, const double T0, const double P0, @@ -144,6 +148,10 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, const double R1, const double R2, const double w, const double rho0, const double Cv); +int init_sg_SimpleMACAW(const int matindex, EOS *eos, const double A, const double B, + const double Cvinf, const double v0, const double T0, + const double Gc); + int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const double s1, const double s2, const double s3, const double G0, const double b, const double rho0, const double T0, const double P0, diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 903983d8051..cf857633975 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -20,6 +20,7 @@ add_executable( test_eos_sap_polynomial.cpp test_eos_noble_abel.cpp test_eos_stiff.cpp + test_eos_simple_macaw.cpp test_eos_carnahan_starling.cpp test_eos_vinet.cpp test_eos_mgusup.cpp diff --git a/test/test_eos_simple_macaw.cpp b/test/test_eos_simple_macaw.cpp new file mode 100644 index 00000000000..f63110a18cb --- /dev/null +++ b/test/test_eos_simple_macaw.cpp @@ -0,0 +1,179 @@ +//------------------------------------------------------------------------------ +// © 2021-2025. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#include +#include +#include +#include +#include +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif +#include + +#include +#include +#include + +using singularity::SimpleMACAW; +using EOS = singularity::Variant; + +constexpr Real REAL_TOL = 1e-12; + +// Only run exception checking when we aren't offloading to the GPUs +SCENARIO("Testing the Simple MACAW EOS", "[SimpleMACAWEOS]") { + GIVEN("Parameters for a Simple MACAW EOS") { + // Parameters for copper (See Table 1 in the paper by Aslam & Lozano) + constexpr Real A = 7.3; + constexpr Real B = 3.9; + constexpr Real Cvinf = 0.000389; + constexpr Real v0 = 1. / 8.952; + constexpr Real T0 = 150.; + constexpr Real Gc = 0.5; + + // Create the EOS + auto host_eos = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); + auto eos = host_eos.GetOnDevice(); + + WHEN("A set of densities is provided") { + Real rho = 0.5; + THEN("The temperatue at this density and an energy on the cold curve should be " + "zero") { + for (int i = 0; i < 10; i++) { + rho += rho + static_cast(i); // cylce through a variety of densities + const Real e = eos.SieColdCurve(rho); + INFO("i: " << i << " rho = " << rho << " e = " << e); + CHECK_THAT(eos.TemperatureFromDensityInternalEnergy(rho, e), + Catch::Matchers::WithinRel(0.0, 1.0e-12)); + } // for + } // Then + } // When + + WHEN("A density and temperature are provided") { + Real rho = 0.56; + Real T = 123.0; + for (int i = 0; i < 10; i++) { + rho += i; // cylce through a variety of densities + T += 15.0 * i; + DYNAMIC_SECTION("For a given density, " << rho << ", and temperature, " << T) { + // Setup thermodynamic derivatives + Real cp = eos.ConstantPressureSpecificHeatFromDensityTemperature(rho, T); + Real cv = eos.SpecificHeatFromDensityTemperature(rho, T); + Real Bs = eos.BulkModulusFromDensityTemperature(rho, T); + Real BT = eos.IsothermalBulkModulusFromDensityTemperature(rho, T); + Real beta = eos.CoefficientThermalExpansionFromDensityTemperature(rho, T); + Real G = eos.GruneisenParamFromDensityTemperature(rho, T); + + // Test a battery of different thermodynamic identity tests + INFO("rho = " << rho << ", T = " << T); + INFO("cp = " << cp << ", cv = " << cv << ", Bs = " << Bs << ", BT = " << BT); + INFO("beta = " << beta << ", G = " << G); + THEN("The thermodynamic stability ensures: cp >= cv and Bs >= BT") { + REQUIRE(Bs >= BT); + REQUIRE(cp >= cv); + } // Then + + Real term1 = cp * BT / (cv * Bs); + INFO("cp * BT / (cv * Bs) = " << term1); + THEN("The thermodynamic equality satisfies: cp * BT / (cv * Bs) = 1") { + REQUIRE(isClose(term1, 1.0, REAL_TOL)); + } // Then + + Real term2 = BT / Bs + beta * beta * T * BT / (rho * cp); + INFO("cv / cp + beta * beta * T * BT / (rho * cp) = " << term2); + THEN("The thermodynamic equality satisfies: cv / cp + beta^2 * T * BT / (rho * " + "cp)") { + REQUIRE(isClose(term2, 1.0, REAL_TOL)); + } // Then + + Real gruneisen = beta * BT / (rho * cv); + INFO("beta * BT / (rho * cv) = " << gruneisen); + THEN("The thermodynamic equality for the Gruneisen parameter satisfies: Gamma " + "= beta * BT / (rho * cv)") { + REQUIRE(isClose(G, gruneisen, REAL_TOL)); + } // Then + + } // Dynamic Section + } // for + } // When + + host_eos.Finalize(); + eos.Finalize(); + } // Given +} // Scenario + +SCENARIO("Testing the Variant API Simple MACAW EOS", "[SimpleMACAWEOS]") { + GIVEN("Parameters for a Simple MACAW EOS") { + // Parameters for copper (See Table 1 in the paper by Aslam & Lozano) + constexpr Real A = 7.3; + constexpr Real B = 3.9; + constexpr Real Cvinf = 0.000389; + constexpr Real v0 = 1. / 8.952; + constexpr Real T0 = 150.; + constexpr Real Gc = 0.5; + + // Create the EOS + using EOS = singularity::Variant; + EOS eos_in_variant = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); + + WHEN("Densities are provided") { + // For large temperatures, the specific heat capacity should approach a constant + // value (The Dulong-Petit Law) + Real rho = 0.1; + Real T = 1e7; + Real Cv = eos_in_variant.SpecificHeatFromDensityTemperature(rho, T); + for (int i = 0; i < 10; i++) { + rho += rho + i; // cylce through a variety of densities + DYNAMIC_SECTION("For a given density " << rho << " and temperature " << T) { + INFO(std::fixed << std::setprecision(15) << "Cv = " << Cv << ", Cvinf = " + << Cvinf << ", T = " << T << ", rho = " << rho); + THEN("The Dulong-Petit law should be satisfied") { + REQUIRE(isClose(Cv, Cvinf, 1e-8 * Cvinf)); + } // Then + } // Dynamic section + } // for + + // Check for zero specific heat at zero temperature + Cv = eos_in_variant.SpecificHeatFromDensityTemperature(rho, 0.); + THEN("The specific heat capacity at constant volume should be zero at zero " + "temperature") { + INFO(std::fixed << std::setprecision(15) << "Cv = " << Cv << ", rho = " << rho); + REQUIRE(Cv == 0.0); + } + } + + // Check that the EOS is inverting correctly + WHEN("A set of densities and energies are provided") { + using EOS = singularity::Variant; + EOS eos_in_variant = SimpleMACAW(A, B, Cvinf, v0, T0, Gc); + + const Real rho_0 = 0.3; + const Real sie_0 = 190.0; + const Real P = eos_in_variant.PressureFromDensityInternalEnergy(rho_0, sie_0); + const Real T = eos_in_variant.TemperatureFromDensityInternalEnergy(rho_0, sie_0); + Real rho, sie; + Real *lambda; + eos_in_variant.DensityEnergyFromPressureTemperature(P, T, lambda, rho, sie); + INFO("rho_0 = " << rho_0 << ", rho = " << rho << ", sie_0 = " << sie_0 + << ", sie = " << sie); + INFO("P = " << P << " T = " << T); + THEN("The inversion back to rho and sie from P and T should be identital.") { + REQUIRE(isClose(rho, rho_0, REAL_TOL * rho_0)); + REQUIRE(isClose(sie, sie_0, 100 * REAL_TOL * sie_0)); + } + } // When + eos_in_variant.Finalize(); + } // Given +} // Scenario