diff --git a/CHANGELOG.md b/CHANGELOG.md index b0c8463243d..6537ac386d1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,7 @@ ### Added (new features/APIs/variables/...) - [[PR556]](https://github.com/lanl/singularity-eos/pull/556) Add introspection into types available in the variant - [[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 +- [[PR526]](https://github.com/lanl/singularity-eos/pull/526) Added the MultiEOS class to wrap multiple EOS models that are linked via PTE (e.g. a multiphase material or reacting material). ### Fixed (Repair bugs, etc) - [[PR567]](https://github.com/lanl/singularity-eos/pull/567) Fixed an OOB array access bug in the Fixed T PTE solver diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 0bc2cf1350f..bc20519d4ec 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -21,6 +21,10 @@ .. _PowerMG: https://www.osti.gov/biblio/1762624 +.. _Sesame: https://www.lanl.gov/org/ddste/aldsc/theoretical/physics-chemistry-materials/sesame-database.php + +.. _EOSPAC: https://laws.lanl.gov/projects/data/eos/eospacReleases.php + EOS Models =========== @@ -2189,12 +2193,13 @@ EOSPAC EOS ```````````` .. warning:: + Entropy is not yet available for this EOS This is a striaghtforward wrapper of the `EOSPAC`_ library for the `Sesame`_ database. The constructor for the ``EOSPAC`` model has several overloads -.. code-block:: +.. code-block:: cpp EOSPAC(int matid, TableSplit split, bool invert_at_setup = false, Real insert_data = 0.0, @@ -2237,8 +2242,240 @@ curve plus Cowan-nuclear model for ions and the final option Note for performance reasons this EOS uses a slightly different vector API. See :ref:`EOSPAC Vector Functions ` for more details. -.. _Sesame: https://www.lanl.gov/org/ddste/aldsc/theoretical/physics-chemistry-materials/sesame-database.php +MultiEOS +````````` + +.. warning:: + + The Serialization/Deserialization capability for this EOS appears bugged + at the moment, especially with EOSPAC. Use with caution. + +This class allows multiple EOS models to be combined *at compile time* with a +PTE closure model into a single material. This is useful for simulating e.g. +reacting systems where reactants and products describe the same material. It +can also be used for multi-phase EOS, but again the number of phases and +crucially the *type* of EOS model must be known at compile time. The advantage +of combining these models at compile-time is that logically the EOS will +behave identically to any other EOS model in ``singularity-eos``. It's possible +that there are performance benefits, especially when all EOS are the same model. + +The primary way of constructing the ``MultiEOS`` model involves simply +providing an arbitrary series of model objects + +.. code-block:: cpp + + MultiEOS(EOSModelsT_ &&...eos_models) + +Note that the type of the ``MultiEOS`` model is deduced from the inputs to the +constructor. In the follwoing example, it would be +``MultiEOS``: + +.. code-block:: cpp + + auto model1 = IdealGas(Gamma_1, Cv_1); + auto model2 = IdealGas(Gamma_2, Cv_2); + auto multi = MultiEOS(model1, model2); + +The ``MultiEOS`` model can also optionally accept a mass fraction cutoff where +materials below the cutoff value are ignored for the PTE caluclation + +.. code-block:: cpp + + MultiEOS(Real mass_frac_cutoff_in, EOSModelsT_ &&...eos_models); + +For all of the standard ``singularity-eos`` member functions provided through +the ``Variant`` class (see +:ref:`The Equation of State API `), bulk +properties are returned for the material and no component densities or energies +are available. + +Mass fractions of the individual components are provided via the ``lambda`` +variable via type-based indexing. + +.. note:: + + Mass fraction information **must** be provided via the ``lambda`` parameter + or a compile-time error will be generated. + +For example, the ``IndexerUtils`` in ``singularity-eos`` provide a way to +generate a type-based lambda indexer. For example, a two component spiner EOS +mixture could use a lambda of the type + +.. code-block:: cpp + + using LambdaT = + VariadicIndexer, IndexableTypes::MassFraction<1>>; + +A helper macro is provided that uses SFINAE to generate a compile time error +when the appropriate number of mass fractions are not provided via type-based +indexing in the ``LambdaIndexer`` type: + +.. code-block:: cpp + + SINGULARITY_INDEXER_HAS_MASS_FRAC(LambdaIndexer, nmat) + +where ``nmat`` is the number of materials in the ``MultiEOS`` class (known at +compile time). + +All derivative quantities (i.e. the Gruneisen parameter, bulk modulus, and heat +capacity) are found via finite differences. The reason for this is that really +there is no accurate way to average the individual material derivatives. For +example, we can attempt to derive an averaging rule for the heat capacity. The +heat capacity is given by + +.. math:: + + C_V = \left( \frac{\partial e}{\partial T} \right)_V, + +i.e. the derivative of the specific internal energy with respect to temperature +at constant *total* volume. The total internal energy is given by + +.. math:: + + e_\mathrm{tot} = \sum\limits_{m=0} \mu_m e_m, + +where the subscript :math:`m` is the material index in the PTE system (in this +case the material models in the ``MultiEOS`` class) and :math:`\mu_m` is the +mass fraction of that material. Taking the derivative of this equation with +respect to temperature at constant volume gives + +.. math:: + + \left( \frac{\partial e_\mathrm{tot}}{\partial T} \right)_{V_\mathrm{tot}} + = \sum\limits_{m=0} \mu_m + \left( \frac{\partial e_m}{\partial T} \right)_{V_\mathrm{tot}} + +But note that the right-hand-side derivatives are at constant *total* volume, +not constant *material* volume. As the temperature increases, a single +material may expand and do work on adjacent materials, causing their energies to +increase more than would be expected for the given temperature increase while +the materials remain in pressure equilibrium. + +The best way to express this behavior analytically would be to mass-average the +derivatives with respect to pressure and temperature, i.e. + +.. math:: + + \left( \frac{\partial e_\mathrm{tot}}{\partial T} \right)_P + = \sum\limits_{m=0} \mu_m + \left( \frac{\partial e_m}{\partial T} \right)_P, + + +and then use thermodynamic identities to combine them to describe the bulk +behavior of the mixture. However, not all models in ``singularity-eos`` provide +pressure and temperature derivatives at present, so these quantities are not +readily accessible. + +Instead, the ``MultiEOS`` class relies on perturbing the mixture in density and +temperature or density and energy space and then combining these finite +difference results via thermodynamic identities. The downside of this approach +is that any finite difference errors can be significantly magnified in the +thermodynamic identities. + +When access to the bare class is available, additional public member functions +are provided that populate indexers (i.e. anything with a square bracket +operator) for the material densities and specific internal energies. These are + +.. code-block:: cpp + + template + PORTABLE_INLINE_FUNCTION SolverStatus GetStatesFromDensityEnergy( + const Real density_tot, const Real sie_tot, Real &pressure, Real &temperature, + RealIndexer &density_mat, RealIndexer &sie_mat, LambdaIndexer &lambdas, + bool const doing_derivs = false, bool small_mass_mat_consistency = false); + +.. code-block:: cpp + + template + PORTABLE_INLINE_FUNCTION SolverStatus GetStatesFromDensityTemperature( + const Real density_tot, Real &sie_tot, Real &pressure, const Real temperature, + RealIndexer &&density_mat, RealIndexer &&sie_mat, LambdaIndexer &lambdas, + bool const doing_derivs = false); + +.. code-block:: cpp + + template + PORTABLE_INLINE_FUNCTION SolverStatus GetStatesFromDensityPressure( + const Real density_tot, Real &sie_tot, const Real pressure, Real &temperature, + RealIndexer &&density_mat, RealIndexer &&sie_mat, LambdaIndexer &lambdas, + bool const doing_derivs = false); + +.. code-block:: cpp + + template + PORTABLE_INLINE_FUNCTION void + GetStatesFromPressureTemperature(Real &density_tot, Real &sie_tot, const Real pressure, + const Real temperature, RealIndexer &density_mat, + RealIndexer &sie_mat, LambdaIndexer &lambdas); + +Note that the ``doing_derivs`` option essentially tightens the tolerance on the +PTE solver in order to ensure that a perturbation to the PTE state via a finite +difference leads to a new PTE solution. The ``small_mass_mat_consistency`` +option for the density-energy lookup populates materials below the mass +fraction cutoff with P-T lookups from the PTE state for the other materials. + +Similarly, there are also public member functions for calculating the mixture +derivatives given a PTE state. For the heat capacity function, + +.. code-block:: cpp + + template + PORTABLE_INLINE_FUNCTION Real CalculateCvFromState(Real const rho, Real const sie, + [[maybe_unused]] Real const pressure, + Real const temperature, + RealIndexer &density_mat, + RealIndexer &sie_mat, + LambdaIndexer &lambdas); + +the state is perturbed in temperature space at constant density. For the bulk +modulus function, + +.. code-block:: cpp + + template + PORTABLE_INLINE_FUNCTION Real CalculateBmodFromState( + Real const rho, Real const sie, Real const pressure, Real const temperature, + RealIndexer &density_mat, RealIndexer &sie_mat, LambdaIndexer &lambdas); + +the state is perturbed in both temperature and density (with all materials in +PTE) where the bulk modulus is found from the formula, + +.. math:: + + B_S = B_T + \frac{T}{\rho} \left( \frac{\partial P}{\partial T} \right)_\rho^2 + \left( \frac{\partial T}{\partial e} \right)_\rho. + +The Gruneisen parameter function is + +.. code-block:: cpp + + template + PORTABLE_INLINE_FUNCTION Real CalculateGruneisenFromState( + Real const rho, Real const sie, Real const pressure, + [[maybe_unused]] Real const temperature, RealIndexer &density_mat, + RealIndexer &sie_mat, LambdaIndexer &lambdas); + +and perturbs the PTE energy at constant density. + +Note that in the ``FillEOS()`` and ``ValuesAtReferenceState()`` functions, the +values are all caluclated from density-temperature perturbations and the +appropriate thermodynamic identities. This can lead to slight inconsistencies, +but fewer PTE perturbations are needed as a result. + +There is also a public member function to directly perturb the PTE state in +density and temperature, which is what the bulk modulus, ``FillEOS()``, and +``ValuesAtReferenceState()`` functions use: + +.. code-block:: cpp + + template + PORTABLE_INLINE_FUNCTION void + PerturbPTEStateRT(Real const rho, Real const sie, Real const pressure, + Real const temperature, RealIndexer const &density_mat, + RealIndexer const &sie_mat, LambdaIndexer &lambdas, Real &dedT_R, + Real &dedR_T, Real &dPdT_R, Real &dPdR_T) + -.. _EOSPAC: https://laws.lanl.gov/projects/data/eos/eospacReleases.php diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index 326c4173367..88f8b9cf784 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1016,8 +1016,6 @@ might look something like this: eos = EOSBuilder::Modify(eos, scale); } -.. _eos methods reference section: - CheckParams ------------ @@ -1030,6 +1028,8 @@ which raise an error and/or print an equation of state specific error message if something has gone wrong. Most EOS constructors and ways of building an EOS call ``CheckParams`` by default. +.. _eos methods reference section: + Equation of State Methods Reference ------------------------------------ diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 81302eb703c..66cd19632d1 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -25,6 +25,7 @@ register_headers( # Normal files base/fast-math/logs.hpp + base/generic_indexer.hpp base/indexable_types.hpp base/robust_utils.hpp base/root-finding-1d/root_finding.hpp @@ -38,6 +39,7 @@ register_headers( base/sp5/singularity_eos_sp5.hpp eos/default_variant.hpp base/hermite.hpp + base/tuple_utils.hpp eos/eos_variant.hpp eos/eos_stellar_collapse.hpp eos/eos_ideal.hpp @@ -72,6 +74,7 @@ register_headers( eos/variant_utils.hpp eos/eos_carnahan_starling.hpp eos/eos_electrons.hpp + eos/eos_multi_eos.hpp ) if (SINGULARITY_BUILD_CLOSURE) diff --git a/singularity-eos/base/generic_indexer.hpp b/singularity-eos/base/generic_indexer.hpp new file mode 100644 index 00000000000..70c35a2ea79 --- /dev/null +++ b/singularity-eos/base/generic_indexer.hpp @@ -0,0 +1,62 @@ +//------------------------------------------------------------------------------ +// © 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_BASE_GENERIC_INDEXER_HPP_ +#define SINGULARITY_EOS_BASE_GENERIC_INDEXER_HPP_ + +#include + +namespace singularity { + +template +struct GenericIndexer { + Arr arr_; + Map map_; + + template + constexpr GenericIndexer(ArrT_ &&arr_in, MapT_ &&map_in) + : arr_(std::forward(arr_in)), map_(std::forward(map_in)) {} + + // & : non-const lvalue + template + constexpr decltype(auto) operator[](I i) & { + return arr_[map_[i]]; + } + + // const& : const lvalue + template + constexpr decltype(auto) operator[](I i) const & { + return arr_[map_[i]]; + } + + // && : rvalue (indexer is a temporary) — forward arr_’s value category + template + constexpr decltype(auto) operator[](I i) && { + return std::forward(arr_)[map_[i]]; // move only if Arr is a value type + } + + // const rvalue indexer + template + constexpr decltype(auto) operator[](I i) const && { + return std::forward(arr_)[map_[i]]; // preserves const + } +}; + +// CTAD: preserve references for lvalues, decay for rvalues +template +GenericIndexer(ArrT_ &&, MapT_ &&) -> GenericIndexer; + +} // namespace singularity + +#endif // #ifndef SINGULARITY_EOS_BASE_GENERIC_INDEXER_HPP_ diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp index 9f6d4737f0e..78143196b49 100644 --- a/singularity-eos/base/indexable_types.hpp +++ b/singularity-eos/base/indexable_types.hpp @@ -23,6 +23,21 @@ #include #include +// SFINAE helper macro that checks if a given indexer object has the requested +// indexable type. +// TODO: The error is a simple template substitution failure error, but there +// are ways short of C++20 concepts to add more readability to the error. We may +// want to add this at some point +#define SINGULARITY_INDEXER_HAS_INDEXABLE_TYPE(Indexer, IndexableType) \ + typename = std::enable_if_t< \ + singularity::variadic_utils::is_indexable_v> + +// Just a small wrapper for the above macro specifically for mass fractions +// NOTE: it's assumed that matnum is a count, not an index (hence the minus 1) +#define SINGULARITY_INDEXER_HAS_MASS_FRAC(Indexer, matnum) \ + SINGULARITY_INDEXER_HAS_INDEXABLE_TYPE( \ + Indexer, singularity::IndexableTypes::MassFraction) + namespace singularity { namespace IndexerUtils { @@ -222,28 +237,27 @@ class VariadicIndexerBase { VariadicIndexerBase() = default; - PORTABLE_FORCEINLINE_FUNCTION - VariadicIndexerBase(const Data_t &data) : data_(data) {} + constexpr VariadicIndexerBase(const Data_t &data) : data_(data) {} template ::value>> - PORTABLE_FORCEINLINE_FUNCTION Real &operator[](const T &t) { + constexpr Real &operator[](const T &t) noexcept { constexpr std::size_t idx = variadic_utils::GetIndexInTL(); return data_[idx]; } - PORTABLE_FORCEINLINE_FUNCTION - Real &operator[](const std::size_t idx) { return data_[idx]; } + constexpr Real &operator[](const std::size_t idx) { return data_[idx]; } template ::value>> - PORTABLE_FORCEINLINE_FUNCTION const Real &operator[](const T &t) const { + constexpr const Real &operator[](const T &t) const noexcept { constexpr std::size_t idx = variadic_utils::GetIndexInTL(); return data_[idx]; } - PORTABLE_FORCEINLINE_FUNCTION - const Real &operator[](const std::size_t idx) const { return data_[idx]; } + constexpr const Real &operator[](const std::size_t idx) const noexcept { + return data_[idx]; + } static inline constexpr std::size_t size() { return sizeof...(Ts); } @@ -268,6 +282,8 @@ struct MeanAtomicNumber {}; struct ElectronFraction {}; struct RootStatus {}; struct TableStatus {}; +template +struct MassFraction {}; } // namespace IndexableTypes } // namespace singularity #endif // SINGULARITY_EOS_BASE_INDEXABLE_TYPES_ diff --git a/singularity-eos/base/math_utils.hpp b/singularity-eos/base/math_utils.hpp index f24fe4dbe81..5993043ba3b 100644 --- a/singularity-eos/base/math_utils.hpp +++ b/singularity-eos/base/math_utils.hpp @@ -15,6 +15,8 @@ #ifndef SINGULARITY_EOS_BASE_MATH_UTILS_HPP_ #define SINGULARITY_EOS_BASE_MATH_UTILS_HPP_ +#include + #include namespace singularity { @@ -58,6 +60,33 @@ PORTABLE_FORCEINLINE_FUNCTION auto pow10(const Real x) { return std::exp(ln10 * x); } +/* + * Kahan summation, with the Neumaier correction + * https://onlinelibrary.wiley.com/doi/10.1002/zamm.19740540106 + */ +struct IdentityOperator { + PORTABLE_FORCEINLINE_FUNCTION Real operator()(const Real x) const { return x; } +}; +template +PORTABLE_FORCEINLINE_FUNCTION Real +sum_neumaier(Data_t &&data, std::size_t n, std::size_t offset = 0, int iskip = -1, + const Operator_t &op = IdentityOperator()) { + Real sum = 0; + Real c = 0; // correction + for (std::size_t i = 0; i < n; ++i) { + if (iskip >= 0 && i == static_cast(iskip)) continue; + Real x = op(data[i + offset]); + Real t = sum + x; + if (std::abs(sum) >= std::abs(x)) { + c += (sum - t) + x; + } else { + c += (x - t) + sum; + } + sum = t; + } + return sum + c; +} + } // namespace math_utils } // namespace singularity diff --git a/singularity-eos/base/tuple_utils.hpp b/singularity-eos/base/tuple_utils.hpp new file mode 100644 index 00000000000..e75a5d8d14b --- /dev/null +++ b/singularity-eos/base/tuple_utils.hpp @@ -0,0 +1,55 @@ +//------------------------------------------------------------------------------ +// © 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_BASE_TUPLE_UTILS_HPP_ +#define SINGULARITY_EOS_BASE_TUPLE_UTILS_HPP_ + +#include + +namespace singularity { +namespace tuple_utils { + +namespace impl { + +// Helper that copies a tuple into a new tuple while also applying a +// transformation function to each element. Wrapped in `impl` namespace due to +// index sequence and called below with a more normal function signature +template +constexpr auto tuple_transform_iteration(TupleT &&tup, FunctionT &&function, + std::index_sequence) { + return std::make_tuple(function(std::get(std::forward(tup)))...); +} +} // namespace impl + +// Returns a copy of the tuple with a transform function called on each element +template +constexpr auto tuple_transform(TupleT &&tup, FunctionT &&function) { + constexpr std::size_t N = std::tuple_size_v>; + return impl::tuple_transform_iteration(std::forward(tup), + std::forward(function), + std::make_index_sequence{}); +} + +// SFINAE helper to dicriminate between packs and tuples +template +struct is_std_tuple : std::false_type {}; +template +struct is_std_tuple> : std::true_type {}; +template +inline constexpr bool is_std_tuple_v = is_std_tuple>::value; + +} // namespace tuple_utils +} // namespace singularity + +#endif // #ifndef SINGULARITY_EOS_BASE_TUPLE_UTILS_HPP_ diff --git a/singularity-eos/base/variadic_utils.hpp b/singularity-eos/base/variadic_utils.hpp index 8b30c0cb5b2..4908e095d88 100644 --- a/singularity-eos/base/variadic_utils.hpp +++ b/singularity-eos/base/variadic_utils.hpp @@ -92,6 +92,12 @@ constexpr bool contains_v() { template struct type_list {}; +// contains over a type_list +template +struct contains_list; +template +struct contains_list> : contains {}; + // variadic list of modifiers template