diff --git a/CHANGELOG.md b/CHANGELOG.md index 1f91381416..23369ea421 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,8 @@ ## Current develop ### Fixed (Repair bugs, etc) +### Added (new features/APIs/variables/...) +- [[PR331]](https://github.com/lanl/singularity-eos/pull/331) Included code and documentation for a full, temperature consistent, Mie-Gruneisen EOS based on a linear Us-up relation. ### Added (new features/APIs/variables/...) - [[PR326]](https://github.com/lanl/singularity-eos/pull/326) Document how to do a release diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 26846724af..e32f2a5e47 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -580,6 +580,9 @@ Given the inconsisetency in the temperature, we have made the choice **not** to expose the entropy for this EOS. **Requesting an entropy value will result in an error.** +If a linear :math:`U_s`-:math:`u_p` relation is enough for your problem, we recommend using the MGUsup +EOS described below. It is a complete EOS with consistent temperature. + Given a reference density, :math:`\rho_0`, we first parameterize the EOS using :math:`\eta` as a measure of compression given by @@ -751,6 +754,148 @@ 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. +Mie-Gruneisen linear :math:`U_s`- :math:`u_p` EOS +````````````````````````````````````````````````` + +One of the most commonly-used EOS is the linear :math:`U_s`- :math:`u_p` version of the Mie-Gruneisen EOS. This EOS +uses the Hugoniot as the reference curve and is extensively used in shock physics. +This version implements the exact thermodynamic temperature on the Hugoniot and also adds an entropy. + +The pressure follows the traditional Mie-Gruneisen form, + +.. math:: + + P(\rho, e) = P_H(\rho) + \rho\Gamma(\rho) \left(e - e_H(\rho) \right), + +Here the subscript :math:`H` is a reminder that the reference curve is a +Hugoniot. :math:`\Gamma` is the Gruneisen parameter and the first approximation +is that :math:`\rho\Gamma(\rho)=\rho_0\Gamma(\rho_0)` +which is the same assumption as in the Gruneisen EOS when :math:`b=0`. + +The above is an incomplete equation of state because it only relates the +pressure to the density and energy, the minimum required in a solution to the +Euler equations. To complete the EOS and determine the temperature and entropy, a constant +heat capacity is assumed so that + +.. math:: + + T(\rho, e) = \frac{\left(e-e_H(\rho)\right)}{C_V} + T_H(\rho) + +Note the difference from the Gruneisen EOS described above. We still use a constant :math:`C_V`, +and it is usually taken at the reference temperature, but +we now extrapolate from the temperature on the Hugoniot, :math:`T_H(\rho)`, and not +from the reference temperature, :math:`T_0`. + +With this consistent temperature we can derive an entropy in a similar way as for the Vinet EOS. Using +thermodynamic derivatives we can show that + +.. math:: + + \Gamma \rho = \frac{\alpha B_T}{C_V} , + +and we arrive at + +.. math:: + + S(\rho,T) = S_0 - \Gamma(\rho_0)C_V \eta + {C_V} \ln \frac{T}{T_0} , + + +where :math:`\eta` is a measure of compression given by + +.. math:: + + \eta = 1 - \frac{\rho_0}{\rho}. + +This is convenient because :math:`\eta = 0` when :math:`\rho = \rho_0`, +:math:`\eta = 1` at the infinite density limit, and :math:`\eta = -\infty` at +the zero density limit. + +The pressure, energy, and temperature, on the Hugoniot are derived from the +shock jump conditions, + +.. math:: + + \rho_0 U_s &= \rho (U_s - u_p) \\ + P_H &= \rho_0 U_s u_p \ , + +assuming a linear :math:`U_s`- :math:`u_p` relation, + +.. math:: + + U_s = C_s + s u_p . + +Here :math:`U_s` is the shock velocity and :math:`u_p` is the particle +velocity. As is pointed out in the description of the Gruneisen EOS, +for many materials, the :math:`U_s`- :math:`u_p` relationship is roughly linear +so only this :math:`s` parameter is needed. The units for :math:`C_s` is velocity while +:math:`s` is unitless. Note that the parameter :math:`s` is related to the +fundamental derivative of shock physics as shown by `Wills `_. + +Solving the jump equations above gives that the reference pressure along the Hugoniot is determined by + +.. math:: + + P_H(\rho) = C_s^2 \rho_0 \frac{\eta}{\left(1 - s \eta \right)^2} . + +Note the singularity at :math:`s \eta = 1` which limits this model's validity to compressions +:math:`\eta << 1/s`. If your problem can be expected to have compressions of this order, you should use the PowerMG +EOS that is explicitely constructed for large compressions. +The assumption of linear :math:`U_s`- :math:`u_p` relation is simply not valid at large compressions. + +The energy along the Hugoniot is given by + +.. math:: + + E_H(\rho) = \frac{P_H \eta }{2 \rho_0} + E_0 . + +The temperature on the Hugoniot is hard to derive explicitely but with the help of Mathematica +we can solve + +.. math:: + + T_H(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{e^{\Gamma(\rho_0) \eta}}{2 C_V \rho_0} + \int_0^\eta e^{-\Gamma(\rho_0) z} z^2 \frac{d}{dz} \left( \frac{P_H}{z}\right) dz + + +into the explicit formula + +.. math:: + + T_H(\rho) &= T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}{2 C_V s^2} + \left[\frac{- s \eta}{(1 - s \eta)^2} + \left( \frac{\Gamma(\rho_0)}{s} - 3 \right) + \left( e^{\Gamma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right)\right. \\ + & \ \left. + e^{-\frac{\Gamma(\rho_0)}{s} (1-s \eta)} + \left( Ei(\frac{\Gamma(\rho_0)}{s}(1-s \eta))-Ei(\frac{\Gamma(\rho_0)}{s}) \right) + \left((\frac{\Gamma(\rho_0)}{s})^2 - 4 \frac{\Gamma(\rho_0)}{s} + 2 \right) \right] + +where :math:`Ei` is the exponential integral function. We replace the :math:`Ei` difference with a sum with cutoff +giving an error less than machine precision. For :math:`s \eta` close to :math:`0`, there are +severe cancellations in this formula and we use the expansion + +.. math:: + + {T_H}_{exp}(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}{2 C_V s^2} + \left[ -2 \ln ( 1- s \eta) + \frac{s \eta}{(1 - s \eta)^2} ( 3 s \eta - 2) \right] \ . + + +The first omitted term in the expansion inside the square brackets is :math:`\Gamma(\rho_0) \eta^4 / 6`. This expansion is +in fact even better than the common approximation of replacing the full temperature on the Hugoniot with the temperature on the +isentrope, that is, the first term :math:`T_0 e^{\Gamma(\rho_0) \eta}`. + +The constructor for the ``MGUsup`` EOS has the signature + +.. code-block:: cpp + + 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) + +where +``rho0`` is :math:`\rho_0`, ``T0`` is :math:`T_0`, +``Cs`` is :math:`C_s`, ``s`` is :math:`s`, +``G0`` is :math:`\Gamma(\rho_0)`, ``Cv0`` is :math:`C_V`, +``E0`` is :math:`E_0`, and ``S0`` is :math:`S_0`. + + JWL EOS `````````` diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp index 69a839235a..f5d45382f6 100644 --- a/singularity-eos/eos/eos.hpp +++ b/singularity-eos/eos/eos.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2024. 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 @@ -39,13 +39,13 @@ #include #include #include +#include #include #include #include #include #include #include - // Modifiers #include #include @@ -67,7 +67,7 @@ using singularity::variadic_utils::transform_variadic_list; // all eos's static constexpr const auto full_eos_list = - tl +#include + +#include + +#include +#include +#include +#include +#include + +namespace singularity { + +using namespace eos_base; +using namespace robust; + +// COMMENT: This is a complete, thermodynamically consistent Mie-Gruneisen +// EOS based on a linear Us-up Hugoniot, Us=Cs+s*up. +class MGUsup : public EosBase { + public: + 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) { + _CheckMGUsup(); + } + + MGUsup GetOnDevice() { return *this; } + PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real + MinInternalEnergyFromDensity(const Real rho, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const { + return _Cv0; + } + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return _Cv0; + } + // added for testing AEM Dec 2023 + PORTABLE_INLINE_FUNCTION Real HugPressureFromDensity(const Real rho) const; + PORTABLE_INLINE_FUNCTION Real HugInternalEnergyFromDensity(const Real rho) const; + PORTABLE_INLINE_FUNCTION Real HugTemperatureFromDensity(const Real rho) const; + // Thermal Bulk Modulus added AEM Dec 2022 + PORTABLE_INLINE_FUNCTION Real TBulkModulusFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; + // Thermal expansion coefficient added AEM 2022 + PORTABLE_INLINE_FUNCTION Real TExpansionCoeffFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const { + return robust::ratio(_G0 * _rho0, rho); + } + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return robust::ratio(_G0 * _rho0, rho); + } + PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, + Real &cv, Real &bmod, const unsigned long output, + Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION + void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, + Real &bmod, Real &dpde, Real &dvdt, + Real *lambda = 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_BASE_CLASS_USINGS(MGUsup) + PORTABLE_INLINE_FUNCTION + int nlambda() const noexcept { return 0; } + static constexpr unsigned long PreferredInput() { return _preferred_input; } + static inline unsigned long scratch_size(std::string method, unsigned int nelements) { + return 0; + } + static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } + PORTABLE_INLINE_FUNCTION void PrintParams() const { + 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); + printf("\n\n"); + } + // Density/Energy from P/T not unique, if used will give error + PORTABLE_INLINE_FUNCTION void + DensityEnergyFromPressureTemperature(const Real press, const Real temp, Real *lambda, + Real &rho, Real &sie) const; + inline void Finalize() {} + static std::string EosType() { return std::string("MGUsup"); } + static std::string EosPyType() { return EosType(); } + + private: + static constexpr const unsigned long _preferred_input = + thermalqs::density | thermalqs::specific_internal_energy; + Real _rho0, _T0, _Cs, _s, _G0, _Cv0, _E0, _S0; + void _CheckMGUsup(); +}; + +inline void MGUsup::_CheckMGUsup() { + + if (_rho0 < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter rho0 < 0"); + } + if (_T0 < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter T0 < 0"); + } + if (_Cs < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter Cs < 0"); + } + if (_s < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter s < 0"); + } + if (_G0 < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter G0 < 0"); + } + if (_Cv0 < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter Cv0 < 0"); + } +} + +PORTABLE_INLINE_FUNCTION Real MGUsup::HugPressureFromDensity(Real rho) const { + Real eta = 1.0 - robust::ratio(_rho0, rho); + return _Cs * _Cs * _rho0 * robust::ratio(eta, (1.0 - _s * eta) * (1.0 - _s * eta)); +} +PORTABLE_INLINE_FUNCTION Real MGUsup::HugInternalEnergyFromDensity(Real rho) const { + Real eta = 1.0 - robust::ratio(_rho0, rho); + return _E0 + robust::ratio(eta, 2.0 * _rho0) * HugPressureFromDensity(rho); +} +PORTABLE_INLINE_FUNCTION Real MGUsup::HugTemperatureFromDensity(Real rho) const { + int sumkmax = 20; + Real cutoff = 1.e-16; + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real f1 = 1.0 - _s * eta; + if (f1 <= 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("MGUsup model parameters s and rho0 together with rho " + "give a negative argument for a logarithm."); + } + Real G0os = robust::ratio(_G0, _s); + // sk, lk, mk, sum, and enough values may change in the loop + Real sk = -_G0 * eta; + Real lk = _G0; + Real mk = 0.0; + Real sum = sk; + int enough = 0; + Real temp; + Real pf = _Cs * _Cs / (2.0 * _Cv0 * _s * _s); + if (eta * eta < 1.e-8) { + temp = _T0 * exp(_G0 * eta) + + pf * (2.0 * std::log(f1) - _s * eta / (f1 * f1) * (3.0 * _s * eta - 2.0)); + } else { + for (int i = 0; ((i < sumkmax) && (enough == 0)); i++) { + mk = G0os / (i + 2) / (i + 2); + sk = (sk * f1 - lk * eta) * (i + 1) * mk; + lk = mk * (i + 1) * lk; + sum = sum + sk; + if (robust::ratio((sk * sk), (sum * sum)) < (cutoff * cutoff)) { + enough = 1; + } + } + if (enough == 0) { +#ifndef NDEBUG + PORTABLE_WARN("Hugoniot Temperature not converged"); +#endif // NDEBUG + } + // printf("sum=%e\n",sum); + temp = _T0 - pf * ((G0os - 3.0) + exp(-G0os) * (G0os * G0os - 4.0 * G0os + 2.0) * + (std::log(f1) + sum)); + temp = temp * exp(_G0 * eta); + temp = temp - pf / f1 * ((4.0 - G0os) - 1.0 / f1); + } + return temp; +} + +PORTABLE_INLINE_FUNCTION Real MGUsup::InternalEnergyFromDensityTemperature( + const Real rho, const Real temp, Real *lambda) const { + Real value = + HugInternalEnergyFromDensity(rho) + _Cv0 * (temp - HugTemperatureFromDensity(rho)); + return value; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::PressureFromDensityTemperature(const Real rho, + const Real temp, + Real *lambda) const { + Real value = HugPressureFromDensity(rho) + + _G0 * _rho0 * _Cv0 * (temp - HugTemperatureFromDensity(rho)); + return value; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::EntropyFromDensityTemperature(const Real rho, + const Real temp, + Real *lambda) const { + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real value = _S0 - _G0 * _Cv0 * eta + _Cv0 * std::log(temp / _T0); + return value; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::TExpansionCoeffFromDensityTemperature( + const Real rho, const Real temp, Real *lambda) const { + Real value = + robust::ratio(_Cv0 * _rho0 * _G0, TBulkModulusFromDensityTemperature(rho, temp)); + return value; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::TBulkModulusFromDensityTemperature( + const Real rho, const Real temp, Real *lambda) const { + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real value = robust::ratio((1.0 + _s * eta - _G0 * _s * eta * eta), (1.0 - _s * eta)); + if (eta == 0.0) { + value = _Cs * _Cs * _rho0; + } else { + value = value * robust::ratio(HugPressureFromDensity(rho), eta); + } + value = value - _G0 * _G0 * _Cv0 * _rho0 * HugTemperatureFromDensity(rho); + value = robust::ratio(_rho0, rho) * value; + return value; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::BulkModulusFromDensityTemperature( + const Real rho, const Real temp, Real *lambda) const { + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real value = + robust::ratio((1.0 + _s * eta - _G0 * _s * eta * eta), eta * (1.0 - _s * eta)); + if (eta == 0.0) { + value = _Cs * _Cs * _rho0; + } else { + value = value * robust::ratio(HugPressureFromDensity(rho), eta); + } + value = value - _G0 * _G0 * _Cv0 * _rho0 * (temp - HugTemperatureFromDensity(rho)); + value = robust::ratio(_rho0, rho) * value; + return value; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::TemperatureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda) const { + Real value = + (sie - HugInternalEnergyFromDensity(rho)) / _Cv0 + HugTemperatureFromDensity(rho); + if (value < 0.0) { +#ifndef NDEBUG + PORTABLE_WARN("Negative temperature"); +#endif // NDEBUG + value = 1.e-12; + } + return value; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::PressureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda) const { + Real value = HugPressureFromDensity(rho) + + _rho0 * _G0 * (sie - HugInternalEnergyFromDensity(rho)); + return value; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::MinInternalEnergyFromDensity(const Real rho, + Real *lambda) const { + MinInternalEnergyIsNotEnabled("MGUsup"); + return 0.0; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda) const { + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real value = std::log(TemperatureFromDensityInternalEnergy(rho, sie) / _T0); + value = _S0 - _G0 * _Cv0 * eta + _Cv0 * value; + if (value < 0.0) { +#ifndef NDEBUG + PORTABLE_WARN("Negative entropy"); +#endif // NDEBUG + value = 1.e-12; + } + return value; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::BulkModulusFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda) const { + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real value = robust::ratio((1.0 + _s * eta - _G0 * _s * eta * eta), (1.0 - _s * eta)); + if (eta == 0.0) { + value = _Cs * _Cs * _rho0; + } else { + value = value * robust::ratio(HugPressureFromDensity(rho), eta); + } + value = value + _G0 * _G0 * _rho0 * (sie - HugInternalEnergyFromDensity(rho)); + value = robust::ratio(_rho0, rho) * value; + return value; +} +// AEM: Give error since function is not well defined +PORTABLE_INLINE_FUNCTION void +MGUsup::DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Real *lambda, Real &rho, Real &sie) const { + EOS_ERROR("MGUsup::DensityEnergyFromPressureTemperature: " + "Not implemented.\n"); +} +// AEM: We should add entropy and Gruneissen parameters here so that it is complete +// If we add also alpha and BT, those should also be in here. +PORTABLE_INLINE_FUNCTION void MGUsup::FillEos(Real &rho, Real &temp, Real &sie, + Real &press, Real &cv, Real &bmod, + const unsigned long output, + Real *lambda) const { + const unsigned long input = ~output; // everything that is not output is input + if (thermalqs::density & output) { + EOS_ERROR("MGUsup FillEos: Density is required input.\n"); + } else if (thermalqs::temperature & output && + thermalqs::specific_internal_energy & output) { + EOS_ERROR("MGUsup FillEos: Density and Internal Energy or Density and Temperature " + "are required input parameters.\n"); + } + if (thermalqs::temperature & input) { + sie = InternalEnergyFromDensityTemperature(rho, temp); + } + if (output & thermalqs::temperature) + temp = TemperatureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::specific_internal_energy) sie = sie; + if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::specific_heat) + cv = SpecificHeatFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::bulk_modulus) + bmod = BulkModulusFromDensityInternalEnergy(rho, sie); +} + +// TODO(JMM): pre-cache these rather than recomputing them each time +PORTABLE_INLINE_FUNCTION +void MGUsup::ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, + Real &cv, Real &bmod, Real &dpde, Real &dvdt, + Real *lambda) const { + // AEM: Added all variables I think should be output eventually + Real tbmod; + // Real entropy, alpha, Gamma; + + rho = _rho0; + temp = _T0; + sie = _E0; + press = PressureFromDensityTemperature(rho, temp, lambda); + // entropy = _S0; + cv = _Cv0; + tbmod = _Cs * _Cs * _rho0 - _G0 * _G0 * _Cv0 * _rho0 * _T0; + // alpha = _A0; + bmod = _Cs * _Cs * _rho0; + // Gamma = _G0; + // AEM: I suggest taking the two following away. + dpde = _G0 * _rho0; + dvdt = robust::ratio(-1.0, _T0 * _G0 * _rho0); +} +} // namespace singularity + +#endif // _SINGULARITY_EOS_EOS_MGUSUP_HPP_ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d63bbbb51b..6add5dbd39 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -21,6 +21,7 @@ add_executable( test_eos_noble_abel.cpp test_eos_stiff.cpp test_eos_vinet.cpp + test_eos_mgusup.cpp ) add_executable( diff --git a/test/test_eos_mgusup.cpp b/test/test_eos_mgusup.cpp new file mode 100644 index 0000000000..df35715eac --- /dev/null +++ b/test/test_eos_mgusup.cpp @@ -0,0 +1,527 @@ +//------------------------------------------------------------------------------ +// © 2021-2024. 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 +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +#include +#include +#include +#include +#include +#include + +using singularity::MGUsup; +using EOS = singularity::Variant; + +SCENARIO("MGUsup EOS rho sie", "[VectorEOS][MGUsupEOS]") { + GIVEN("Parameters for a MGUsup EOS") { + // Unit conversions + constexpr Real Mbcc_per_g = 1e12; + // MGUsup parameters for tin beta phase + constexpr Real rho0 = 7.285; + constexpr Real T0 = 298.0; + constexpr Real Cs = 2766.0e2; + constexpr Real s = 1.5344; + constexpr Real G0 = 2.4659; + constexpr Real Cv0 = 0.2149e-05 * Mbcc_per_g; + constexpr Real E0 = 0.658e-03 * Mbcc_per_g; + constexpr Real S0 = 0.4419e-05 * Mbcc_per_g; + // Create the EOS + EOS host_eos = MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0); + EOS eos = host_eos.GetOnDevice(); + MGUsup host_eos2 = MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0); + MGUsup eos2 = host_eos2.GetOnDevice(); + + eos.PrintParams(); + + GIVEN("Densities and energies") { + constexpr int num = 4; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_density("density"); + // AEM: I think it should not be "density" again below + Kokkos::View v_energy("density"); +#else + // Otherwise just create arrays to contain values and create pointers to + // be passed to the functions in place of the Kokkos views + std::array density; + std::array energy; + auto v_density = density.data(); + auto v_energy = energy.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize density and energy", 0, 1, PORTABLE_LAMBDA(int i) { + v_density[0] = 7.0; + v_density[1] = 8.0; + v_density[2] = 9.0; + v_density[3] = 7.285; + v_energy[0] = 4.6e8; + v_energy[1] = 1.146e9; + v_energy[2] = 3.28e9; + v_energy[3] = 0.658e9; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto density = Kokkos::create_mirror_view(v_density); + auto energy = Kokkos::create_mirror_view(v_energy); + Kokkos::deep_copy(density, v_density); + Kokkos::deep_copy(energy, v_energy); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Gold standard values for a subset of lookups + constexpr std::array temperature_true{ + 1.502141159804009e+02, 4.26645103516999e+02, 7.08167099601959e+02, T0}; + constexpr std::array hugtemperature_true{ + 2.684894520258222e+02, 3.90542065829047e+02, 7.78960970777667e+02, T0}; + constexpr std::array pressure_true{ + -2.466829656715215e+10, 6.82999362849669e+10, 2.09379236091689e+11, 0.}; + constexpr std::array hugpressure_true{ + -2.010229955618467e+10, 6.690618533331688e+10, 2.121122201185450e+11, 0.}; + constexpr std::array hugenergy_true{ + 7.141736971616102e+8, 1.068414572008593e+09, 3.432136029156598e+09, E0}; + constexpr std::array entropy_true{ + 3.1626206459885668e+06, 4.7165703738323096e+06, 5.2693499546078807e+06, S0}; + constexpr std::array cv_true{Cv0, Cv0, Cv0, Cv0}; + constexpr std::array bulkmodulus_true{ + 4.38665321162636e+11, 8.776332438691490e+11, 1.465222098683647e+12, + Cs * Cs * rho0}; + constexpr std::array gruneisen_true{G0 * rho0 / 7.0, G0 * rho0 / 8.0, + G0 * rho0 / 9.0, G0}; + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create device views for outputs and mirror those views on the host + Kokkos::View v_temperature("Temperature"); + Kokkos::View v_hugtemperature("HugTemperature"); + Kokkos::View v_pressure("Pressure"); + Kokkos::View v_hugpressure("HugPressure"); + Kokkos::View v_hugenergy("HugEnergy"); + Kokkos::View v_entropy("Entropy"); + Kokkos::View v_cv("Cv"); + Kokkos::View v_bulkmodulus("bmod"); + Kokkos::View v_gruneisen("Gamma"); + auto h_temperature = Kokkos::create_mirror_view(v_temperature); + auto h_hugtemperature = Kokkos::create_mirror_view(v_hugtemperature); + auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_hugpressure = Kokkos::create_mirror_view(v_hugpressure); + auto h_hugenergy = Kokkos::create_mirror_view(v_hugenergy); + auto h_entropy = Kokkos::create_mirror_view(v_entropy); + auto h_cv = Kokkos::create_mirror_view(v_cv); + auto h_bulkmodulus = Kokkos::create_mirror_view(v_bulkmodulus); + auto h_gruneisen = Kokkos::create_mirror_view(v_gruneisen); +#else + // Create arrays for the outputs and then pointers to those arrays that + // will be passed to the functions in place of the Kokkos views + std::array h_temperature; + std::array h_hugtemperature; + std::array h_pressure; + std::array h_hugpressure; + std::array h_hugenergy; + std::array h_entropy; + std::array h_cv; + std::array h_bulkmodulus; + std::array h_gruneisen; + // Just alias the existing pointers + auto v_temperature = h_temperature.data(); + auto v_hugtemperature = h_hugtemperature.data(); + auto v_pressure = h_pressure.data(); + auto v_hugpressure = h_hugpressure.data(); + auto v_hugenergy = h_hugenergy.data(); + auto v_entropy = h_entropy.data(); + auto v_cv = h_cv.data(); + auto v_bulkmodulus = h_bulkmodulus.data(); + auto v_gruneisen = h_gruneisen.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // temporary + + WHEN("A PH(rho) lookup is performed") { + portableFor( + "Test Hugoniot pressure", 0, num, PORTABLE_LAMBDA(const int i) { + v_hugpressure[i] = eos2.HugPressureFromDensity(v_density[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_hugpressure, v_hugpressure); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned PH(rho) should be equal to the true value") { + array_compare(num, density, energy, h_hugpressure, hugpressure_true, "Density", + "Energy"); + } + } + + WHEN("A EH(rho) lookup is performed") { + portableFor( + "Test Hugoniot internal energy", 0, num, PORTABLE_LAMBDA(const int i) { + v_hugenergy[i] = eos2.HugInternalEnergyFromDensity(v_density[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_hugenergy, v_hugenergy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned EH(rho) should be equal to the true value") { + array_compare(num, density, energy, h_hugenergy, hugenergy_true, "Density", + "Energy"); + } + } + + WHEN("A TH(rho) lookup is performed") { + portableFor( + "Test Hugoniot temperature", 0, num, PORTABLE_LAMBDA(const int i) { + v_hugtemperature[i] = eos2.HugTemperatureFromDensity(v_density[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_hugtemperature, v_hugtemperature); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned TH(rho) should be equal to the true value") { + array_compare(num, density, energy, h_hugtemperature, hugtemperature_true, + "Density", "Energy"); + } + } + // end temporary + + WHEN("A T(rho, e) lookup is performed") { + eos.TemperatureFromDensityInternalEnergy(v_density, v_energy, v_temperature, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_temperature, v_temperature); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned T(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_temperature, temperature_true, "Density", + "Energy"); + } + } + + WHEN("A P(rho, e) lookup is performed") { + eos.PressureFromDensityInternalEnergy(v_density, v_energy, v_pressure, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_pressure, v_pressure); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned P(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_pressure, pressure_true, "Density", + "Energy"); + } + } + + WHEN("A S(rho, e) lookup is performed") { + portableFor( + "Test entropy", 0, num, PORTABLE_LAMBDA(const int i) { + v_entropy[i] = + eos2.EntropyFromDensityInternalEnergy(v_density[i], v_energy[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_entropy, v_entropy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned S(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_entropy, entropy_true, "Density", + "Energy"); + } + } + + WHEN("A B_S(rho, e) lookup is performed") { + eos.BulkModulusFromDensityInternalEnergy(v_density, v_energy, v_bulkmodulus, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_bulkmodulus, v_bulkmodulus); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned B_S(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_bulkmodulus, bulkmodulus_true, "Density", + "Energy"); + } + } + + WHEN("A Cv(rho, e) lookup is performed") { + eos.SpecificHeatFromDensityInternalEnergy(v_density, v_energy, v_cv, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_cv, v_cv); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned T(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_cv, cv_true, "Density", "Energy"); + } + } + + WHEN("A Gamma(rho, e) lookup is performed") { + eos.GruneisenParamFromDensityInternalEnergy(v_density, v_energy, v_gruneisen, + num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_gruneisen, v_gruneisen); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned Gamma(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_gruneisen, gruneisen_true, "Density", + "Energy"); + } + } + } + } +} +SCENARIO("MGUsup EOS rho T", "[VectorEOS][MGUsupEOS]") { + GIVEN("Parameters for a MGUsup EOS") { + // Unit conversions + constexpr Real Mbcc_per_g = 1e12; + // MGUsup parameters for copper + constexpr Real rho0 = 7.285; + constexpr Real T0 = 298.0; + constexpr Real Cs = 2766.0e2; + constexpr Real s = 1.5344; + constexpr Real G0 = 2.4659; + constexpr Real Cv0 = 0.2149e-05 * Mbcc_per_g; + constexpr Real E0 = 0.658e-03 * Mbcc_per_g; + constexpr Real S0 = 0.4419e-05 * Mbcc_per_g; + // Create the EOS + EOS host_eos = MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0); + EOS eos = host_eos.GetOnDevice(); + MGUsup host_eos2 = MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0); + MGUsup eos2 = host_eos2.GetOnDevice(); + + eos.PrintParams(); + + GIVEN("Densities and temperatures") { + constexpr int num = 4; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_density("density"); + Kokkos::View v_temperature("temperature"); +#else + // Otherwise just create arrays to contain values and create pointers to + // be passed to the functions in place of the Kokkos views + std::array density; + std::array temperature; + auto v_density = density.data(); + auto v_temperature = temperature.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize density and temperature", 0, 1, PORTABLE_LAMBDA(int i) { + v_density[0] = 7.0; + v_density[1] = 8.0; + v_density[2] = 9.0; + v_density[3] = 7.285; + v_temperature[0] = 150.215; + v_temperature[1] = 426.645; + v_temperature[2] = 708.165; + v_temperature[3] = 298.0; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto density = Kokkos::create_mirror_view(v_density); + auto temperature = Kokkos::create_mirror_view(v_temperature); + Kokkos::deep_copy(density, v_density); + Kokkos::deep_copy(temperature, v_temperature); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Gold standard values for a subset of lookups + constexpr std::array sie_true{4.6000189975811839e+08, + 1.1459997775419691e+09, + 3.2799954879553909e+09, 6.58e+08}; + constexpr std::array pressure_true{ + -2.466826243974248e+10, 6.82999322887127e+10, 2.09379155036952e+11, 0.}; + constexpr std::array entropy_true{ + 3.1626332929526418e+06, 4.7165698524198849e+06, 5.2693435831578374e+06, S0}; + constexpr std::array cv_true{Cv0, Cv0, Cv0, Cv0}; + constexpr std::array tbulkmodulus_true{ + 4.2378339466579810e+11, 8.4064844786200464e+11, 1.4106538914691243e+12, + (Cs * Cs - G0 * G0 * T0 * Cv0) * rho0}; + constexpr std::array alpha_true{ + 9.1095620143267597e-05, 4.5922657969193220e-05, 2.7366607342141916e-05, + G0 * Cv0 / (Cs * Cs - G0 * G0 * T0 * Cv0)}; + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create device views for outputs and mirror those views on the host + Kokkos::View v_energy("Energy"); + Kokkos::View v_pressure("Pressure"); + Kokkos::View v_entropy("Entropy"); + Kokkos::View v_cv("Cv"); + Kokkos::View v_tbulkmodulus("tbmod"); + Kokkos::View v_alpha("alpha"); + auto h_energy = Kokkos::create_mirror_view(v_energy); + auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_entropy = Kokkos::create_mirror_view(v_entropy); + auto h_cv = Kokkos::create_mirror_view(v_cv); + auto h_tbulkmodulus = Kokkos::create_mirror_view(v_tbulkmodulus); + auto h_alpha = Kokkos::create_mirror_view(v_alpha); +#else + // Create arrays for the outputs and then pointers to those arrays that + // will be passed to the functions in place of the Kokkos views + std::array h_energy; + std::array h_pressure; + std::array h_entropy; + std::array h_cv; + std::array h_tbulkmodulus; + std::array h_alpha; + // Just alias the existing pointers + auto v_energy = h_energy.data(); + auto v_pressure = h_pressure.data(); + auto v_entropy = h_entropy.data(); + auto v_cv = h_cv.data(); + auto v_tbulkmodulus = h_tbulkmodulus.data(); + auto v_alpha = h_alpha.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + WHEN("A E(rho, T) lookup is performed") { + eos.InternalEnergyFromDensityTemperature(v_density, v_temperature, v_energy, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_energy, v_energy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned E(rho, T) should be equal to the true value") { + array_compare(num, density, temperature, h_energy, sie_true, "Density", + "Temperature"); + } + } + + WHEN("A P(rho, T) lookup is performed") { + eos.PressureFromDensityTemperature(v_density, v_temperature, v_pressure, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_pressure, v_pressure); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned P(rho, T) should be equal to the true value") { + array_compare(num, density, temperature, h_pressure, pressure_true, "Density", + "Temperature"); + } + } + + portableFor( + "entropy", 0, num, PORTABLE_LAMBDA(const int i) { + v_entropy[i] = + eos2.EntropyFromDensityInternalEnergy(v_density[i], v_energy[i]); + }); + WHEN("A S(rho, T) lookup is performed") { + portableFor( + "entropy", 0, num, PORTABLE_LAMBDA(const int i) { + v_entropy[i] = + eos2.EntropyFromDensityTemperature(v_density[i], v_temperature[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_entropy, v_entropy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned S(rho, T) should be equal to the true value") { + array_compare(num, density, temperature, h_entropy, entropy_true, "Density", + "Temperature"); + } + } + + WHEN("A B_T(rho, T) lookup is performed") { + portableFor( + "bulk modulus", 0, num, PORTABLE_LAMBDA(const int i) { + v_tbulkmodulus[i] = + eos2.TBulkModulusFromDensityTemperature(v_density[i], v_temperature[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_tbulkmodulus, v_tbulkmodulus); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned B_T(rho, T) should be equal to the true value") { + array_compare(num, density, temperature, h_tbulkmodulus, tbulkmodulus_true, + "Density", "Temperature"); + } + } + + WHEN("A Cv(rho, T) lookup is performed") { + eos.SpecificHeatFromDensityTemperature(v_density, v_temperature, v_cv, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_cv, v_cv); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned T(rho, e) should be equal to the true value") { + array_compare(num, density, temperature, h_cv, cv_true, "Density", + "Temperature"); + } + } + + WHEN("A alpha(rho, T) lookup is performed") { + portableFor( + "alpha", 0, num, PORTABLE_LAMBDA(const int i) { + v_alpha[i] = eos2.TExpansionCoeffFromDensityTemperature(v_density[i], + v_temperature[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_alpha, v_alpha); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned alpha(rho, T) should be equal to the true value") { + array_compare(num, density, temperature, h_alpha, alpha_true, "Density", + "Temperature"); + } + } + } + } +} +SCENARIO("MGUsup EOS SetUp", "[VectorEOS][MGUsupEOS]") { + GIVEN("Parameters for a MGUsup EOS") { + // Unit conversions + constexpr Real Mbcc_per_g = 1e12; + // MGUsup parameters for copper + Real rho0 = 7.285; + Real T0 = 298.0; + Real Cs = 2766.0e2; + Real s = 1.5344; + Real G0 = 2.4659; + Real Cv0 = 0.2149e-05 * Mbcc_per_g; + constexpr Real E0 = 0.658e-03 * Mbcc_per_g; + constexpr Real S0 = 0.4419e-05 * Mbcc_per_g; + // Create the EOS + EOS host_eos = MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0); + EOS eos = host_eos.GetOnDevice(); + eos.Finalize(); + + WHEN("Faulty/not set parameter rho0 is given") { + rho0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); + THEN("An error message should be written out") {} + } + WHEN("Faulty/not set parameter T0 is given") { + T0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); + THEN("An error message should be written out") {} + } + WHEN("Faulty/not set parameter Cs is given") { + Cs = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); + THEN("An error message should be written out") {} + } + WHEN("Faulty/not set parameter s is given") { + s = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); + THEN("An error message should be written out") {} + } + WHEN("Faulty/not set parameter G0 is given") { + G0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); + THEN("An error message should be written out") {} + } + WHEN("Faulty/not set parameter Cv0 is given") { + Cv0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); + THEN("An error message should be written out") {} + } + } +}