diff --git a/.kessel/workflows/cmake.py b/.kessel/workflows/cmake.py index dead5310b8d..f5b52a8d78f 100644 --- a/.kessel/workflows/cmake.py +++ b/.kessel/workflows/cmake.py @@ -61,6 +61,7 @@ def test(self, args): self.exec(f""" pushd {self.build_dir} if [[ -f sesame2spiner/sesame2spiner ]]; then + echo "/usr/projects/data/eos/eos-developmental/Sn2162/v01/sn2162-v01.bin" > sesameFilesDir.txt sesame2spiner/sesame2spiner -s materials.sp5 ../sesame2spiner/examples/unit_tests/*.dat sesame2spiner/sesame2spiner -s duplicates.sp5 ../sesame2spiner/examples/duplicate-test/*.dat fi diff --git a/CHANGELOG.md b/CHANGELOG.md index c318bdf681e..63e5dbf91fb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR600]](https://github.com/lanl/singularity-eos/pull/600) Added mass fractions for multiphase tables to sesame2spiner and added MassFractionsFromDensityTemperture/InternalEnergy to the SpinerDependsRhoSIE and SpinerDependsRhoT equations of state. - [[PR589]](https://github.com/lanl/singularity-eos/pull/589) InternalEnergyFromDensityPressure ### Fixed (Repair bugs, etc) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 948abc1fd05..d2c56b764e3 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1720,6 +1720,39 @@ constructor for ``SpinerEOSDependsRhoSie`` is identical. ``SpinerEOS`` model does **not** support the ``MeanAtomicProperties`` struct. +Additionally Spiner EOS models support mass fraction lookups of the form + +.. code-block:: cpp + + template + PORTABLE_INLINE_FUNCTION + void SpinerEOSDependsRhoT::MassFractionsFromDensityTemperature( + const Real rho, const Real temp, Real *mass_fracs, Indexer_t &&lambda) const; + + template + PORTABLE_INLINE_FUNCTION void + SpinerEOSDependsRhoT::MassFractionsFromDensityInternalEnergy( + const Real rho, const Real sie, Real *mass_fracs, Indexer_t &&lambda) const; + +which sets the members of the ``*mass_fracs`` array the mass fractions of multiple phases as a function of density and either temperature or internal energy in a table if they are available. If they are not available, the ``mass_fracs`` array is assumed to be of at least length 1 and is set to unity. Alternatively, the functions + +.. code-block:: cpp + + template + PORTABLE_INLINE_FUNCTION + void SpinerEOSDependsRhoT::MassFractionsFromDensityTemperature( + const Real rho, const Real temp, Indexer_t &&lambda) const; + + template + PORTABLE_INLINE_FUNCTION + void SpinerEOSDependsRhoT::MassFractionsFromDensityInternalEnergy( + const Real rho, const Real sie, Indexer_t &&lambda) const; + +assume the mass fractions array is contained in the lambda indexer and +set them accordingly. They are assumed to be indexable via the +``MassFractions`` indexable type, discussed above, or to be appended +at the end of the lambda after ``nlambda()`` elements. + ``sp5`` files and ``sesame2spiner`` ````````````````````````````````````` diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index 7105f884852..faae9c7fbe7 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -694,7 +694,29 @@ The available types currently supported by default are: struct singularity::IndexableTypes::MeanAtomicNumber; struct singularity::IndexableTypes::ElectronFraction; -However if you are not limited to these types. Any type will do and +Additionally, indexable types may be "vector" quantities. The type + +.. code-block:: cpp + + struct singularity::IndexableTypes::MassFractions; + +must be constructed with an index, e.g., + +.. code-block:: cpp + + eos.MassFractionsFromDensityTemperature(rho, T, lambda); + x0 = lambda[MassFractions(0)]; + +and the default constructor is not supported. Finally, the structs + +.. code-block:: cpp + + struct singularity::RootStatus; + struct singularity::TableStatus; + +are mostly used for internal/debugging purposes. + +However you are not limited to these types. Any type will do and you can define your own as you like. For example: .. code-block:: @@ -736,17 +758,48 @@ which might be used as where ``MeanIonizationState`` is shorthand for index 2, since you defined that overload. Note that the ``operator[]`` must be marked -``const``. To more easily enable mixing and matching integer-based +``const``. "Vectorized" types such as ``MassFractions`` are expected +to expose a public member field named ``n``, which can be utilized by a +custom indexer class. For example: + +.. code-block:: cpp + + class MyLambda_t { + public: + constexpr static bool is_type_indexable = true; + MyLambda_t() = default; + PORTABLE_FORCEINLINE_FUNCTION + Real &operator[](const MeanIonizationState &zbar) const { + return data_[0]; + } + // x.n used as on offset within the underlying container + PORTABLE_FORCEINLINE_FUNCTION + Real &operator[](const MassFractions &x) const { + return data_[1 + x.n]; + } + private: + std::array data_; + }; + +To more easily enable mixing and matching integer-based indexing with type-based indexing, the functions .. code-block:: cpp + template + PORTABLE_FORCEINLINE_FUNCTION + bool SafeGet(Indexer_t const &lambda, const T &t, std::size_t const idx, Real &out); + template PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t const idx, Real &out); .. code-block:: cpp + template + PORTABLE_FORCEINLINE_FUNCTION + bool SafeGet(Indexer_t const &lambda, const T &t, Real &out); + template PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, Real &out); @@ -755,13 +808,21 @@ indexing with type-based indexing, the functions template PORTABLE_FORCEINLINE_FUNCTION - Real SafeMustGet(Indexer_t const &lambda, std::size_t const idx) + Real SafeMustGet(Indexer_t const &lambda, const T &t, std::size_t const idx); + + template + PORTABLE_FORCEINLINE_FUNCTION + Real SafeMustGet(Indexer_t const &lambda, std::size_t const idx); .. code-block:: cpp template PORTABLE_FORCEINLINE_FUNCTION - Real SafeMustGet(Indexer_t const &lambda) + Real SafeMustGet(Indexer_t const &lambda, const T &t); + + template + PORTABLE_FORCEINLINE_FUNCTION + Real SafeMustGet(Indexer_t const &lambda); will update the value of ``out`` with the value at either the appropriate type-based index, ``T``, or the numerical index, ``idx``, if the ``Indexer_t`` @@ -788,24 +849,40 @@ Similarly, the functions .. code-block:: cpp + template + PORTABLE_FORCEINLINE_FUNCTION + inline bool SafeSet(Indexer_t &lambda, const T &t, std::size_t const idx, Real const in); + template PORTABLE_FORCEINLINE_FUNCTION inline bool SafeSet(Indexer_t &lambda, std::size_t const idx, Real const in); .. code-block:: cpp + template + PORTABLE_FORCEINLINE_FUNCTION + inline bool SafeSet(Indexer_t &lambda, const T &t, Real const in) + template PORTABLE_FORCEINLINE_FUNCTION inline bool SafeSet(Indexer_t &lambda, Real const in) .. code-block:: cpp + template + PORTABLE_FORCEINLINE_FUNCTION + inline bool SafeMustSet(Indexer_t &lambda, const T &t, std::size_t const idx, Real const in); + template PORTABLE_FORCEINLINE_FUNCTION inline bool SafeMustSet(Indexer_t &lambda, std::size_t const idx, Real const in); .. code-block:: cpp + template + PORTABLE_FORCEINLINE_FUNCTION + inline bool SafeMustSet(Indexer_t &lambda, const T &t, Real const in) + template PORTABLE_FORCEINLINE_FUNCTION inline bool SafeMustSet(Indexer_t &lambda, Real const in) diff --git a/eospac-wrapper/eospac_wrapper.cpp b/eospac-wrapper/eospac_wrapper.cpp index daaa0f3886e..8a5e64b652d 100644 --- a/eospac-wrapper/eospac_wrapper.cpp +++ b/eospac-wrapper/eospac_wrapper.cpp @@ -335,6 +335,14 @@ void eosSafeTableCmnts(EOS_INTEGER *table, EOS_CHAR *comments, Verbosity eospacW eosCheckError(errorCode, "eos_GetTableCmnts", eospacWarn); } +void eosSafeTableMetaData(EOS_INTEGER *table, EOS_INTEGER infoItem, EOS_CHAR *infoString, + Verbosity eospacWarn) { + EOS_INTEGER errorCode = EOS_OK; + EOS_INTEGER infoTypes[1] = {infoItem}; + eos_GetTableMetaData(table, infoTypes, infoString, &errorCode); + eosCheckError(errorCode, "eos_GetTableMetaData", eospacWarn); +} + void eosCheckError(EOS_INTEGER errorCode, const std::string &name, Verbosity eospacWarn) { EOS_CHAR errorMessage[EOS_MaxErrMsgLen]; if (errorCode != EOS_OK && eospacWarn != Verbosity::Quiet) { @@ -343,6 +351,31 @@ void eosCheckError(EOS_INTEGER errorCode, const std::string &name, Verbosity eos } } +int eosCheckTableExistence(EOS_INTEGER tableType_, int matid_, Verbosity eospacWarn) { + // Note this opens and closes the table + int exists = 1; + int NT = 1; + EOS_INTEGER matid[1] = {matid_}; + EOS_INTEGER tableType[1] = {tableType_}; + EOS_INTEGER tableHandle[1]; + EOS_INTEGER errorCode = EOS_OK; + + eos_CreateTables(&NT, tableType, matid, tableHandle, &errorCode); + eosCheckError(errorCode, "eos_CreateTables", eospacWarn); + eos_LoadTables(&NT, &tableHandle[0], &errorCode); + if (errorCode != EOS_OK) { + if (eosErrorCodesEqual(EOS_NO_DATA_TABLE, errorCode)) { + exists = 0; + } else { + // something else happended + eosCheckError(errorCode, "eos_LoadTables", eospacWarn); + exists = -1; + } + } + eosSafeDestroy(NT, tableHandle, eospacWarn); + return exists; +} + std::string eosErrorString(EOS_INTEGER errorCode) { // Algorithmicallly generated by parsing the EOSPAC docs // I'm sorry. It's gross. ~JMM diff --git a/eospac-wrapper/eospac_wrapper.hpp b/eospac-wrapper/eospac_wrapper.hpp index 45079f47159..637fe5b0170 100644 --- a/eospac-wrapper/eospac_wrapper.hpp +++ b/eospac-wrapper/eospac_wrapper.hpp @@ -122,13 +122,14 @@ void eosSafeTableInfo(EOS_INTEGER *table, EOS_INTEGER numInfoItems, EOS_INTEGER infoItems[], EOS_REAL infoVals[], Verbosity eospacWarn); void eosSafeTableCmnts(EOS_INTEGER *table, EOS_CHAR *comments, Verbosity eospacWarn); - +void eosSafeTableMetaData(EOS_INTEGER *table, EOS_INTEGER infoItem, EOS_CHAR *infoString, + Verbosity eospacWarn); void eosCheckError(EOS_INTEGER errorCode, const std::string &name, Verbosity eospacWarn); std::string eosErrorString(EOS_INTEGER errorCode); void eosSafeDestroy(int ntables, EOS_INTEGER tableHandles[], Verbosity eospacWarn); std::string getName(std::string comment); bool eosErrorCodesEqual(EOS_INTEGER lhs, EOS_INTEGER rhs); - +int eosCheckTableExistence(EOS_INTEGER tableType_, int matid_, Verbosity eospacWarn); } // namespace EospacWrapper #endif // _EOSPAC_WRAPPER_EOSPAC_WRAPPER_HPP_ diff --git a/sesame2spiner/examples/unit_tests/tin.dat b/sesame2spiner/examples/unit_tests/tin.dat new file mode 100644 index 00000000000..523f138b48b --- /dev/null +++ b/sesame2spiner/examples/unit_tests/tin.dat @@ -0,0 +1,18 @@ +# ====================================================================== +# Example input deck for sesame2spiner, +# a tool for converting eospac to spiner +# © 2026. 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. +# ====================================================================== + +matid=2162 +name=tin diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index eeafb752936..c3b5d48ae52 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -53,7 +53,7 @@ herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRh std::string sMatid = std::to_string(matid); herr_t status = 0; - hid_t matGroup, lTGroup, leGroup, coldGroup; + hid_t matGroup, lTGroup, leGroup, coldGroup, mfGroup; matGroup = H5Gcreate(loc, sMatid.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); status += H5Lcreate_soft(sMatid.c_str(), loc, name.c_str(), H5P_DEFAULT, H5P_DEFAULT); @@ -116,6 +116,25 @@ herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRh // status += transitionMask.saveHDF(coldGroup, SP5::Fields::transitionMask); } + { + DataBox mf, mask; + Bounds nphBounds; + std::string phase_names; + + if (eosMassFraction(matid, lRhoBounds, lTBounds, nphBounds, mf, mask, phase_names, + eospacWarn)) { + mfGroup = H5Gcreate(matGroup, SP5::Depends::massFrac, H5P_DEFAULT, H5P_DEFAULT, + H5P_DEFAULT); + + status += mf.saveHDF(mfGroup, SP5::Fields::massFrac); + const int numphases = mf.dim(1); + status += H5LTset_attribute_int(mfGroup, ".", "numphases", &numphases, 1); + status += + H5LTset_attribute_string(mfGroup, ".", "phase names", phase_names.c_str()); + status += H5Gclose(mfGroup); + } + } + if (addSubtables) { int i = 0; std::vector splits = {TableSplit::ElectronOnly, TableSplit::IonCold}; diff --git a/sesame2spiner/io_eospac.cpp b/sesame2spiner/io_eospac.cpp index 2eac41d536a..19fb3e81ffd 100644 --- a/sesame2spiner/io_eospac.cpp +++ b/sesame2spiner/io_eospac.cpp @@ -328,6 +328,82 @@ void eosColdCurveMask(int matid, const Bounds &lRhoBounds, const int numSie, eosSafeDestroy(NT, tableHandle, eospacWarn); } +bool eosMassFraction(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds, + Bounds &nphBounds, DataBox &Ms, DataBox &mask, + std::string &phase_names, Verbosity eospacWarn) { + using namespace EospacWrapper; + + EOS_INTEGER errorCode = EOS_OK; + constexpr int NT = 2; + EOS_INTEGER tableHandle[NT]; + EOS_INTEGER tableType[NT] = {EOS_M_DT, EOS_Comment}; + std::vector matid_v(NT, matid); + std::vector table_names = {"EOS_M_DT", "Eos_Material_Phases"}; + + const int exists = eosCheckTableExistence(EOS_M_DT, matid, eospacWarn); + if (exists > 0) { + eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_M_DT", "Eos_Material_Phases"}, + eospacWarn); + } else { + phase_names = std::string(""); + return false; + } + // indep vars + // Reuses the EOS log temp and log rho bounds + // TODO(@adempsey): Use a different grid for mass fractions: + // - Linear in rho and T + // - Windowed grids for each phase + std::vector rhos, Ts, phs; + makeInterpPoints(rhos, lRhoBounds); + makeInterpPoints(Ts, lTBounds); + + EOS_REAL infoVals[1]; + EOS_INTEGER infoItems[1] = {EOS_NUM_PHASES}; + eosSafeTableInfo(&tableHandle[0], 1, infoItems, infoVals, eospacWarn); + const int nph = static_cast(infoVals[0]); + + EOS_CHAR infoString[EOS_META_DATA_STRLEN]; + infoString[0] = '\0'; + eosSafeTableMetaData(&tableHandle[1], EOS_Material_Phases, infoString, eospacWarn); + phase_names = std::string(infoString); + + DataBox dMdt, dMdr; + Ms.resize(rhos.size(), Ts.size(), nph); + Ms.setRange(1, lTBounds.grid); + Ms.setRange(2, lRhoBounds.grid); + + dMdr.copyMetadata(Ms); + dMdt.copyMetadata(Ms); + mask.copyMetadata(Ms); + + // Interpolatable vars + EOS_INTEGER nXYPairs = rhos.size() * Ts.size(); + std::vector M_pack(nXYPairs * nph), DMDR_T(nXYPairs), DMDT_R(nXYPairs), + rho_flat(nXYPairs), T_flat(nXYPairs); + + // prepare flat data structures + for (std::size_t j = 0, iflat = 0; j < rhos.size(); ++j) { + for (std::size_t i = 0; i < Ts.size(); ++i, iflat++) { + rho_flat[iflat] = densityToSesame(rhos[j]); + T_flat[iflat] = temperatureToSesame(Ts[i]); + } + } + + const bool no_errors = eosSafeInterpolate(&tableHandle[0], nXYPairs, rho_flat.data(), + T_flat.data(), M_pack.data(), DMDR_T.data(), + DMDT_R.data(), "EOS_M_DT", eospacWarn); + + for (size_t j = 0, iflat = 0; j < rhos.size(); j++) { + for (size_t i = 0; i < Ts.size(); i++, iflat++) { + for (size_t k = 0; k < nph; k++) { + Ms(j, i, k) = M_pack[k * nXYPairs + iflat]; + } + } + } + eosSafeDestroy(NT, tableHandle, eospacWarn); + return true; +} + void makeInterpPoints(std::vector &v, const Bounds &b) { v.resize(b.grid.nPoints()); for (size_t i = 0; i < v.size(); i++) { diff --git a/sesame2spiner/io_eospac.hpp b/sesame2spiner/io_eospac.hpp index cd73390cb54..2a2c166da2d 100644 --- a/sesame2spiner/io_eospac.hpp +++ b/sesame2spiner/io_eospac.hpp @@ -76,6 +76,9 @@ void eosColdCurves(int matid, const Bounds &lRhoBounds, DataBox &Ps, DataBox &si void eosColdCurveMask(int matid, const Bounds &lRhoBounds, const int numSie, const DataBox &sieColdCurve, DataBox &mask, Verbosity eospacWarn = Verbosity::Quiet); +bool eosMassFraction(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds, + Bounds &nphBounds, DataBox &mf, DataBox &mask, + std::string &phase_names, Verbosity eospacWarn = Verbosity::Quiet); void makeInterpPoints(std::vector &v, const Bounds &b); diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp index 9f6d4737f0e..6e188fa5fe9 100644 --- a/singularity-eos/base/indexable_types.hpp +++ b/singularity-eos/base/indexable_types.hpp @@ -45,8 +45,8 @@ enum class AllowedIndexing { Numeric, TypeOnly }; // type-based index is present in the Indexer OR if the Indexer doesn't support // type-based indexing. template -PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t const idx, - Real &out) { +PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, const T &t, + std::size_t const idx, Real &out) { // If null then nothing happens if (variadic_utils::is_nullptr(lambda)) { return false; @@ -54,7 +54,7 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t // Return value if type index is available if constexpr (variadic_utils::is_indexable_v) { - out = lambda[T{}]; + out = lambda[t]; return true; } @@ -74,13 +74,18 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t // Something else... return false; } +template +PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t const idx, + Real &out) { + return SafeGet(lambda, T{}, idx, out); +} // Break out "Set" functionality from "Get". The original "Get()" did both, but // the "safe" version needs to separate that functionality for setting the // values in a lambda template -PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const idx, - Real const in) { +PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, const T &t, + std::size_t const idx, Real const in) { // If null then nothing happens if (variadic_utils::is_nullptr(lambda)) { return false; @@ -88,7 +93,7 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const // Return value if type index is available if constexpr (variadic_utils::is_indexable_v) { - lambda[T{}] = in; + lambda[t] = in; return true; } @@ -108,6 +113,11 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const // Something else... return false; } +template +PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const idx, + Real const in) { + return SafeSet(lambda, T{}, idx, in); +} // Same as above but causes an error condition (static or dynamic) if the value // can't be obtained. Note that the `decltype(auto)` is intended to preserve the @@ -117,8 +127,8 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const // for either setting or getting values. This should also allow for `const` // correctness downstream in the wrappers where `lambda` is `const` template -PORTABLE_FORCEINLINE_FUNCTION decltype(auto) SafeMustGetSet(Indexer_t &&lambda, - std::size_t const idx) { +PORTABLE_FORCEINLINE_FUNCTION decltype(auto) +SafeMustGetSet(Indexer_t &&lambda, const T &t, std::size_t const idx) { // Error on null pointer PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda), "Indexer can't be nullptr"); @@ -128,7 +138,7 @@ PORTABLE_FORCEINLINE_FUNCTION decltype(auto) SafeMustGetSet(Indexer_t &&lambda, static_assert(variadic_utils::is_indexable_v); // Use std::forward to maintain value category for lambda, and use // parentheses to do the same for the output of the lambda[] operation - return (std::forward(lambda)[T{}]); + return (std::forward(lambda)[t]); } else if constexpr (AI == AllowedIndexing::Numeric) { // Fall back to numerical indexing if allowed static_assert(variadic_utils::has_int_index_v); @@ -143,7 +153,11 @@ PORTABLE_FORCEINLINE_FUNCTION decltype(auto) SafeMustGetSet(Indexer_t &&lambda, "must be called with a numerical index"); } } - +template +PORTABLE_FORCEINLINE_FUNCTION decltype(auto) SafeMustGetSet(Indexer_t &&lambda, + std::size_t const idx) { + return SafeMustGetSet(lambda, T{}, idx); +} } // namespace impl // Overload when numerical index is provided @@ -152,6 +166,11 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t Real &out) { return impl::SafeGet(lambda, idx, out); } +template +PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, const T &t, + std::size_t const idx, Real &out) { + return impl::SafeGet(lambda, t, idx, out); +} // Overload when numerical index isn't provided template @@ -159,6 +178,12 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, Real &out) { std::size_t idx = 0; return impl::SafeGet(lambda, idx, out); } +template +PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, const T &t, + Real &out) { + std::size_t idx = 0; + return impl::SafeGet(lambda, t, idx, out); +} // Overload when numerical index is provided template @@ -166,6 +191,11 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const Real const in) { return impl::SafeSet(lambda, idx, in); } +template +PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, const T &t, + std::size_t const idx, Real const in) { + return impl::SafeSet(lambda, t, idx, in); +} // Overload when numerical index isn't provided template @@ -173,6 +203,11 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, Real const in) { std::size_t idx = 0; return impl::SafeSet(lambda, idx, in); } +template +PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, const T &t, Real const in) { + std::size_t idx = 0; + return impl::SafeSet(lambda, t, idx, in); +} // Overload when numerical index is provided template @@ -180,6 +215,11 @@ PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda, std::size_t const idx) { return impl::SafeMustGetSet(lambda, idx); } +template +PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda, const T &t, + std::size_t const idx) { + return impl::SafeMustGetSet(lambda, t, idx); +} // Overload when numerical index isn't provided template @@ -187,6 +227,11 @@ PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda) { std::size_t idx = 0; return impl::SafeMustGetSet(lambda, idx); } +template +PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda, const T &t) { + std::size_t idx = 0; + return impl::SafeMustGetSet(lambda, t, idx); +} // Overload when numerical index is provided template @@ -194,6 +239,11 @@ PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t co Real const in) { impl::SafeMustGetSet(lambda, idx) = in; } +template +PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, const T &t, + std::size_t const idx, Real const in) { + impl::SafeMustGetSet(lambda, t, idx) = in; +} // Overload when numerical index isn't provided template @@ -201,6 +251,12 @@ PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, Real const in) std::size_t idx = 0; impl::SafeMustGetSet(lambda, idx) = in; } +template +PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, const T &t, + Real const in) { + std::size_t idx = 0; + impl::SafeMustGetSet(lambda, t, idx) = in; +} // This is a convenience struct to easily build a small indexer with // a set of indexable types. @@ -268,6 +324,11 @@ struct MeanAtomicNumber {}; struct ElectronFraction {}; struct RootStatus {}; struct TableStatus {}; +struct MassFractions { + std::size_t n; + PORTABLE_FORCEINLINE_FUNCTION + MassFractions(const std::size_t index_) : n(index_) {} +}; } // namespace IndexableTypes } // namespace singularity #endif // SINGULARITY_EOS_BASE_INDEXABLE_TYPES_ diff --git a/singularity-eos/base/sp5/singularity_eos_sp5.hpp b/singularity-eos/base/sp5/singularity_eos_sp5.hpp index 002da520eac..67eedf64da9 100644 --- a/singularity-eos/base/sp5/singularity_eos_sp5.hpp +++ b/singularity-eos/base/sp5/singularity_eos_sp5.hpp @@ -26,6 +26,7 @@ namespace Depends { constexpr char logRhoLogSie[] = "dependsLogRhoLogSie"; constexpr char logRhoLogT[] = "dependsLogRhoLogT"; constexpr char coldCurve[] = "coldCurve"; +constexpr char massFrac[] = "massFrac"; } // namespace Depends namespace SubTable { @@ -68,6 +69,7 @@ constexpr char dEdRho[] = "dEdRho"; constexpr char dEdT[] = "dEdT"; constexpr char mask[] = "mask"; constexpr char transitionMask[] = "transition mask"; +constexpr char massFrac[] = "mass fractions"; } // namespace Fields } // namespace SP5 diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index 246144083ed..82f2f09e880 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -128,6 +128,7 @@ inline void SetUpOutputScalingOption(EOS_INTEGER options[], EOS_REAL values[], } // namespace impl_eospac +// TODO(@adempsey): Add mass fractions + better multiphase support class EOSPAC : public EosBase { public: diff --git a/singularity-eos/eos/eos_spiner_common.hpp b/singularity-eos/eos/eos_spiner_common.hpp index df04ff9d3cd..adec55d5164 100644 --- a/singularity-eos/eos/eos_spiner_common.hpp +++ b/singularity-eos/eos/eos_spiner_common.hpp @@ -151,6 +151,38 @@ inline hid_t h5_safe_gopen(hid_t &file, const char *name, hid_t property) { return H5Gopen(file, name, property); } +inline char *h5_safe_read_attr_string(hid_t &file, const char *grp, const char *name, + std::size_t &len) { + // NOTE THIS FUNCTION ALLOCATES THE BUFFER STRING!!! + hsize_t dims[1] = {0}; + H5T_class_t attr_type; + if (H5LTget_attribute_info(file, grp, name, dims, &attr_type, &len) != H5_SUCCESS) { + std::string msg = "Could not read the string " + std::string(name); + EOS_ERROR(msg); + return nullptr; + } + if (len <= 0) { + std::string msg = "Length of string" + std::string(name) + " is 0"; + EOS_ERROR(msg); + return nullptr; + } + if (attr_type != H5T_STRING) { + std::string msg = "Attribute type is NOT a string for " + std::string(name); + EOS_ERROR(msg); + return nullptr; + } + len += 1; + char *buffer = (char *)malloc(len); + if (H5LTget_attribute_string(file, grp, name, buffer) != H5_SUCCESS) { + std::string msg = "Error reading string attribute for " + std::string(name); + free(buffer); + EOS_ERROR(msg); + return nullptr; + } + buffer[len - 1] = '\0'; + return buffer; +} + template inline void h5_safe_get_attribute(hid_t &file, const char *grp, const char *name, T *output, const bool optional = false) { diff --git a/singularity-eos/eos/eos_spiner_rho_sie.hpp b/singularity-eos/eos/eos_spiner_rho_sie.hpp index 1a4041608ee..f0c2986fa4d 100644 --- a/singularity-eos/eos/eos_spiner_rho_sie.hpp +++ b/singularity-eos/eos/eos_spiner_rho_sie.hpp @@ -196,6 +196,22 @@ class SpinerEOSDependsRhoSieTransformable const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const; template + PORTABLE_INLINE_FUNCTION void MassFractionsFromDensityTemperature( + const Real rho, const Real T, Real *mass_frac, + Indexer_t &&lambda = static_cast(nullptr)) const; + template + PORTABLE_INLINE_FUNCTION void + MassFractionsFromDensityTemperature(const Real rho, const Real T, + Indexer_t &&lambda) const; + template + PORTABLE_INLINE_FUNCTION void MassFractionsFromDensityInternalEnergy( + const Real rho, const Real sie, Real *mass_frac, + Indexer_t &&lambda = static_cast(nullptr)) const; + template + PORTABLE_INLINE_FUNCTION void + MassFractionsFromDensityInternalEnergy(const Real rho, const Real sie, + Indexer_t &&lambda) const; + template PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const; @@ -244,8 +260,13 @@ class SpinerEOSDependsRhoSieTransformable Real RhoPmin(const Real temp) const { return rho_at_pmin_.interpToReal(spiner_common::to_log(temp, lTOffset_)); } + PORTABLE_FORCEINLINE_FUNCTION int GetNumberofPhases() const { return numphases; } + const char *GetPhaseNames() const { return phase_names; } + // TODO(JMM): Should nlambda be made non-static so it can report the + // number of phases too? constexpr static inline int nlambda() noexcept { return _n_lambda; } + template static inline constexpr bool NeedsLambda() { return std::is_same::value; @@ -275,7 +296,7 @@ class SpinerEOSDependsRhoSieTransformable private: inline herr_t loadDataboxes_(const std::string &matid_str, hid_t file, hid_t lTGroup, - hid_t lEGroup, hid_t coldGroupd); + hid_t lEGroup, hid_t coldGroupd, hid_t mfGroup); inline void calcBMod_(SP5Tables &tables); template @@ -295,6 +316,7 @@ class SpinerEOSDependsRhoSieTransformable DataBox rho_at_pmin_; SP5Tables dependsRhoT_; SP5Tables dependsRhoSie_; + DataBox mF_; DataBox PlRhoMax_, dPdRhoMax_; DataBox PCold_, sieCold_, bModCold_, dPdRhoCold_; int numRho_, numT_; @@ -309,7 +331,7 @@ class SpinerEOSDependsRhoSieTransformable &(dependsRhoT_.dPdRho), &(dependsRhoT_.dPdE), &(dependsRhoT_.dTdRho), \ &(dependsRhoT_.dTdE), &(dependsRhoT_.dEdRho), &(dependsRhoSie_.P), \ &(dependsRhoSie_.bMod), &(dependsRhoSie_.dPdRho), &(dependsRhoSie_.dPdE), \ - &(dependsRhoSie_.dTdRho), &(dependsRhoSie_.dTdE), &(dependsRhoSie_.dEdRho), \ + &(dependsRhoSie_.dTdRho), &(dependsRhoSie_.dTdE), &(dependsRhoSie_.dEdRho), &mF_, \ &PlRhoMax_, &dPdRhoMax_, &PCold_, &sieCold_, &bModCold_, &dPdRhoCold_ std::vector GetDataBoxPointers_() const { return std::vector{DBLIST}; @@ -325,6 +347,14 @@ class SpinerEOSDependsRhoSieTransformable MeanAtomicProperties AZbar_; bool reproducible_ = false; bool pmin_vapor_dome_ = false; + bool has_mf = false; + // Need to hold the phase names for multiphase EOS + // This isn't great, but the class needs to be trivially copyable + // I've chosen something reasonable, e.g., 15 phases with 32 character names + char *phase_names = nullptr; + std::size_t len_phase_names = 0; + DataStatus phase_names_status = DataStatus::Deallocated; + int numphases = 1; // only used to exclude vapor dome static constexpr const Real VAPOR_DPDR_THRESH = 1e-8; static constexpr const int _n_lambda = 1; @@ -373,11 +403,20 @@ inline SpinerEOSDependsRhoSieTransformable< spiner_common::h5_safe_get_attribute(file, materialName.c_str(), SP5::Material::matid, &matid_); matid_str = std::to_string(matid_); + // mass fractions + has_mf = H5Lexists(matGroup, SP5::Depends::massFrac, H5P_DEFAULT); + hid_t mfGroup = -1; + if (has_mf) { + mfGroup = spiner_common::h5_safe_gopen(matGroup, SP5::Depends::massFrac, H5P_DEFAULT); + } - status += loadDataboxes_(matid_str, file, lTGroup, lEGroup, coldGroup); + status += loadDataboxes_(matid_str, file, lTGroup, lEGroup, coldGroup, mfGroup); InitializeTransformer(); + if (has_mf) { + spiner_common::h5_safe_gclose(mfGroup); + } spiner_common::h5_safe_gclose(lTGroup); spiner_common::h5_safe_gclose(lEGroup); spiner_common::h5_safe_gclose(matGroup); @@ -387,7 +426,7 @@ inline SpinerEOSDependsRhoSieTransformable< template