diff --git a/CHANGELOG.md b/CHANGELOG.md index 2434898253..c4140b03b4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR459]](https://github.com/lanl/singularity-eos/pull/459) Add electron and ion tables to EOSPAC and SpinerEOS backends - [[PR453]](https://github.com/lanl/singularity-eos/pull/453) A PT space PTE solver - [[PR444]](https://github.com/lanl/singularity-eos/pull/444) Add Z split modifier and electron ideal gas EOS diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 4fa44a29ef..28d4a6dbb8 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -218,9 +218,26 @@ density. ``singularity-eos`` assumes that, when 3T physics is active, electrons and ions are each described by a separate equation of state object. Several models are specifically designed to represent, e.g., -the electron equation of state or ion equation of state. The tabulated -models may also support loading tables specifically for electron or -ion equations of state. +the electron equation of state or ion equation of state. + +The tabulated models may also support loading tables specifically for +electron or ion equations of state. In these cases, an ``enum class`` +specifies which component of the material is being requst: + +.. code-block:: cpp + + enum class TableSplit { Total, ElectronOnly, IonCold }; + +where here ``ElectronOnly`` is the free electrons of the material, +corresponding to Sesame 304 tables, ``IonCold`` is the ions plus cold +curve (i.e., the ions if you don't know what a cold curve is), +corresponding to the Sesame 303 tables, and ``Total`` is the sum of +free electrons and ions, correspodning to the Sesame 301 tables. + +.. note:: + + The tabulated equations of state have an implicitly assumed + ionization fraction that depends on the electron temperature. Available EOS Information and Nomenclature ------------------------------------------ @@ -1665,15 +1682,19 @@ The constructor for ``SpinerEOSDependsRhoT`` is given by two overloads: .. code-block:: cpp SpinerEOSDependsRhoT(const std::string &filename, int matid, + const singularity::TableSplit = singularity::TableSplit::Total, bool reproduciblity_mode = false); SpinerEOSDependsRhoT(const std::string &filename, const std::string &materialName, + const singularity::TableSplit = singularity::TableSplit::Total, bool reproducibility_mode = false); where here ``filename`` is the input file, ``matid`` is the unique material ID in the database in the file, ``materialName`` is the name of the material in the file, and ``reproducability_mode`` is a boolean which slightly changes how initial guesses for root finds are -computed. The constructor for ``SpinerEOSDependsRhoSie`` is identical. +computed. The ``TableSplit`` option selects electron, ion, or total +equation of state tables for partial ionization. The constructor for +``SpinerEOSDependsRhoSie`` is identical. .. note:: @@ -1736,6 +1757,12 @@ and files to locate files in the `sesame`_ database, and database need not be provided by the command line. For how to specify `sesame`_ file locations, see the `eospac`_ manual. +.. note:: + + To enable 3T subtables with ``SpinerEOS``, you must set + ``ionization=true`` in the input file for ``sesame2spiner`` for the + desired material. + Piecewise Spiner Grids ```````````````````````` @@ -2134,27 +2161,40 @@ EOSPAC EOS 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 looks like +`Sesame`_ database. The constructor for the ``EOSPAC`` model has several overloads .. code-block:: - EOSPAC(int matid, bool invert_at_setup = false, Real insert_data = 0.0, eospacMonotonicity monotonicity = eospacMonotonicity::none, bool apply_smoothing = false, eospacSplit apply_splitting = eospacSplit::none, bool linear_interp = false) + EOSPAC(int matid, TableSplit split, bool invert_at_setup = false, + Real insert_data = 0.0, + eospacMonotonicity monotonicity = eospacMonotonicity::none, + bool apply_smoothing = false, + eospacSplit apply_splitting = eospacSplit::none, + bool linear_interp = false); + EOSPAC(int matid, bool invert_at_setup = false, Real insert_data = 0.0, + eospacMonotonicity monotonicity = eospacMonotonicity::none, + bool apply_smoothing = false, + eospacSplit apply_splitting = eospacSplit::none, + bool linear_interp = false); + EOSPAC(int matid, const Options &opts); where ``matid`` is the unique material number in the database, -``invert_at_setup`` specifies whether or not pre-compute tables of -temperature as a function of density and energy, ``insert_data`` -inserts specified number of grid points between original grid points -in the `Sesame`_ table, ``monotonicity` enforces monotonicity in x, -y or both (:math:`monotonicityX/Y/XY`), ``apply_smoothing`` enables -data table smoothing that imposes a linear floor on temperature dependence, -forces linear temperature dependence for low temperature, and forces -linear density dependence for low and high density, ``apply_splitting`` -has the following options for ion data tables not found in the `Sesame`_ -database :. :math:`splitNumProp` uses the cold curve plus number-proportional -model, :math:`splitIdealGas` uses the cold curve plus ideal gas model -and :math:`splitCowan` uses the cold curve plus Cowan-nuclear model -for ions and the final option ``linear_interp`` uses linear instead of -bilinear interpolation. +``split`` is the ``TableSplit`` parameter for electron, ion, or total +tables used in partial ionization. ``invert_at_setup`` specifies +whether or not pre-compute tables of temperature as a function of +density and energy, ``insert_data`` inserts specified number of grid +points between original grid points in the `Sesame`_ table, +``monotonicity` enforces monotonicity in x, y or both +(:math:`monotonicityX/Y/XY`), ``apply_smoothing`` enables data table +smoothing that imposes a linear floor on temperature dependence, +forces linear temperature dependence for low temperature, and forces +linear density dependence for low and high density, +``apply_splitting`` has the following options for ion data tables not +found in the `Sesame`_ database :. :math:`splitNumProp` uses the cold +curve plus number-proportional model, :math:`splitIdealGas` uses the +cold curve plus ideal gas model and :math:`splitCowan` uses the cold +curve plus Cowan-nuclear model for ions and the final option +``linear_interp`` uses linear instead of bilinear interpolation. .. note:: diff --git a/sesame2spiner/examples/unit_tests/steel.dat b/sesame2spiner/examples/unit_tests/steel.dat index cb8e4887d3..ebc8874636 100644 --- a/sesame2spiner/examples/unit_tests/steel.dat +++ b/sesame2spiner/examples/unit_tests/steel.dat @@ -2,7 +2,7 @@ # Example input deck for sesame2spiner, # a tool for converting eospac to spiner # Author: Jonah Miller (jonahm@lanl.gov) -# © 2021-2023. Triad National Security, LLC. All rights reserved. This +# © 2021-2025. Triad National Security, LLC. All rights reserved. This # program was produced under U.S. Government contract 89233218CNA000001 # for Los Alamos National Laboratory (LANL), which is operated by Triad # National Security, LLC for the U.S. Department of Energy/National @@ -17,3 +17,4 @@ matid=4272 name=steel +ionization = true diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/generate_files.cpp index 3966a9dc90..eeafb75293 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/generate_files.cpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -46,7 +46,8 @@ using namespace EospacWrapper; herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRhoBounds, const Bounds &lTBounds, const Bounds &leBounds, - const std::string &name, Verbosity eospacWarn) { + const std::string &name, const bool addSubtables, + Verbosity eospacWarn) { const int matid = metadata.matid; std::string sMatid = std::to_string(matid); @@ -94,49 +95,47 @@ herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRh coldGroup = H5Gcreate(matGroup, SP5::Depends::coldCurve, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - { - DataBox P, sie, T, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho, mask; - eosDataOfRhoSie(matid, lRhoBounds, leBounds, P, T, bMod, dPdRho, dPdE, dTdRho, dTdE, - dEdRho, mask, eospacWarn); - status += P.saveHDF(leGroup, SP5::Fields::P); - status += T.saveHDF(leGroup, SP5::Fields::T); - status += bMod.saveHDF(leGroup, SP5::Fields::bMod); - status += dPdRho.saveHDF(leGroup, SP5::Fields::dPdRho); - status += dPdE.saveHDF(leGroup, SP5::Fields::dPdE); - status += dTdRho.saveHDF(leGroup, SP5::Fields::dTdRho); - status += dTdE.saveHDF(leGroup, SP5::Fields::dTdE); - status += dEdRho.saveHDF(leGroup, SP5::Fields::dEdRho); - status += mask.saveHDF(leGroup, SP5::Fields::mask); - } - - { - DataBox P, sie, T, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho, dEdT, mask; - eosDataOfRhoT(matid, lRhoBounds, lTBounds, P, sie, bMod, dPdRho, dPdE, dTdRho, dTdE, - dEdRho, dEdT, mask, eospacWarn); - status += P.saveHDF(lTGroup, SP5::Fields::P); - status += sie.saveHDF(lTGroup, SP5::Fields::sie); - status += bMod.saveHDF(lTGroup, SP5::Fields::bMod); - status += dPdRho.saveHDF(lTGroup, SP5::Fields::dPdRho); - status += dPdE.saveHDF(lTGroup, SP5::Fields::dPdE); - status += dTdRho.saveHDF(lTGroup, SP5::Fields::dTdRho); - status += dTdE.saveHDF(lTGroup, SP5::Fields::dTdE); - status += dEdRho.saveHDF(lTGroup, SP5::Fields::dEdRho); - status += dEdT.saveHDF(lTGroup, SP5::Fields::dEdT); - status += mask.saveHDF(lTGroup, SP5::Fields::mask); - } + status += saveTablesRhoSie(leGroup, matid, TableSplit::Total, lRhoBounds, leBounds, + eospacWarn); + status += + saveTablesRhoT(lTGroup, matid, TableSplit::Total, lRhoBounds, lTBounds, eospacWarn); { DataBox P, sie, dPdRho, dEdRho, bMod, mask, transitionMask; eosColdCurves(matid, lRhoBounds, P, sie, dPdRho, dEdRho, bMod, mask, eospacWarn); - eosColdCurveMask(matid, lRhoBounds, leBounds.grid.nPoints(), sie, transitionMask, - eospacWarn); + // currently unused + // eosColdCurveMask(matid, lRhoBounds, leBounds.grid.nPoints(), sie, transitionMask, + // eospacWarn); status += P.saveHDF(coldGroup, SP5::Fields::P); status += sie.saveHDF(coldGroup, SP5::Fields::sie); status += bMod.saveHDF(coldGroup, SP5::Fields::bMod); status += dPdRho.saveHDF(coldGroup, SP5::Fields::dPdRho); status += dEdRho.saveHDF(coldGroup, SP5::Fields::dEdRho); - status += mask.saveHDF(coldGroup, SP5::Fields::mask); - status += transitionMask.saveHDF(coldGroup, SP5::Fields::transitionMask); + // currently unused + // status += mask.saveHDF(coldGroup, SP5::Fields::mask); + // status += transitionMask.saveHDF(coldGroup, SP5::Fields::transitionMask); + } + + if (addSubtables) { + int i = 0; + std::vector splits = {TableSplit::ElectronOnly, TableSplit::IonCold}; + std::vector grpnames = {SP5::SubTable::electronOnly, + SP5::SubTable::ionCold}; + for (auto split : splits) { + std::string grpname = grpnames[i++]; + { + hid_t grp = + H5Gcreate(leGroup, grpname.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + status += saveTablesRhoSie(grp, matid, split, lRhoBounds, leBounds, eospacWarn); + status += H5Gclose(grp); + } + { + hid_t grp = + H5Gcreate(lTGroup, grpname.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + status += saveTablesRhoT(grp, matid, split, lRhoBounds, lTBounds, eospacWarn); + status += H5Gclose(grp); + } + } } status += H5Gclose(leGroup); @@ -219,8 +218,14 @@ herr_t saveAllMaterials(const std::string &savename, << lRhoBounds << lTBounds << leBounds << std::endl; } - status += - saveMaterial(file, metadata, lRhoBounds, lTBounds, leBounds, name, eospacWarn); + const bool add_subtables = params[i].Get("ionization", false); + if (eospacWarn == Verbosity::Debug) { + std::cout << "Adding subtables for partial ionization? " << add_subtables + << std::endl; + } + + status += saveMaterial(file, metadata, lRhoBounds, lTBounds, leBounds, name, + add_subtables, eospacWarn); if (status != H5_SUCCESS) { std::cerr << "WARNING: problem with HDf5" << std::endl; } @@ -234,6 +239,45 @@ herr_t saveAllMaterials(const std::string &savename, return status; } +herr_t saveTablesRhoSie(hid_t loc, int matid, TableSplit split, const Bounds &lRhoBounds, + const Bounds &leBounds, Verbosity eospacWarn) { + herr_t status = 0; + DataBox P, T, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho, mask; + eosDataOfRhoSie(matid, split, lRhoBounds, leBounds, P, T, bMod, dPdRho, dPdE, dTdRho, + dTdE, dEdRho, mask, eospacWarn); + status += P.saveHDF(loc, SP5::Fields::P); + status += T.saveHDF(loc, SP5::Fields::T); + status += bMod.saveHDF(loc, SP5::Fields::bMod); + status += dPdRho.saveHDF(loc, SP5::Fields::dPdRho); + status += dPdE.saveHDF(loc, SP5::Fields::dPdE); + status += dTdRho.saveHDF(loc, SP5::Fields::dTdRho); + status += dTdE.saveHDF(loc, SP5::Fields::dTdE); + status += dEdRho.saveHDF(loc, SP5::Fields::dEdRho); + // currently unused + // status += mask.saveHDF(loc, SP5::Fields::mask); + return status; +} + +herr_t saveTablesRhoT(hid_t loc, int matid, TableSplit split, const Bounds &lRhoBounds, + const Bounds &lTBounds, Verbosity eospacWarn) { + herr_t status = 0; + DataBox P, sie, bMod, dPdRho, dPdE, dTdRho, dTdE, dEdRho, dEdT, mask; + eosDataOfRhoT(matid, split, lRhoBounds, lTBounds, P, sie, bMod, dPdRho, dPdE, dTdRho, + dTdE, dEdRho, dEdT, mask, eospacWarn); + status += P.saveHDF(loc, SP5::Fields::P); + status += sie.saveHDF(loc, SP5::Fields::sie); + status += bMod.saveHDF(loc, SP5::Fields::bMod); + status += dPdRho.saveHDF(loc, SP5::Fields::dPdRho); + status += dPdE.saveHDF(loc, SP5::Fields::dPdE); + status += dTdRho.saveHDF(loc, SP5::Fields::dTdRho); + status += dTdE.saveHDF(loc, SP5::Fields::dTdE); + status += dEdRho.saveHDF(loc, SP5::Fields::dEdRho); + status += dEdT.saveHDF(loc, SP5::Fields::dEdT); + // Currently unused + // status += mask.saveHDF(loc, SP5::Fields::mask); + return status; +} + void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params ¶ms, Bounds &lRhoBounds, Bounds &lTBounds, Bounds &leBounds) { diff --git a/sesame2spiner/generate_files.hpp b/sesame2spiner/generate_files.hpp index 914ec62ed8..543a60fb7f 100644 --- a/sesame2spiner/generate_files.hpp +++ b/sesame2spiner/generate_files.hpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -42,12 +42,25 @@ constexpr Real T_SPLIT_POINT_DEFAULT = 1e4; herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, const Bounds &lRhoBounds, const Bounds &lTBounds, const Bounds &leBounds, - const std::string &name, Verbosity eospacWarn = Verbosity::Quiet); + const std::string &name, const bool addSubtables, + Verbosity eospacWarn = Verbosity::Quiet); +inline herr_t saveMaterial(hid_t loc, const SesameMetadata &metadata, + const Bounds &lRhoBounds, const Bounds &lTBounds, + const Bounds &leBounds, const std::string &name, + Verbosity eospacWarn = Verbosity::Quiet) { + return saveMaterial(loc, metadata, lRhoBounds, lTBounds, leBounds, name, false, + eospacWarn); +} herr_t saveAllMaterials(const std::string &savename, const std::vector &filenames, bool printMetadata, Verbosity eospacWarn); +herr_t saveTablesRhoSie(hid_t loc, int matid, TableSplit split, const Bounds &lRhoBounds, + const Bounds &leBounds, Verbosity eospacWarn = Verbosity::Quiet); +herr_t saveTablesRhoT(hid_t loc, int matid, TableSplit split, const Bounds &lRhoBounds, + const Bounds &lTBounds, Verbosity eospacWarn = Verbosity::Quiet); + void getMatBounds(int i, int matid, const SesameMetadata &metadata, const Params ¶ms, Bounds &lRhoBounds, Bounds &lTBounds, Bounds &leBounds); diff --git a/sesame2spiner/io_eospac.cpp b/sesame2spiner/io_eospac.cpp index 0953029223..2eac41d536 100644 --- a/sesame2spiner/io_eospac.cpp +++ b/sesame2spiner/io_eospac.cpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -24,29 +24,32 @@ #include "io_eospac.hpp" // TODO: more error checking of bounds? -void eosDataOfRhoSie(int matid, const Bounds &lRhoBounds, const Bounds &leBounds, - DataBox &Ps, DataBox &Ts, DataBox &bMods, DataBox &dPdRho, - DataBox &dPde, DataBox &dTdRho, DataBox &dTde, DataBox &dEdRho, - DataBox &mask, Verbosity eospacWarn) { +void eosDataOfRhoSie(int matid, const TableSplit split, const Bounds &lRhoBounds, + const Bounds &leBounds, DataBox &Ps, DataBox &Ts, DataBox &bMods, + DataBox &dPdRho, DataBox &dPde, DataBox &dTdRho, DataBox &dTde, + DataBox &dEdRho, DataBox &mask, Verbosity eospacWarn) { using namespace EospacWrapper; constexpr int NT = 3; - constexpr EOS_INTEGER nXYPairs = 1; EOS_INTEGER tableHandle[NT]; EOS_INTEGER eospacPofRT, eospacTofRE, eospacEofRT; - EOS_INTEGER tableType[NT] = {EOS_Pt_DT, EOS_T_DUt, EOS_Ut_DT}; - - // Interpolatable vars - EOS_REAL var[1], dx[1], dy[1]; + EOS_INTEGER tableType[NT] = {impl::select(split, EOS_Pt_DT, EOS_Pe_DT, EOS_Pic_DT), + impl::select(split, EOS_T_DUt, EOS_T_DUe, EOS_T_DUic), + impl::select(split, EOS_Ut_DT, EOS_Ue_DT, EOS_Uic_DT)}; // indep vars std::vector rhos, sies; makeInterpPoints(rhos, lRhoBounds); makeInterpPoints(sies, leBounds); + EospacWrapper::eospacSplit apply_splitting = + (split == TableSplit::Total) ? eospacSplit::none : eospacSplit::splitNumProp; + // Load tables - eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Pt_DT", "EOS_T_DUt", "EOS_Ut_DT"}, - eospacWarn); + std::vector names = {"EOS_Pt_DT", "EOS_T_DUt", "EOS_Ut_DT"}; + impl::modifyNames(split, names); + eosSafeLoad(NT, matid, tableType, tableHandle, names, eospacWarn, false, 0.0, + eospacMonotonicity::none, false, apply_splitting, false); eospacPofRT = tableHandle[0]; eospacTofRE = tableHandle[1]; eospacEofRT = tableHandle[2]; @@ -65,62 +68,68 @@ void eosDataOfRhoSie(int matid, const Bounds &lRhoBounds, const Bounds &leBounds dEdRho.copyMetadata(Ps); mask.copyMetadata(Ps); + // Interpolatable vars + EOS_INTEGER nXYPairs = rhos.size() * sies.size(); + std::vector T_pack(nXYPairs), P_pack(nXYPairs), sie_pack(nXYPairs), + DTDR_E(nXYPairs), DTDE_R(nXYPairs), DPDR_T(nXYPairs), DPDT_R(nXYPairs), + DEDR_T(nXYPairs), DEDT_R(nXYPairs), rho_flat(nXYPairs), sie_flat(nXYPairs); + std::size_t iflat = 0; + for (std::size_t j = 0; j < rhos.size(); ++j) { + for (std::size_t i = 0; i < sies.size(); ++i) { + rho_flat[iflat] = densityToSesame(rhos[j]); + sie_flat[iflat] = sieToSesame(sies[i]); + iflat++; + } + } + const bool no_errors_tofre = eosSafeInterpolate( + &eospacTofRE, nXYPairs, rho_flat.data(), sie_flat.data(), T_pack.data(), + DTDR_E.data(), DTDE_R.data(), "TofRE", eospacWarn); + const bool no_errors_pofrt = eosSafeInterpolate( + &eospacPofRT, nXYPairs, rho_flat.data(), T_pack.data(), P_pack.data(), + DPDR_T.data(), DPDT_R.data(), "PofRT", eospacWarn); + const bool no_errors_eofrt = eosSafeInterpolate( + &eospacEofRT, nXYPairs, rho_flat.data(), T_pack.data(), sie_pack.data(), + DEDR_T.data(), DEDT_R.data(), "EofRT", eospacWarn); + const bool no_errors = no_errors_tofre && no_errors_pofrt && no_errors_eofrt; + // Loop by hand to ensure ordering ordering of independent // variables is under our control. + iflat = 0; for (size_t j = 0; j < rhos.size(); j++) { Real rho = densityToSesame(rhos[j]); for (size_t i = 0; i < sies.size(); i++) { - Real sie = sieToSesame(sies[i]); - // T - bool no_errors = true; - no_errors = no_errors && eosSafeInterpolate(&eospacTofRE, nXYPairs, &rho, &sie, var, - dx, dy, "TofRE", eospacWarn); - Real T = var[0]; - Real DTDR_E = dx[0]; - Real DTDE_R = dy[0]; - // P - no_errors = no_errors && eosSafeInterpolate(&eospacPofRT, nXYPairs, &rho, &T, var, - dx, dy, "PofRT", eospacWarn); - Real P = var[0]; - Real DPDR_T = dx[0]; - Real DPDT_R = dy[0]; - // Bulk modulus - no_errors = no_errors && eosSafeInterpolate(&eospacEofRT, nXYPairs, &rho, &T, var, - dx, dy, "EofRT", eospacWarn); - Real DEDR_T = dx[0]; - Real DEDT_R = dy[0]; - Real DPDE_R = DPDT_R / DEDT_R; - Real bMod = getBulkModulus(rho, P, DPDR_T, DPDE_R, DEDR_T); + Real DPDE_R = DPDT_R[iflat] / DEDT_R[iflat]; + Real bMod = + getBulkModulus(rho, P_pack[iflat], DPDR_T[iflat], DPDE_R, DEDR_T[iflat]); // Fill DataBoxes - Ts(j, i) = temperatureFromSesame(T); - Ps(j, i) = pressureFromSesame(P); + Ts(j, i) = temperatureFromSesame(T_pack[iflat]); + Ps(j, i) = pressureFromSesame(P_pack[iflat]); bMods(j, i) = bulkModulusFromSesame(std::max(bMod, 0.0)); - dPdRho(j, i) = pressureFromSesame(DPDR_T + DTDR_E * DPDT_R); - dPde(j, i) = sieToSesame(pressureFromSesame(DPDT_R * DTDE_R)); - dTdRho(j, i) = temperatureFromSesame(DTDR_E); - dTde(j, i) = sieToSesame(temperatureFromSesame(DTDE_R)); - dEdRho(j, i) = densityToSesame(sieFromSesame(DEDR_T)); + dPdRho(j, i) = pressureFromSesame(DPDR_T[iflat] + DTDR_E[iflat] * DPDT_R[iflat]); + dPde(j, i) = sieToSesame(pressureFromSesame(DPDT_R[iflat] * DTDE_R[iflat])); + dTdRho(j, i) = temperatureFromSesame(DTDR_E[iflat]); + dTde(j, i) = sieToSesame(temperatureFromSesame(DTDE_R[iflat])); + dEdRho(j, i) = densityToSesame(sieFromSesame(DEDR_T[iflat])); mask(j, i) = no_errors ? 1.0 : 0.0; + iflat++; } } eosSafeDestroy(NT, tableHandle, eospacWarn); } -void eosDataOfRhoT(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds, - DataBox &Ps, DataBox &sies, DataBox &bMods, DataBox &dPdRho, - DataBox &dPdE, DataBox &dTdRho, DataBox &dTde, DataBox &dEdRho, - DataBox &dEdT, DataBox &mask, Verbosity eospacWarn) { +void eosDataOfRhoT(int matid, const TableSplit split, const Bounds &lRhoBounds, + const Bounds &lTBounds, DataBox &Ps, DataBox &sies, DataBox &bMods, + DataBox &dPdRho, DataBox &dPdE, DataBox &dTdRho, DataBox &dTde, + DataBox &dEdRho, DataBox &dEdT, DataBox &mask, Verbosity eospacWarn) { using namespace EospacWrapper; constexpr int NT = 3; - constexpr EOS_INTEGER nXYPairs = 1; EOS_INTEGER tableHandle[NT]; EOS_INTEGER eospacPofRT, eospacTofRE, eospacEofRT; - EOS_INTEGER tableType[NT] = {EOS_Pt_DT, EOS_T_DUt, EOS_Ut_DT}; - - // Interpolatable vars - EOS_REAL var[1], dx[1], dy[1]; + EOS_INTEGER tableType[NT] = {impl::select(split, EOS_Pt_DT, EOS_Pe_DT, EOS_Pic_DT), + impl::select(split, EOS_T_DUt, EOS_T_DUe, EOS_T_DUic), + impl::select(split, EOS_Ut_DT, EOS_Ue_DT, EOS_Uic_DT)}; // indep vars std::vector rhos, Ts; @@ -128,8 +137,13 @@ void eosDataOfRhoT(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds, makeInterpPoints(Ts, lTBounds); // Load tables - eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Pt_DT", "EOS_T_DUt", "EOS_Ut_DT"}, - eospacWarn); + EospacWrapper::eospacSplit apply_splitting = + (split == TableSplit::Total) ? eospacSplit::none : eospacSplit::splitNumProp; + + std::vector names = {"EOS_Pt_DT", "EOS_T_DUt", "EOS_Ut_DT"}; + impl::modifyNames(split, names); + eosSafeLoad(NT, matid, tableType, tableHandle, names, eospacWarn, false, 0.0, + eospacMonotonicity::none, false, apply_splitting, false); eospacPofRT = tableHandle[0]; eospacTofRE = tableHandle[1]; eospacEofRT = tableHandle[2]; @@ -148,43 +162,55 @@ void eosDataOfRhoT(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds, dEdT.copyMetadata(Ps); mask.copyMetadata(Ps); - // Loop by hand to ensure ordering ordering of independent - // variables is under our control. + // Interpolatable vars + EOS_INTEGER nXYPairs = rhos.size() * Ts.size(); + std::vector P_pack(nXYPairs), DPDR_T(nXYPairs), DPDT_R(nXYPairs), + E_pack(nXYPairs), DEDR_T(nXYPairs), DEDT_R(nXYPairs), T_pack(nXYPairs), + DTDR_E(nXYPairs), DTDE_R(nXYPairs), rho_flat(nXYPairs), T_flat(nXYPairs); + + // prepare flat data structures + std::size_t iflat = 0; + for (std::size_t j = 0; j < rhos.size(); ++j) { + for (std::size_t i = 0; i < Ts.size(); ++i) { + rho_flat[iflat] = densityToSesame(rhos[j]); + T_flat[iflat] = temperatureToSesame(Ts[i]); + iflat++; + } + } + + // Pressure + const bool no_errors_prt = eosSafeInterpolate( + &eospacPofRT, nXYPairs, rho_flat.data(), T_flat.data(), P_pack.data(), + DPDR_T.data(), DPDT_R.data(), "PofRT", eospacWarn); + // Energy + const bool no_errors_ert = eosSafeInterpolate( + &eospacEofRT, nXYPairs, rho_flat.data(), T_flat.data(), E_pack.data(), + DEDR_T.data(), DEDT_R.data(), "EofRT", eospacWarn); + // T derivatives + const bool no_errors_tre = eosSafeInterpolate( + &eospacTofRE, nXYPairs, rho_flat.data(), E_pack.data(), T_pack.data(), + DTDR_E.data(), DTDE_R.data(), "TofRE", eospacWarn); + const bool no_errors = no_errors_prt && no_errors_ert && no_errors_tre; + + // fill databoxes + iflat = 0; for (size_t j = 0; j < rhos.size(); j++) { Real rho = densityToSesame(rhos[j]); for (size_t i = 0; i < Ts.size(); i++) { - Real T = temperatureToSesame(Ts[i]); - bool no_errors = true; - // Pressure - no_errors = no_errors && eosSafeInterpolate(&eospacPofRT, nXYPairs, &rho, &T, var, - dx, dy, "PofRT", eospacWarn); - Real P = var[0]; - Real DPDR_T = dx[0]; - Real DPDT_R = dy[0]; - // Energy - no_errors = no_errors && eosSafeInterpolate(&eospacEofRT, nXYPairs, &rho, &T, var, - dx, dy, "EofRT", eospacWarn); - Real E = var[0]; - Real DEDR_T = dx[0]; - Real DEDT_R = dy[0]; - // T derivatives - no_errors = no_errors && eosSafeInterpolate(&eospacTofRE, nXYPairs, &rho, &E, var, - dx, dy, "TofRE", eospacWarn); - Real DTDR_E = dx[0]; - Real DTDE_R = dy[0]; - Real DPDE_R = DPDT_R / DEDT_R; - Real bMod = getBulkModulus(rho, P, DPDR_T, DPDE_R, DEDR_T); - // Fill DataBoxes - Ps(j, i) = pressureFromSesame(P); - sies(j, i) = sieFromSesame(E); + Real DPDE_R = DPDT_R[iflat] / DEDT_R[iflat]; + Real bMod = + getBulkModulus(rho, P_pack[iflat], DPDR_T[iflat], DPDE_R, DEDR_T[iflat]); + Ps(j, i) = pressureFromSesame(P_pack[iflat]); + sies(j, i) = sieFromSesame(E_pack[iflat]); bMods(j, i) = bulkModulusFromSesame(std::max(bMod, 0.0)); - dPdRho(j, i) = pressureFromSesame(DPDR_T + DTDR_E * DPDT_R); - dPdE(j, i) = sieToSesame(pressureFromSesame(DPDT_R * DTDE_R)); - dTdRho(j, i) = temperatureFromSesame(DTDR_E); - dTde(j, i) = sieToSesame(temperatureFromSesame(DTDE_R)); - dEdRho(j, i) = densityToSesame(sieFromSesame(DEDR_T)); - dEdT(j, i) = sieFromSesame(temperatureToSesame(DEDT_R)); + dPdRho(j, i) = pressureFromSesame(DPDR_T[iflat] + DTDR_E[iflat] * DPDT_R[iflat]); + dPdE(j, i) = sieToSesame(pressureFromSesame(DPDT_R[iflat] * DTDE_R[iflat])); + dTdRho(j, i) = temperatureFromSesame(DTDR_E[iflat]); + dTde(j, i) = sieToSesame(temperatureFromSesame(DTDE_R[iflat])); + dEdRho(j, i) = densityToSesame(sieFromSesame(DEDR_T[iflat])); + dEdT(j, i) = sieFromSesame(temperatureToSesame(DEDT_R[iflat])); mask(j, i) = no_errors ? 1.0 : 0.0; + iflat++; } } eosSafeDestroy(NT, tableHandle, eospacWarn); @@ -196,17 +222,14 @@ void eosColdCurves(int matid, const Bounds &lRhoBounds, DataBox &Ps, DataBox &si using namespace EospacWrapper; constexpr int NT = 2; - constexpr EOS_INTEGER nXYPairs = 1; EOS_INTEGER tableHandle[NT]; EOS_INTEGER eospacPColdCurve, eospacSieColdCurve; EOS_INTEGER tableType[NT] = {EOS_Pc_D, EOS_Uc_D}; - // Interpolatable vars - EOS_REAL var[1], dx[1], dy[1]; - // indep vars - std::vector rhos; + std::vector rhos, Ts; makeInterpPoints(rhos, lRhoBounds); + Ts.resize(rhos.size()); // Load tables eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Pc_D", "EOS_Uc_D"}, eospacWarn); @@ -222,31 +245,40 @@ void eosColdCurves(int matid, const Bounds &lRhoBounds, DataBox &Ps, DataBox &si bMods.copyMetadata(Ps); mask.copyMetadata(Ps); - // loop by hand to ensure ordering - // TODO(JMM): this is probably not necessary - for (size_t i = 0; i < rhos.size(); i++) { - Real rho = densityToSesame(rhos[i]); - Real T = temperatureToSesame(0); - bool no_errors = true; - // pressure cold curve - no_errors = no_errors && eosSafeInterpolate(&eospacPColdCurve, nXYPairs, &rho, &T, - var, dx, dy, "PCold", eospacWarn); - Real P = var[0]; - Real DPDR_T = dx[0]; - // energy cold curve - no_errors = no_errors && eosSafeInterpolate(&eospacSieColdCurve, nXYPairs, &rho, &T, - var, dx, dy, "sieCold", eospacWarn); - Real sie = var[0]; - Real DEDR_T = dx[0]; - // bulk modulus cold curev - Real bMod = getBulkModulus(rho, P, DPDR_T, 0, 0); + for (std::size_t i = 0; i < rhos.size(); i++) { + rhos[i] = densityToSesame(rhos[i]); + Ts[i] = temperatureToSesame(0); + } + EOS_INTEGER nXYPairs = rhos.size(); + + // Interpolatable vars + std::vector P_pack, DPDR_T, sie_pack, DEDR_T, dy; + P_pack.resize(rhos.size()); + DPDR_T.resize(rhos.size()); + sie_pack.resize(rhos.size()); + DEDR_T.resize(rhos.size()); + dy.resize(rhos.size()); + + // Vector EOSPAC calls + bool no_errors = true; + no_errors = no_errors && eosSafeInterpolate(&eospacPColdCurve, nXYPairs, rhos.data(), + Ts.data(), P_pack.data(), DPDR_T.data(), + dy.data(), "PCold", eospacWarn); + no_errors = no_errors && eosSafeInterpolate(&eospacSieColdCurve, nXYPairs, rhos.data(), + Ts.data(), sie_pack.data(), DEDR_T.data(), + dy.data(), "sieCold", eospacWarn); + + // fill in vals + for (std::size_t i = 0; i < rhos.size(); i++) { + // bulk modulus cold curve + Real bMod = getBulkModulus(rhos[i], P_pack[i], DPDR_T[i], 0, 0); // fill DataBoxes - Ps(i) = pressureFromSesame(P); - sies(i) = sieFromSesame(sie); + Ps(i) = pressureFromSesame(P_pack[i]); + sies(i) = sieFromSesame(sie_pack[i]); bMods(i) = bulkModulusFromSesame(std::max(bMod, 0.0)); - dPdRho(i) = pressureFromSesame(DPDR_T); // on cold curve DTDR_E = 0 - dEdRho(i) = sieFromSesame(DEDR_T); - mask(i) = no_errors ? 1.0 : 0.0; + dPdRho(i) = pressureFromSesame(DPDR_T[i]); // on cold curve DTDR_E = 0 + dEdRho(i) = sieFromSesame(DEDR_T[i]); + mask(i) = no_errors ? 1.0 : 0.0; // TODO(JMM): Currently unused. } } @@ -279,7 +311,7 @@ void eosColdCurveMask(int matid, const Bounds &lRhoBounds, const int numSie, Bounds coldBounds(-1, 1, numSie, false); mask.setRange(0, coldBounds.grid); - // loop and fill dummy variable just to see if EOSPAC errors out + // loop and fill dummy variable just to se.e if EOSPAC errors out for (size_t j = 0; j < rhos.size(); j++) { Real rho = densityToSesame(rhos[j]); Real sieCold = sieColdCurve(j); @@ -302,3 +334,16 @@ void makeInterpPoints(std::vector &v, const Bounds &b) { v[i] = b.i2lin(i); } } +namespace impl { +void modifyNames(TableSplit split, std::vector &names) { + if (split != TableSplit::Total) { + auto rt = std::regex("t"); + std::string newstr = (split == TableSplit::ElectronOnly) ? "e" : "ic"; + for (auto &s : names) { + if (s != "EOS_D_PtT") { + s = std::regex_replace(s, rt, newstr); + } + } + } +} +} // namespace impl diff --git a/sesame2spiner/io_eospac.hpp b/sesame2spiner/io_eospac.hpp index 454d44ec0f..cd73390cb5 100644 --- a/sesame2spiner/io_eospac.hpp +++ b/sesame2spiner/io_eospac.hpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -17,6 +17,7 @@ #ifndef _SESAME2SPINER_IO_EOSPAC_HPP_ #define _SESAME2SPINER_IO_EOSPAC_HPP_ +#include #include #include // eospac API @@ -26,6 +27,7 @@ #endif #include +#include #include #include #include @@ -34,19 +36,38 @@ using EospacWrapper::Verbosity; constexpr int NGRIDS = 3; +using singularity::TableSplit; using Bounds = singularity::table_utils::Bounds; using Grid_t = Spiner::PiecewiseGrid1D; using DataBox = Spiner::DataBox; -void eosDataOfRhoSie(int matid, const Bounds &lRhoBounds, const Bounds &leBounds, - DataBox &P, DataBox &T, DataBox &bMods, DataBox &dPdRho, - DataBox &dPdE, DataBox &dTdRho, DataBox &dTdE, DataBox &dEdRho, - DataBox &mask, Verbosity eospacWarn = Verbosity::Quiet); +void eosDataOfRhoSie(int matid, const TableSplit split, const Bounds &lRhoBounds, + const Bounds &leBounds, DataBox &P, DataBox &T, DataBox &bMods, + DataBox &dPdRho, DataBox &dPdE, DataBox &dTdRho, DataBox &dTdE, + DataBox &dEdRho, DataBox &mask, + Verbosity eospacWarn = Verbosity::Quiet); +inline void eosDataOfRhoSie(int matid, const Bounds &lRhoBounds, const Bounds &leBounds, + DataBox &P, DataBox &T, DataBox &bMods, DataBox &dPdRho, + DataBox &dPdE, DataBox &dTdRho, DataBox &dTdE, + DataBox &dEdRho, DataBox &mask, + Verbosity eospacWarn = Verbosity::Quiet) { + eosDataOfRhoSie(matid, TableSplit::Total, lRhoBounds, leBounds, P, T, bMods, dPdRho, + dPdE, dTdRho, dTdE, dEdRho, mask, eospacWarn); +} -void eosDataOfRhoT(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds, - DataBox &Ps, DataBox &sies, DataBox &bMods, DataBox &dPdRho, - DataBox &dPdE, DataBox &dTdRho, DataBox &dTdE, DataBox &dEdRho, - DataBox &dEdT, DataBox &mask, Verbosity eospacWarn = Verbosity::Quiet); +void eosDataOfRhoT(int matid, const TableSplit split, const Bounds &lRhoBounds, + const Bounds &lTBounds, DataBox &Ps, DataBox &sies, DataBox &bMods, + DataBox &dPdRho, DataBox &dPdE, DataBox &dTdRho, DataBox &dTdE, + DataBox &dEdRho, DataBox &dEdT, DataBox &mask, + Verbosity eospacWarn = Verbosity::Quiet); +inline void eosDataOfRhoT(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds, + DataBox &Ps, DataBox &sies, DataBox &bMods, DataBox &dPdRho, + DataBox &dPdE, DataBox &dTdRho, DataBox &dTdE, DataBox &dEdRho, + DataBox &dEdT, DataBox &mask, + Verbosity eospacWarn = Verbosity::Quiet) { + eosDataOfRhoT(matid, TableSplit::Total, lRhoBounds, lTBounds, Ps, sies, bMods, dPdRho, + dPdE, dTdRho, dTdE, dEdRho, dEdT, mask, eospacWarn); +} void eosColdCurves(int matid, const Bounds &lRhoBounds, DataBox &Ps, DataBox &sies, DataBox &dPdRho, DataBox &dEdRho, DataBox &bMod, DataBox &mask, @@ -58,4 +79,16 @@ void eosColdCurveMask(int matid, const Bounds &lRhoBounds, const int numSie, void makeInterpPoints(std::vector &v, const Bounds &b); +namespace impl { +template +T select(TableSplit split, T a, T b, T c) { + if (split == TableSplit::Total) return a; + if (split == TableSplit::ElectronOnly) + return b; + else // if (split == TableSplit::IonCold) + return c; +} +void modifyNames(TableSplit split, std::vector &names); +} // namespace impl + #endif // _SESAME2SPINER_IO_EOSPAC_HPP_ diff --git a/singularity-eos/base/constants.hpp b/singularity-eos/base/constants.hpp index 73934daf31..b1ab8cfc47 100644 --- a/singularity-eos/base/constants.hpp +++ b/singularity-eos/base/constants.hpp @@ -34,6 +34,7 @@ constexpr unsigned long all_values = (1 << 7) - 1; constexpr size_t MAX_NUM_LAMBDAS = 3; enum class DataStatus { Deallocated = 0, OnDevice = 1, OnHost = 2, UnManaged = 3 }; enum class TableStatus { OnTable = 0, OffBottom = 1, OffTop = 2 }; +enum class TableSplit { Total = 0, ElectronOnly = 1, IonCold = 2 }; constexpr Real ROOM_TEMPERATURE = 293; // K constexpr Real ATMOSPHERIC_PRESSURE = 1e6; diff --git a/singularity-eos/base/sp5/singularity_eos_sp5.hpp b/singularity-eos/base/sp5/singularity_eos_sp5.hpp index a7f2832045..002da520ea 100644 --- a/singularity-eos/base/sp5/singularity_eos_sp5.hpp +++ b/singularity-eos/base/sp5/singularity_eos_sp5.hpp @@ -1,7 +1,7 @@ //------------------------------------------------------------------------------ // Extra definitions to the SP5 file format required for singularity-eos // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -28,6 +28,11 @@ constexpr char logRhoLogT[] = "dependsLogRhoLogT"; constexpr char coldCurve[] = "coldCurve"; } // namespace Depends +namespace SubTable { +constexpr char electronOnly[] = "electronOnly"; +constexpr char ionCold[] = "ionCold"; +} // namespace SubTable + namespace Offsets { constexpr char messageName[] = "interpretation"; constexpr char message[] = diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index 0e6ef06e0c..82c579d5bd 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -129,14 +130,39 @@ inline void SetUpOutputScalingOption(EOS_INTEGER options[], EOS_REAL values[], class EOSPAC : public EosBase { public: + struct Options { + TableSplit split = TableSplit::Total; + bool invert_at_setup = false; + Real insert_data = 0.0; + eospacMonotonicity monotonicity = eospacMonotonicity::none; + bool apply_smoothing = false; + eospacSplit apply_splitting = eospacSplit::none; + bool linear_interp = false; + }; + inline EOSPAC() = default; ////add options here... and elsewhere - inline EOSPAC(int matid, bool invert_at_setup = false, Real insert_data = 0.0, + // Complete constructor + inline EOSPAC(int matid, TableSplit split, bool invert_at_setup = false, + Real insert_data = 0.0, eospacMonotonicity monotonicity = eospacMonotonicity::none, bool apply_smoothing = false, eospacSplit apply_splitting = eospacSplit::none, bool linear_interp = false); + // Overloads with fewer positional arguments + inline EOSPAC(int matid, bool invert_at_setup = false, Real insert_data = 0.0, + eospacMonotonicity monotonicity = eospacMonotonicity::none, + bool apply_smoothing = false, + eospacSplit apply_splitting = eospacSplit::none, + bool linear_interp = false) + : EOSPAC(matid, TableSplit::Total, invert_at_setup, insert_data, monotonicity, + apply_smoothing, apply_splitting, linear_interp) {} + inline EOSPAC(int matid, const Options &opts) + : EOSPAC(matid, opts.split, opts.invert_at_setup, opts.insert_data, + opts.monotonicity, opts.apply_smoothing, opts.apply_splitting, + opts.linear_interp) {} + PORTABLE_INLINE_FUNCTION void CheckParams() const { // TODO(JMM): More validation checks? PORTABLE_ALWAYS_REQUIRE(rho_min_ >= 0, "Non-negative minimum density"); @@ -1110,6 +1136,7 @@ class EOSPAC : public EosBase { static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::temperature; int matid_; + TableSplit split_; static constexpr int NT = 7; EOS_INTEGER PofRT_table_; EOS_INTEGER TofRE_table_; @@ -1158,21 +1185,41 @@ class EOSPAC : public EosBase { // Implementation details below // ====================================================================== -inline EOSPAC::EOSPAC(const int matid, bool invert_at_setup, Real insert_data, - EospacWrapper::eospacMonotonicity monotonicity, +inline EOSPAC::EOSPAC(const int matid, TableSplit split, bool invert_at_setup, + Real insert_data, EospacWrapper::eospacMonotonicity monotonicity, bool apply_smoothing, EospacWrapper::eospacSplit apply_splitting, bool linear_interp) - : matid_(matid) { + : matid_(matid), split_(split) { using namespace EospacWrapper; - EOS_INTEGER tableType[NT] = {EOS_Pt_DT, EOS_T_DUt, EOS_Ut_DT, EOS_D_PtT, - EOS_T_DPt, EOS_Pt_DUt, EOS_Uc_D}; - eosSafeLoad( - NT, matid, tableType, tablehandle, - std::vector({"EOS_Pt_DT", "EOS_T_DUt", "EOS_Ut_DT", "EOS_D_PtT", - "EOS_T_DPt", "EOS_Pt_DUt", "EOS_Uc_D"}), - Verbosity::Debug, invert_at_setup = invert_at_setup, insert_data = insert_data, - monotonicity = monotonicity, apply_smoothing = apply_smoothing, - apply_splitting = apply_splitting, linear_interp = linear_interp); + auto TableSelect = [&](auto a, auto b, auto c) { + if (split == TableSplit::Total) return a; + if (split == TableSplit::ElectronOnly) + return b; + else // if (split == TableSplit::IonCold) + return c; + }; + EOS_INTEGER tableType[NT] = {TableSelect(EOS_Pt_DT, EOS_Pe_DT, EOS_Pic_DT), + TableSelect(EOS_T_DUt, EOS_T_DUe, EOS_T_DUic), + TableSelect(EOS_Ut_DT, EOS_Ue_DT, EOS_Uic_DT), + EOS_D_PtT, + TableSelect(EOS_T_DPt, EOS_T_DPe, EOS_T_DPic), + TableSelect(EOS_Pt_DUt, EOS_Pe_DUe, EOS_Pic_DUic), + EOS_Uc_D}; + std::vector tableNames = {"EOS_Pt_DT", "EOS_T_DUt", "EOS_Ut_DT", + "EOS_D_PtT", "EOS_T_DPt", "EOS_Pt_DUt", + "EOS_Uc_D"}; + if (split != TableSplit::Total) { + auto rt = std::regex("t"); + std::string newstr = (split == TableSplit::ElectronOnly) ? "e" : "ic"; + for (auto &s : tableNames) { + if (s != "EOS_D_PtT") { + s = std::regex_replace(s, rt, newstr); + } + } + } + eosSafeLoad(NT, matid, tableType, tablehandle, tableNames, Verbosity::Debug, + invert_at_setup, insert_data, monotonicity, apply_smoothing, + apply_splitting, linear_interp); PofRT_table_ = tablehandle[0]; TofRE_table_ = tablehandle[1]; EofRT_table_ = tablehandle[2]; @@ -1325,6 +1372,9 @@ EOSPAC::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &b } if (output & thermalqs::density) { if (input & thermalqs::pressure && input & thermalqs::temperature) { + PORTABLE_REQUIRE(split_ == TableSplit::Total, + "Density of pressure and temperature only supported for total " + "tables at this time"); EOS_INTEGER table = RofPT_table_; eosSafeInterpolate(&table, nxypairs, P, T, R, dx, dy, "RofPT", Verbosity::Quiet); rho = R[0]; @@ -1560,6 +1610,10 @@ PORTABLE_INLINE_FUNCTION void EOSPAC::DensityEnergyFromPressureTemperature( #if SINGULARITY_ON_DEVICE PORTABLE_ALWAYS_ABORT("EOSPAC calls not supported on device\n"); #else + PORTABLE_REQUIRE(split_ == TableSplit::Total, + "Density of pressure and temperature only supported for total " + "tables at this time"); + using namespace EospacWrapper; EOS_REAL P[1] = {pressureToSesame(press)}; EOS_REAL T[1] = {temperatureToSesame(temp)}; diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 6af2e644d1..b9d5d9fd05 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -82,13 +82,22 @@ class SpinerEOSDependsRhoT : public EosBase { }; SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_) SG_ADD_BASE_CLASS_USINGS(SpinerEOSDependsRhoT); + inline SpinerEOSDependsRhoT(const std::string &filename, int matid, TableSplit split, + bool reproducibility_mode = false); inline SpinerEOSDependsRhoT(const std::string &filename, int matid, - bool reproduciblity_mode = false); + bool reproducibility_mode = false) + : SpinerEOSDependsRhoT(filename, matid, TableSplit::Total, reproducibility_mode) {} inline SpinerEOSDependsRhoT(const std::string &filename, - const std::string &materialName, + const std::string &materialName, TableSplit split, bool reproducibility_mode = false); + inline SpinerEOSDependsRhoT(const std::string &filename, + const std::string &materialName, + bool reproducibility_mode = false) + : SpinerEOSDependsRhoT(filename, materialName, TableSplit::Total, + reproducibility_mode) {} PORTABLE_INLINE_FUNCTION - SpinerEOSDependsRhoT() : memoryStatus_(DataStatus::Deallocated) {} + SpinerEOSDependsRhoT() + : memoryStatus_(DataStatus::Deallocated), split_(TableSplit::Total) {} inline SpinerEOSDependsRhoT GetOnDevice(); @@ -315,6 +324,7 @@ class SpinerEOSDependsRhoT : public EosBase { Real lRhoOffset_, lTOffset_; // offsets must be non-negative MeanAtomicProperties AZbar_; int matid_; + TableSplit split_; bool reproducible_; // whereAmI_ and status_ used only for reporting. They are not thread-safe. mutable TableStatus whereAmI_ = TableStatus::OnTable; @@ -360,11 +370,20 @@ class SpinerEOSDependsRhoSie : public EosBase { SG_ADD_BASE_CLASS_USINGS(SpinerEOSDependsRhoSie); PORTABLE_INLINE_FUNCTION SpinerEOSDependsRhoSie() : memoryStatus_(DataStatus::Deallocated) {} + inline SpinerEOSDependsRhoSie(const std::string &filename, int matid, TableSplit split, + bool reproducibility_mode = false); inline SpinerEOSDependsRhoSie(const std::string &filename, int matid, + bool reproducibility_mode = false) + : SpinerEOSDependsRhoSie(filename, matid, TableSplit::Total, reproducibility_mode) { + } + inline SpinerEOSDependsRhoSie(const std::string &filename, + const std::string &materialName, TableSplit split, bool reproducibility_mode = false); inline SpinerEOSDependsRhoSie(const std::string &filename, const std::string &materialName, - bool reproducibility_mode = false); + bool reproducibility_mode = false) + : SpinerEOSDependsRhoSie(filename, materialName, TableSplit::Total, + reproducibility_mode) {} inline SpinerEOSDependsRhoSie GetOnDevice(); PORTABLE_INLINE_FUNCTION void CheckParams() const { @@ -555,6 +574,7 @@ class SpinerEOSDependsRhoSie : public EosBase { thermalqs::density | thermalqs::temperature; // static constexpr const char _eos_type[] = "SpinerEOSDependsRhoSie"; int matid_; + TableSplit split_; MeanAtomicProperties AZbar_; bool reproducible_; mutable RootFinding1D::Status status_; @@ -631,8 +651,9 @@ class interp { } // namespace callable_interp inline SpinerEOSDependsRhoT::SpinerEOSDependsRhoT(const std::string &filename, int matid, + TableSplit split, bool reproducibility_mode) - : matid_(matid), reproducible_(reproducibility_mode), + : matid_(matid), split_(split), reproducible_(reproducibility_mode), status_(RootFinding1D::Status::SUCCESS), memoryStatus_(DataStatus::OnHost) { std::string matid_str = std::to_string(matid); @@ -649,7 +670,14 @@ inline SpinerEOSDependsRhoT::SpinerEOSDependsRhoT(const std::string &filename, i "Log mode used at runtime must be identical to the one used to generate the file!"); matGroup = H5Gopen(file, matid_str.c_str(), H5P_DEFAULT); - lTGroup = H5Gopen(matGroup, SP5::Depends::logRhoLogT, H5P_DEFAULT); + + std::string lTGroupName = SP5::Depends::logRhoLogT; + if (split == TableSplit::ElectronOnly) { + lTGroupName += (std::string("/") + SP5::SubTable::electronOnly); + } else if (split == TableSplit::IonCold) { + lTGroupName += (std::string("/") + SP5::SubTable::ionCold); + } + lTGroup = H5Gopen(matGroup, lTGroupName.c_str(), H5P_DEFAULT); coldGroup = H5Gopen(matGroup, SP5::Depends::coldCurve, H5P_DEFAULT); status += loadDataboxes_(matid_str, file, lTGroup, coldGroup); @@ -668,17 +696,25 @@ inline SpinerEOSDependsRhoT::SpinerEOSDependsRhoT(const std::string &filename, i inline SpinerEOSDependsRhoT::SpinerEOSDependsRhoT(const std::string &filename, const std::string &materialName, + TableSplit split, bool reproducibility_mode) - : reproducible_(reproducibility_mode), status_(RootFinding1D::Status::SUCCESS), - memoryStatus_(DataStatus::OnHost) { + : split_(split), reproducible_(reproducibility_mode), + status_(RootFinding1D::Status::SUCCESS), memoryStatus_(DataStatus::OnHost) { std::string matid_str; - hid_t file, matGroup, lTGroup, coldGroup; + hid_t file, matGroup, lTGroup, subGroup, coldGroup; herr_t status = H5_SUCCESS; file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); matGroup = H5Gopen(file, materialName.c_str(), H5P_DEFAULT); - lTGroup = H5Gopen(matGroup, SP5::Depends::logRhoLogT, H5P_DEFAULT); + + std::string lTGroupName = SP5::Depends::logRhoLogT; + if (split == TableSplit::ElectronOnly) { + lTGroupName += (std::string("/") + SP5::SubTable::electronOnly); + } else if (split == TableSplit::IonCold) { + lTGroupName += (std::string("/") + SP5::SubTable::ionCold); + } + lTGroup = H5Gopen(matGroup, lTGroupName.c_str(), H5P_DEFAULT); coldGroup = H5Gopen(matGroup, SP5::Depends::coldCurve, H5P_DEFAULT); status += @@ -1499,9 +1535,9 @@ SpinerEOSDependsRhoT::getLocDependsRhoT_(const Real lRho, const Real lT) const { } inline SpinerEOSDependsRhoSie::SpinerEOSDependsRhoSie(const std::string &filename, - int matid, + int matid, TableSplit split, bool reproducibility_mode) - : matid_(matid), reproducible_(reproducibility_mode), + : matid_(matid), split_(split), reproducible_(reproducibility_mode), status_(RootFinding1D::Status::SUCCESS), memoryStatus_(DataStatus::OnHost) { std::string matid_str = std::to_string(matid); @@ -1510,8 +1546,19 @@ inline SpinerEOSDependsRhoSie::SpinerEOSDependsRhoSie(const std::string &filenam file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); matGroup = H5Gopen(file, matid_str.c_str(), H5P_DEFAULT); - lTGroup = H5Gopen(matGroup, SP5::Depends::logRhoLogT, H5P_DEFAULT); - lEGroup = H5Gopen(matGroup, SP5::Depends::logRhoLogSie, H5P_DEFAULT); + + std::string lTGroupName = SP5::Depends::logRhoLogT; + std::string lEGroupName = SP5::Depends::logRhoLogSie; + if (split == TableSplit::ElectronOnly) { + lTGroupName += (std::string("/") + SP5::SubTable::electronOnly); + lEGroupName += (std::string("/") + SP5::SubTable::electronOnly); + } else if (split == TableSplit::IonCold) { + lTGroupName += (std::string("/") + SP5::SubTable::ionCold); + lEGroupName += (std::string("/") + SP5::SubTable::ionCold); + } + + lTGroup = H5Gopen(matGroup, lTGroupName.c_str(), H5P_DEFAULT); + lEGroup = H5Gopen(matGroup, lEGroupName.c_str(), H5P_DEFAULT); status += loadDataboxes_(matid_str, file, lTGroup, lEGroup); @@ -1527,9 +1574,10 @@ inline SpinerEOSDependsRhoSie::SpinerEOSDependsRhoSie(const std::string &filenam inline SpinerEOSDependsRhoSie::SpinerEOSDependsRhoSie(const std::string &filename, const std::string &materialName, + TableSplit split, bool reproducibility_mode) - : reproducible_(reproducibility_mode), status_(RootFinding1D::Status::SUCCESS), - memoryStatus_(DataStatus::OnHost) { + : split_(split), reproducible_(reproducibility_mode), + status_(RootFinding1D::Status::SUCCESS), memoryStatus_(DataStatus::OnHost) { std::string matid_str; hid_t file, matGroup, lTGroup, lEGroup; @@ -1537,8 +1585,19 @@ inline SpinerEOSDependsRhoSie::SpinerEOSDependsRhoSie(const std::string &filenam file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); matGroup = H5Gopen(file, materialName.c_str(), H5P_DEFAULT); - lTGroup = H5Gopen(matGroup, SP5::Depends::logRhoLogT, H5P_DEFAULT); - lEGroup = H5Gopen(matGroup, SP5::Depends::logRhoLogSie, H5P_DEFAULT); + + std::string lTGroupName = SP5::Depends::logRhoLogT; + std::string lEGroupName = SP5::Depends::logRhoLogSie; + if (split == TableSplit::ElectronOnly) { + lTGroupName += (std::string("/") + SP5::SubTable::electronOnly); + lEGroupName += (std::string("/") + SP5::SubTable::electronOnly); + } else if (split == TableSplit::IonCold) { + lTGroupName += (std::string("/") + SP5::SubTable::ionCold); + lEGroupName += (std::string("/") + SP5::SubTable::ionCold); + } + + lTGroup = H5Gopen(matGroup, lTGroupName.c_str(), H5P_DEFAULT); + lEGroup = H5Gopen(matGroup, lEGroupName.c_str(), H5P_DEFAULT); status += H5LTget_attribute_int(file, materialName.c_str(), SP5::Material::matid, &matid_); diff --git a/test/profile_eospac.cpp b/test/profile_eospac.cpp index 52daeb01e9..4a051d87fb 100644 --- a/test/profile_eospac.cpp +++ b/test/profile_eospac.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -34,8 +34,9 @@ #include #include -#include #include +#include +#include using namespace singularity; @@ -57,6 +58,9 @@ constexpr Real T_MIN = 1.8e2; // Kelvin constexpr Real T_MAX = 1e3; constexpr Real E_MIN = 1e9; constexpr Real E_MAX = 1e13; +using EOS = Variant, ShiftedEOS, ScaledEOS, + ScaledEOS>, BilinearRampEOS, + BilinearRampEOS>>; inline auto now() { return std::chrono::high_resolution_clock::now(); } diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp index 6a31ed9451..72d9761d73 100644 --- a/test/test_eos_tabulated.cpp +++ b/test/test_eos_tabulated.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -381,6 +381,68 @@ SCENARIO("SpinerEOS and EOSPAC Serialization", } } +SCENARIO("Tabulated EOS with partial ionization", + "[SpinerEOS][DependsRhoT][DependsRhoSie][EOSPAC][SplitTables]") { + using singularity::TableSplit; + GIVEN("Spiner EOSes for total, electron, and ion EOS") { + SpinerEOSDependsRhoT rt_tot(eosName, steelID, TableSplit::Total); + SpinerEOSDependsRhoT rt_e(eosName, steelID, TableSplit::ElectronOnly); + SpinerEOSDependsRhoT rt_ic(eosName, steelID, TableSplit::IonCold); + std::vector rt = {rt_tot, rt_e, rt_ic}; + + SpinerEOSDependsRhoSie re_tot(eosName, steelID, TableSplit::Total); + SpinerEOSDependsRhoSie re_e(eosName, steelID, TableSplit::ElectronOnly); + SpinerEOSDependsRhoSie re_ic(eosName, steelID, TableSplit::IonCold); + std::vector re = {re_tot, re_e, re_ic}; + + EOSPAC pac_tot(steelID, TableSplit::Total); + EOSPAC pac_e(steelID, TableSplit::ElectronOnly); + EOSPAC pac_ic(steelID, TableSplit::IonCold); + std::vector pac = {pac_tot, pac_e, pac_ic}; + + WHEN("We evaluate pressure from density and temperature") { + constexpr Real rho = 1e1; + constexpr Real T = 1e4; + const std::size_t N = rt.size(); + std::vector P_rt(N), P_re(N), P_pac(N); + for (std::size_t i = 0; i < N; ++i) { + P_rt[i] = rt[i].PressureFromDensityTemperature(rho, T); + P_re[i] = re[i].PressureFromDensityTemperature(rho, T); + P_pac[i] = pac[i].PressureFromDensityTemperature(rho, T); + } + THEN("The subtables add to the total for DependsRhoT") { + REQUIRE(isClose(P_rt[0], P_rt[1] + P_rt[2])); + } + THEN("The subtables add to the total for DependsRhoSie") { + REQUIRE(isClose(P_re[0], P_re[1] + P_re[2])); + } + THEN("The subtables add to the total for EOSPAC") { + REQUIRE(isClose(P_pac[0], P_pac[1] + P_pac[2])); + } + THEN("DependsRhoT agrees with EOSPAC") { + for (std::size_t i = 0; i < N; ++i) { + REQUIRE(isClose(P_rt[i], P_pac[i])); + } + } + THEN("DependsRhoSie agrees with EOSPAC") { + for (std::size_t i = 0; i < N; ++i) { + REQUIRE(isClose(P_re[i], P_pac[i])); + } + } + } + + for (auto &eos : rt) { + eos.Finalize(); + } + for (auto &eos : re) { + eos.Finalize(); + } + for (auto &eos : pac) { + eos.Finalize(); + } + } +} + #endif // SINGULARITY_USE_EOSPAC #endif // SINGULARITY_TEST_SESAME #endif // SPINER_USE_HDF