From df010781afa63969f3492db5beae5b76da9cf053 Mon Sep 17 00:00:00 2001 From: "Alex R. Long" Date: Thu, 9 Jul 2026 11:20:36 -0600 Subject: [PATCH 1/5] Initial commit of ipcress2spiner capability --- CMakeLists.txt | 4 + ipcress2spiner/CMakeLists.txt | 52 ++++ ipcress2spiner/generate_files.cpp | 457 ++++++++++++++++++++++++++++++ ipcress2spiner/generate_files.hpp | 72 +++++ ipcress2spiner/io_eospac.cpp | 424 +++++++++++++++++++++++++++ ipcress2spiner/io_eospac.hpp | 97 +++++++ ipcress2spiner/main.cpp | 291 +++++++++++++++++++ ipcress2spiner/parse_cli.cpp | 80 ++++++ ipcress2spiner/parse_cli.hpp | 71 +++++ ipcress2spiner/parser.cpp | 161 +++++++++++ ipcress2spiner/parser.hpp | 52 ++++ 11 files changed, 1761 insertions(+) create mode 100644 ipcress2spiner/CMakeLists.txt create mode 100644 ipcress2spiner/generate_files.cpp create mode 100644 ipcress2spiner/generate_files.hpp create mode 100644 ipcress2spiner/io_eospac.cpp create mode 100644 ipcress2spiner/io_eospac.hpp create mode 100644 ipcress2spiner/main.cpp create mode 100644 ipcress2spiner/parse_cli.cpp create mode 100644 ipcress2spiner/parse_cli.hpp create mode 100644 ipcress2spiner/parser.cpp create mode 100644 ipcress2spiner/parser.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 55119d1..e5c083f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -109,6 +109,10 @@ if(SINGULARITY_BUILD_TESTS) add_subdirectory(test) endif() +if(SINGULARITY_BUILD_IPCRESS2SPINER) + add_subdirectory(ipcress2spiner) +endif() + set(CPACK_RESOURCE_FILE_LICENSE "${PROJECT_SOURCE_DIR}/LICENSE") include(CPack) diff --git a/ipcress2spiner/CMakeLists.txt b/ipcress2spiner/CMakeLists.txt new file mode 100644 index 0000000..0329aac --- /dev/null +++ b/ipcress2spiner/CMakeLists.txt @@ -0,0 +1,52 @@ +#------------------------------------------------------------------------------ +# © 2021-2023. 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. +#------------------------------------------------------------------------------ +project(ipcress2spiner LANGUAGES CXX Fortran) + +add_executable(ipcress2spiner + #io_eospac.cpp + #io_eospac.hpp + #generate_files.cpp + #generate_files.hpp + parse_cli.cpp + parse_cli.hpp + #parser.cpp + #parser.hpp + main.cpp +) + +# Turn on verbosity right before your package search +#set(CMAKE_FIND_DEBUG_MODE ON) +find_package(gandolf) +#set(CMAKE_FIND_DEBUG_MODE OFF) + +target_include_directories(ipcress2spiner + PUBLIC + $ + $ + ${GANDOLF_INCLUDE_DIRECTORY} +) + + +target_link_libraries(ipcress2spiner + PRIVATE + singularity-opac::singularity-opac + ${GANDOLF_LIBRARY} + PortsofCall::PortsofCall +) +set_property(TARGET ipcress2spiner PROPERTY LINKER_LANGUAGE Fortran) + +install(TARGETS ipcress2spiner DESTINATION ${CMAKE_INSTALL_BINDIR}) + + +# TODO: Add tests for sesame2spiner here. diff --git a/ipcress2spiner/generate_files.cpp b/ipcress2spiner/generate_files.cpp new file mode 100644 index 0000000..e5453f9 --- /dev/null +++ b/ipcress2spiner/generate_files.cpp @@ -0,0 +1,457 @@ +//====================================================================== +// sesame2spiner tool for converting eospac to spiner +// Author: Jonah Miller (jonahm@lanl.gov) +// © 2021-2025. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//====================================================================== + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#ifndef SPINER_USE_HDF +#error "HDF5 must be enabled" +#endif // SPINER_USE_HDF + +#include +#include +#include +#include +#include +#include +#include + +#include "generate_files.hpp" +#include "io_eospac.hpp" +#include "parse_cli.hpp" +#include "parser.hpp" + +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, const bool addSubtables, + Verbosity eospacWarn) { + + const int matid = metadata.matid; + std::string sMatid = std::to_string(matid); + + herr_t status = 0; + 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); + + // Dependent variables metadata + status += H5LTset_attribute_string(loc, sMatid.c_str(), SP5::Offsets::messageName, + SP5::Offsets::message); + status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Offsets::rho, + &lRhoBounds.offset, 1); + status += + H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Offsets::T, &lTBounds.offset, 1); + status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Offsets::sie, + &leBounds.offset, 1); + + // Material metadata + status += + H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Material::exchangeCoefficient, + &metadata.exchangeCoefficient, 1); + status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Material::meanAtomicMass, + &metadata.meanAtomicMass, 1); + status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Material::meanAtomicNumber, + &metadata.meanAtomicNumber, 1); + status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Material::solidBulkModulus, + &metadata.solidBulkModulus, 1); + status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Material::normalDensity, + &metadata.normalDensity, 1); + status += H5LTset_attribute_string(loc, sMatid.c_str(), SP5::Material::comments, + metadata.comments.c_str()); + status += H5LTset_attribute_int(loc, sMatid.c_str(), SP5::Material::matid, + &metadata.matid, 1); + status += H5LTset_attribute_string(loc, sMatid.c_str(), SP5::Material::name, + metadata.name.c_str()); + + lTGroup = H5Gcreate(matGroup, SP5::Depends::logRhoLogT, H5P_DEFAULT, H5P_DEFAULT, + H5P_DEFAULT); + + leGroup = H5Gcreate(matGroup, SP5::Depends::logRhoLogSie, H5P_DEFAULT, H5P_DEFAULT, + H5P_DEFAULT); + coldGroup = + H5Gcreate(matGroup, SP5::Depends::coldCurve, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + + 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); + // 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); + // currently unused + // status += mask.saveHDF(coldGroup, SP5::Fields::mask); + // 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}; + 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); + status += H5Gclose(lTGroup); + status += H5Gclose(coldGroup); + status += H5Gclose(matGroup); + + return status; +} + +herr_t saveAllMaterials(const std::string &savename, + const std::vector &filenames, bool printMetadata, + Verbosity eospacWarn) { + std::vector params; + std::vector matids; + std::unordered_map used_names; + std::unordered_set used_matids; + SesameMetadata metadata; + hid_t file; + herr_t status = H5_SUCCESS; + + for (auto const &filename : filenames) { + AddMaterials(params, matids, filename); + } + + std::cout << "Saving to file " << savename << std::endl; + file = H5Fcreate(savename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + + // singularity version + H5LTset_attribute_string(file, "/", "singularity_version", SINGULARITY_VERSION); + // log type. 0 for true, 1 for NQT1, 2 for NQT2, -1 for single precision true + int log_type = singularity::FastMath::Settings::log_type; + H5LTset_attribute_int(file, "/", SP5::logType, &log_type, 1); + + std::cout << "Processing " << matids.size() << " materials..." << std::endl; + + for (size_t i = 0; i < matids.size(); i++) { + int matid = matids[i]; + if (used_matids.count(matid) > 0) { + std::cerr << "...Duplicate matid " << matid << " detected. Skipping." << std::endl; + continue; + } + used_matids.insert(matid); + + std::cout << "..." << matid << std::endl; + + eosGetMetadata(matid, metadata, Verbosity::Debug); + if (printMetadata) std::cout << metadata << std::endl; + + std::string name = params[i].Get("name", metadata.name); + if (name == "-1" || name == "") { + std::string new_name = "material_" + std::to_string(i); + std::cerr << "...WARNING: no reasonable name found. " + << "Using a default name: " << new_name << std::endl; + name = new_name; + } + if (used_names.count(name) > 0) { + used_names[name] += 1; + std::string new_name = name + "_" + std::to_string(used_names[name]); + std::cerr << "...WARNING: Name " << name << " already used. " + << "Using name: " << new_name << std::endl; + name = new_name; + } else { + used_names[name] = 1; + } + + Bounds lRhoBounds, lTBounds, leBounds; + getMatBounds(i, matid, metadata, params[i], lRhoBounds, lTBounds, leBounds); + + if (eospacWarn == Verbosity::Debug) { + std::cout << "bounds for log(rho), log(T), log(sie) are:\n" + << lRhoBounds << lTBounds << leBounds << std::endl; + } + + 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; + } + } + + std::cout << "Cleaning up." << std::endl; + status += H5Fclose(file); + if (status != H5_SUCCESS) { + std::cerr << "WARNING: problem with HDf5" << std::endl; + } + 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) { + + // The "epsilon" shifts here are required to avoid eospac + // extrapolation errors at table bounds + constexpr Real TINY = std::numeric_limits::epsilon(); + auto TinyShift = [=](Real val, int sign) { + Real shift = std::abs(std::min(10 * val * TINY, TINY)); + return val + sign * shift; + }; + + Real rhoMin = params.Get("rhomin", TinyShift(metadata.rhoMin, 1)); + Real rhoMax = params.Get("rhomax", metadata.rhoMax); + Real TMin = params.Get("Tmin", TinyShift(metadata.TMin, 1)); + Real TMax = params.Get("Tmax", metadata.TMax); + Real sieMin = params.Get("siemin", TinyShift(metadata.sieMin, 1)); + Real sieMax = params.Get("siemax", metadata.sieMax); + + checkValInMatBounds(matid, "rhoMin", rhoMin, metadata.rhoMin, metadata.rhoMax); + checkValInMatBounds(matid, "rhoMax", rhoMax, metadata.rhoMin, metadata.rhoMax); + checkValInMatBounds(matid, "TMin", TMin, metadata.TMin, metadata.TMax); + checkValInMatBounds(matid, "TMax", TMax, metadata.TMin, metadata.TMax); + checkValInMatBounds(matid, "sieMin", sieMin, metadata.sieMin, metadata.sieMax); + checkValInMatBounds(matid, "sieMax", sieMax, metadata.sieMin, metadata.sieMax); + + Real shrinklRhoBounds = params.Get("shrinklRhoBounds", 0.); + Real shrinklTBounds = params.Get("shrinklTBounds", 0.); + Real shrinkleBounds = params.Get("shrinkleBounds", 0.); + + shrinklRhoBounds = std::min(1., std::max(shrinklRhoBounds, 0.)); + shrinklTBounds = std::min(1., std::max(shrinklTBounds, 0.)); + shrinkleBounds = std::min(1., std::max(shrinkleBounds, 0.)); + + if (shrinklRhoBounds > 0 && (params.Contains("rhomin") || params.Contains("rhomax"))) { + std::cerr << "WARNING [" << matid << "]: " + << "shrinklRhoBounds > 0 and rhomin or rhomax set" << std::endl; + } + if (shrinklTBounds > 0 && (params.Contains("Tmin") || params.Contains("Tmax"))) { + std::cerr << "WARNING [" << matid << "]: " + << "shrinklTBounds > 0 and Tmin or Tmax set" << std::endl; + } + if (shrinkleBounds > 0 && (params.Contains("siemin") || params.Contains("siemax"))) { + std::cerr << "WARNING [" << matid << "]: " + << "shrinkleBounds > 0 and siemin or siemax set" << std::endl; + } + + int ppdRho = params.Get("numrho/decade", PPD_DEFAULT_RHO); + int numRhoDefault = Bounds::getNumPointsFromPPD(rhoMin, rhoMax, ppdRho); + + int ppdT = params.Get("numT/decade", PPD_DEFAULT_T); + int numTDefault = Bounds::getNumPointsFromPPD(TMin, TMax, ppdT); + + int ppdSie = params.Get("numSie/decade", PPD_DEFAULT_T); + int numSieDefault = Bounds::getNumPointsFromPPD(sieMin, sieMax, ppdSie); + + int numRho = params.Get("numrho", numRhoDefault); + int numT = params.Get("numT", numTDefault); + int numSie = params.Get("numsie", numSieDefault); + + constexpr Real TAnchor = 298.15; + Real rhoAnchor = params.Get("rho_fine_center", metadata.normalDensity); + if (std::isnan(rhoAnchor) || rhoAnchor <= 0 || rhoAnchor > 1e8) { + std::cerr << "WARNING [" << matid << "] " + << "normal density ill defined. Setting it to a sensible default." + << std::endl; + rhoAnchor = 1; + } + + // Piecewise grids stuff + const bool piecewiseRho = params.Get("piecewiseRho", true); + const bool piecewiseT = params.Get("piecewiseT", true); + const bool piecewiseSie = params.Get("piecewiseSie", true); + + const Real ppd_factor_rho_lo = + params.Get("rhoCoarseFactorLo", COARSE_FACTOR_DEFAULT_RHO_LO); + const Real ppd_factor_rho_hi = + params.Get("rhoCoarseFactorHi", COARSE_FACTOR_DEFAULT_RHO_HI); + const Real ppd_factor_T = params.Get("TCoarseFactor", COARSE_FACTOR_DEFAULT_T); + const Real ppd_factor_sie = params.Get("sieCoarseFactor", COARSE_FACTOR_DEFAULT_T); + const Real rho_fine_diameter = + params.Get("rhoFineDiameterDecades", RHO_FINE_DIAMETER_DEFAULT); + const Real TSplitPoint = params.Get("TSplitPoint", T_SPLIT_POINT_DEFAULT); + const Real rho_fine_center = rhoAnchor; + + // These override the rho center/diameter settings + Real rho_fine_min = params.Get("rhoFineMin", -1); + Real rho_fine_max = params.Get("rhoFineMax", -1); + if (rho_fine_min * rho_fine_max < 0) { + std::cerr << "WARNING [" << matid << "]: " + << "Either rhoFineMin or rhoFineMax is set while the other is still unset." + << " Both must be set to be sensible. Ignoring." << std::endl; + rho_fine_min = rho_fine_max = -1; + } + + // Forces density and temperature to be in a region where an offset + // is not needed. This improves resolution at low densities and + // temperatures. + + // Extrapolation and other resolution tricks will be explored in the + // future. + if (rhoMin < STRICTLY_POS_MIN_RHO) rhoMin = STRICTLY_POS_MIN_RHO; + if (TMin < STRICTLY_POS_MIN_T) TMin = STRICTLY_POS_MIN_T; + + if (piecewiseRho) { + if (rho_fine_min > 0) { + lRhoBounds = Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, + rho_fine_min, rho_fine_max, ppdRho, ppd_factor_rho_lo, + ppd_factor_rho_hi, true, shrinklRhoBounds); + } else { + lRhoBounds = + Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, rho_fine_diameter, + ppdRho, ppd_factor_rho_lo, ppd_factor_rho_hi, true, shrinklRhoBounds); + } + } else { + lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); + } + if (piecewiseT) { + lTBounds = Bounds(Bounds::TwoGrids(), TMin, TMax, TAnchor, TSplitPoint, ppdT, + ppd_factor_T, true, shrinklTBounds); + } else { + lTBounds = Bounds(TMin, TMax, numT, true, shrinklTBounds, TAnchor); + } + if (piecewiseSie) { + // compute temperature as a reasonable anchor point + constexpr int NT = 1; + constexpr EOS_INTEGER nXYPairs = 2; + EOS_INTEGER tableHandle[NT]; + EOS_INTEGER tableType[NT] = {EOS_Ut_DT}; + EOS_REAL rho[2], T[2], sie[2], dx[2], dy[2]; + { + eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Ut_DT"}, Verbosity::Quiet); + EOS_INTEGER eospacEofRT = tableHandle[0]; + rho[0] = rho[1] = densityToSesame(rhoAnchor); + T[0] = temperatureToSesame(TAnchor); + T[1] = temperatureToSesame(TSplitPoint); + eosSafeInterpolate(&eospacEofRT, nXYPairs, rho, T, sie, dx, dy, "EofRT", + Verbosity::Quiet); + eosSafeDestroy(NT, tableHandle, Verbosity::Quiet); + } + const Real sieAnchor = sie[0]; + const Real sieSplitPoint = sie[1]; + leBounds = Bounds(Bounds::TwoGrids(), sieMin, sieMax, sieAnchor, sieSplitPoint, + ppdSie, ppd_factor_sie, true, shrinkleBounds); + } else { + leBounds = Bounds(sieMin, sieMax, numSie, true, shrinkleBounds); + } + + std::cout << "lRho bounds are\n" + << lRhoBounds << "lT bounds are\n" + << lTBounds << "lSie bounds are \n" + << leBounds << std::endl; + + return; +} + +bool checkValInMatBounds(int matid, const std::string &name, Real val, Real vmin, + Real vmax) { + if (val < vmin || val > vmax) { + std::cerr << "WARNING [" << matid << "]: " << name + << " out of sesame table bounds. Consider changing this.\n" + << "\t" << name << ", [bounds] = " << val << ", [" << vmin << ", " << vmax + << "]" << std::endl; + return false; + } + return true; +} diff --git a/ipcress2spiner/generate_files.hpp b/ipcress2spiner/generate_files.hpp new file mode 100644 index 0000000..543a60f --- /dev/null +++ b/ipcress2spiner/generate_files.hpp @@ -0,0 +1,72 @@ +//====================================================================== +// sesame2spiner tool for converting eospac to spiner +// Author: Jonah Miller (jonahm@lanl.gov) +// © 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 _SESAME2SPINER_GENERATE_FILES_HPP_ +#define _SESAME2SPINER_GENERATE_FILES_HPP_ + +#include +#include + +#include +#include + +#include + +#include "io_eospac.hpp" +#include "parser.hpp" + +using namespace EospacWrapper; + +constexpr int PPD_DEFAULT_RHO = 350; +constexpr int PPD_DEFAULT_T = 100; +constexpr Real STRICTLY_POS_MIN_RHO = 1e-8; +constexpr Real STRICTLY_POS_MIN_T = 1e-2; +constexpr Real COARSE_FACTOR_DEFAULT_RHO_LO = 3; +constexpr Real COARSE_FACTOR_DEFAULT_RHO_HI = 5; +constexpr Real COARSE_FACTOR_DEFAULT_T = 1.5; +constexpr Real RHO_FINE_DIAMETER_DEFAULT = 1.5; +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, 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); + +bool checkValInMatBounds(int matid, const std::string &name, Real val, Real vmin, + Real vmax); + +int getNumPointsFromPPD(Real min, Real max, int ppd); + +#endif // _SESAME2SPINER_GENERATE_FILES_HPP_ diff --git a/ipcress2spiner/io_eospac.cpp b/ipcress2spiner/io_eospac.cpp new file mode 100644 index 0000000..2329072 --- /dev/null +++ b/ipcress2spiner/io_eospac.cpp @@ -0,0 +1,424 @@ +//====================================================================== +// sesame2spiner tool for converting eospac to spiner +// Author: Jonah Miller (jonahm@lanl.gov) +// © 2021-2025. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//====================================================================== + +#include +#include +#include +#include + +#include + +#include "io_eospac.hpp" + +// TODO: more error checking of bounds? +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; + EOS_INTEGER tableHandle[NT]; + EOS_INTEGER eospacPofRT, eospacTofRE, eospacEofRT; + 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 + 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]; + + // Make sie faster moving index, for easy inversion + // TODO: is this right? + Ps.resize(rhos.size(), sies.size()); + Ps.setRange(0, leBounds.grid); + Ps.setRange(1, lRhoBounds.grid); + Ts.copyMetadata(Ps); + bMods.copyMetadata(Ps); + dPdRho.copyMetadata(Ps); + dPde.copyMetadata(Ps); + dTdRho.copyMetadata(Ps); + dTde.copyMetadata(Ps); + 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 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_pack[iflat]); + Ps(j, i) = pressureFromSesame(P_pack[iflat]); + bMods(j, i) = bulkModulusFromSesame(std::max(bMod, 0.0)); + 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 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; + EOS_INTEGER tableHandle[NT]; + EOS_INTEGER eospacPofRT, eospacTofRE, eospacEofRT; + 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; + makeInterpPoints(rhos, lRhoBounds); + makeInterpPoints(Ts, lTBounds); + + // Load tables + 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]; + + // Make T faster moving index, for easy inversion + Ps.resize(rhos.size(), Ts.size()); + Ps.setRange(0, lTBounds.grid); + Ps.setRange(1, lRhoBounds.grid); + sies.copyMetadata(Ps); + bMods.copyMetadata(Ps); + dPdRho.copyMetadata(Ps); + dPdE.copyMetadata(Ps); + dTdRho.copyMetadata(Ps); + dTde.copyMetadata(Ps); + dEdRho.copyMetadata(Ps); + dEdT.copyMetadata(Ps); + mask.copyMetadata(Ps); + + // 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 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[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); +} + +void eosColdCurves(int matid, const Bounds &lRhoBounds, DataBox &Ps, DataBox &sies, + DataBox &dPdRho, DataBox &dEdRho, DataBox &bMods, DataBox &mask, + Verbosity eospacWarn) { + using namespace EospacWrapper; + + constexpr int NT = 2; + EOS_INTEGER tableHandle[NT]; + EOS_INTEGER eospacPColdCurve, eospacSieColdCurve; + EOS_INTEGER tableType[NT] = {EOS_Pc_D, EOS_Uc_D}; + + // indep vars + 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); + eospacPColdCurve = tableHandle[0]; + eospacSieColdCurve = tableHandle[1]; + + // vars + Ps.resize(rhos.size()); + Ps.setRange(0, lRhoBounds.grid); + sies.copyMetadata(Ps); + dPdRho.copyMetadata(Ps); + dEdRho.copyMetadata(Ps); + bMods.copyMetadata(Ps); + mask.copyMetadata(Ps); + + 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_pack[i]); + sies(i) = sieFromSesame(sie_pack[i]); + bMods(i) = bulkModulusFromSesame(std::max(bMod, 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. + } +} + +void eosColdCurveMask(int matid, const Bounds &lRhoBounds, const int numSie, + const DataBox &sieColdCurve, DataBox &mask, Verbosity eospacWarn) { + using namespace EospacWrapper; + + constexpr int NT = 1; + constexpr EOS_INTEGER nXYPairs = 1; + EOS_INTEGER tableHandle[NT]; + EOS_INTEGER eospacTofRE; + EOS_INTEGER tableType[NT] = {EOS_T_DUt}; + + // Load tables + eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_T_DUt"}, eospacWarn); + eospacTofRE = tableHandle[0]; + + // Interpolatable vars + EOS_REAL var[1], dx[1], dy[1]; + + // rho is as expected. sie will be filled in the loop. + std::vector rhos; + makeInterpPoints(rhos, lRhoBounds); + + // vars + mask.resize(rhos.size(), numSie); + mask.setRange(1, lRhoBounds.grid); + + // sie + Bounds coldBounds(-1, 1, numSie, false); + mask.setRange(0, coldBounds.grid); + + // 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); + Real sieColdAbs = std::abs(sieCold); + Bounds sieBounds(sieCold - sieColdAbs, sieCold + sieColdAbs, numSie, false); + for (int i = 0; i < numSie; i++) { + Real sie = sieToSesame(sieBounds.grid.x(i)); + // T + bool no_errors = eosSafeInterpolate(&eospacTofRE, nXYPairs, &rho, &sie, var, dx, dy, + "TofRE", eospacWarn); + mask(j, i) = no_errors ? 1.0 : 0.0; + } + } + 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; + + 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++) { + 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/ipcress2spiner/io_eospac.hpp b/ipcress2spiner/io_eospac.hpp new file mode 100644 index 0000000..2a2c166 --- /dev/null +++ b/ipcress2spiner/io_eospac.hpp @@ -0,0 +1,97 @@ +//====================================================================== +// sesame2spiner tool for converting eospac to spiner +// Author: Jonah Miller (jonahm@lanl.gov) +// © 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 _SESAME2SPINER_IO_EOSPAC_HPP_ +#define _SESAME2SPINER_IO_EOSPAC_HPP_ + +#include +#include + +#include // eospac API + +#ifndef SINGULARITY_USE_SPINER_WITH_HDF5 +#error "HDF5 must be enabled" +#endif + +#include +#include +#include +#include +#include + +#include + +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 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 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, + Verbosity eospacWarn = Verbosity::Quiet); + +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); + +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/ipcress2spiner/main.cpp b/ipcress2spiner/main.cpp new file mode 100644 index 0000000..12c4390 --- /dev/null +++ b/ipcress2spiner/main.cpp @@ -0,0 +1,291 @@ +//====================================================================== +// ipcress2spiner tool for converting eospac to spiner +// Author: Alex Long (along@lanl.gov) +// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//====================================================================== + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#ifndef SPINER_USE_HDF +#error "HDF5 must be enabled" +#endif + +#include +#include +#include +#include +#include +#include + +//#include "generate_files.hpp" +//#include "io_eospac.hpp" +#include "parse_cli.hpp" + +/* + * Function prototypes for Gandolf C function equivalents + */ + +extern "C" void c_gmatids( char*, int*, int*, int*, int*); +extern "C" void c_gkeys( char*, int*,char*[], int*, int*, int*); +extern void c_gsizeof( char*, int*, char*, int*, int*); +extern "C" void c_gchgrids( char*, int*, int*, int*, int*, int*, int*, int*); +extern void c_ggetgray( char*, int*, char*,double*, int*, int*, + double*, int*, int*,double*, int*, int*, int*); +extern "C" void c_ggetmg( char*, int*, char*,double*, int*, int*, + double*, int*, int*, + double*, int*, int*,double*, int*, int*, int*); +extern void c_ggetdata( char*, int*, char*,double*, int*, int*, int*); +extern void c_ggetchar( char*, int*, char*,double*, int*, int*, int*); + +/* + * Function prototypes for Gandolf Fortran subroutines + */ + +extern void f_gmatids( char*, int*, int*, int*, int*); +extern void f_gkeys( char*, int*, char*, int*, int*, int*); +extern void f_gsizeof( char*, int*, char*, int*, int*); +extern void f_gchgrids( char*, int*, int*, int*, int*, int*, int*, int*); +extern void f_ggetgray( char*, int*, char*,double*, int*, int*,double*, + int*, int*,double*, int*, int*, int*); +extern void f_ggetmg( char*, int*, char*,double*, int*, int*,double*, + int*, int*,double*, int*, int*,double*, int*, int*, int*); +extern void f_ggetdata( char*, int*, char*,double*, int*, int*, int*); +extern void f_ggetchar( char*, int*, char*,double*, int*, int*, int*); + + + +//constexpr int H5_SUCCESS = 0; // older HDF5's declare this, just declare it here + +void print_gchrids_output(const std::string &filename, const int mat_ID, const int nt, const int nrho, const int nhnu, const int ngray, const int nmg) { + std::cout<<"---- Gandolf information for material "< &temperature_points, + const std::vector &density_points, + const std::vector &group_bounds, + const std::vector &ramg) { + + std::cout<<"Temperature points: "<(nmg); ++idata) { + std::cout< &T_grid, const std::vector &rho_grid, const std::vector &group_bounds, const std::vector &op_data, + double target_T, double target_rho, double target_hnu) { + + size_t n_T = T_grid.size(); + size_t n_rho = rho_grid.size(); + size_t n_groups = group_bounds.size()-1; + + auto T_L = std::distance(T_grid.begin(), std::lower_bound(T_grid.begin(), T_grid.end(), target_T)); + T_L = (target_T < T_grid[T_L]) ? T_L-1 : T_L; + T_L = (T_L == (T_grid.size()-1)) ? T_L-1: T_L; + auto T_R = T_L + 1; + auto log_T_L = log(T_grid[T_L]); + auto log_T_R = log(T_grid[T_R]); + auto log_T_target = log(target_T); + // fraction of T_2 to use + auto T_frac = (log_T_target - log_T_L) / (log_T_R - log_T_L); + + auto rho_L = std::distance(rho_grid.begin(), std::lower_bound(rho_grid.begin(), rho_grid.end(), target_rho)); + rho_L = (target_rho < rho_grid[rho_L]) ? rho_L-1 : rho_L; + rho_L = (rho_L == (rho_grid.size()-1)) ? rho_L-1: rho_L; + auto rho_R = rho_L + 1; + auto log_rho_L = log(rho_grid[rho_L]); + auto log_rho_R = log(rho_grid[rho_R]); + auto log_rho_target = log(target_rho); + // fraction of rho_2 to use + auto rho_frac = (log_rho_target - log_rho_L) / (log_rho_R - log_rho_L); + + auto group = std::distance(group_bounds.begin(), std::lower_bound(group_bounds.begin(), group_bounds.end(), target_hnu)); + group = (target_hnu < group_bounds[group]) ? group-1 : group; + group = (group == (group_bounds.size()-1)) ? group-1 : group; + + //std::cout<<"Targets--T: "< build_opacity_spiner_databox(const std::vector &temperature_points, const std::vector &density_points, const std::vector &group_bounds, const std::vector &ramg) { + + Spiner::RegularGrid1D new_temperature(temperature_points.front(), temperature_points.back(), temperature_points.size()); + Spiner::RegularGrid1D new_density(density_points.front(), density_points.back(), density_points.size()); + Spiner::RegularGrid1D new_group_bounds(group_bounds.front(), group_bounds.back(), 2*group_bounds.size()); + Spiner::DataBox opacity_databox(temperature_points.size(), density_points.size(), 2*group_bounds.size()); + opacity_databox.setRange(0, new_temperature); + opacity_databox.setRange(1, new_density); + opacity_databox.setRange(2, new_group_bounds); + + std::cout<<"T rho hnu opac(sq_cm/g)"< mat_ids(256, -1); + int kmats = static_cast(mat_ids.size()); + int nmats = -1; + c_gmatids(filename_char, mat_ids.data(), &kmats, &nmats, &ier); + + std::cout<<"nmats: "< keys; + for(int ikey = 0;ikey<32;++ikey) { + keys.push_back(new char(16)); + } + int kkeys = static_cast(keys.size()); + int nkeys = -1; + c_gkeys(filename_char, &mat_ID, keys.data(), &kkeys, &nkeys, &ier); + std::cout<<"nkeys: "< nice_keys; + for (int j=0;j temperature_points(nt, 0.0); + std::vector group_bounds(nhnu, 0.0); + std::vector density_points(nrho, 0.0); + int post_nt = nt; + int post_nrho = nrho; + int post_nhnu = nhnu; + int post_nmg = nmg; + + std::string ramg_keyword = "ramg"; + std::string rsmg_keyword = "rsmg"; + std::string pmg_keyword = "pmg"; + + std::array mg_opac_keywords{ramg_keyword, rsmg_keyword, pmg_keyword}; + + for (auto key : mg_opac_keywords) { + + std::cout<<"Interpolating and building databox for "< opacity_data(nmg, 0.0); + c_ggetmg(filename_char, &mat_ID, key.data(), + temperature_points.data(), &nt, &post_nt, + density_points.data(), &nrho, &post_nrho, + group_bounds.data(), &nhnu, &post_nhnu, + opacity_data.data(), &nmg, &post_nmg, &ier); + + auto spiner_opacity_databox = build_opacity_spiner_databox(temperature_points, density_points, group_bounds, opacity_data); + // save here with keyword? + } + } + + std::cout << "Done." << std::endl; + + return (status == H5_SUCCESS) ? 0 : 1; +} diff --git a/ipcress2spiner/parse_cli.cpp b/ipcress2spiner/parse_cli.cpp new file mode 100644 index 0000000..e4e189a --- /dev/null +++ b/ipcress2spiner/parse_cli.cpp @@ -0,0 +1,80 @@ +//====================================================================== +// sesame2spiner tool for converting eospac to spiner +// Author: Jonah Miller (jonahm@lanl.gov) +// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//====================================================================== + +#include +#include +#include +#include +#include +#include +#include + +//#include "io_eospac.hpp" +#include "parse_cli.hpp" + +void parseCLI(int argc, char *argv[], std::string &savename, + std::string &filename, bool &printMetadata, + std::string &helpMessage) { + + std::stringstream helpStream; + helpStream << "Usage: " << argv[0] + << "[-p] [-w] [-h] [-v] [-vv] [-d] [-s ] \n\n" + << "\t : input ipcress file\n" + << "\t-s : filename to save to. Defaults to " << DEFAULT_SAVENAME + << "\n" + << "\t-p: print metadata associated with materials " + << "in parameter files\n" + << "\t-v: print ipcress warnings\n" + << "\t-vv: print debug information\n" + << "\t-w: same as -v\n" + << "\t-d: same as -vv\n" + << "\t-h: print this message\n" + << "\n" + // ARL: Don't have input filenames + //<< "Several example input files:\n" + //<< EXAMPLESTRING << "\n" + << std::endl; + helpMessage = helpStream.str(); + + savename = DEFAULT_SAVENAME; + + if (argc < 2) { + std::cerr << helpMessage << std::endl; + std::exit(1); + } + if (argc == 2 && std::strcmp(argv[1], "-h") == 0) { + std::cout << helpMessage << std::endl; + std::exit(0); + } + for (int i = 1; i < argc; i++) { + if (std::strcmp(argv[i], "-h") == 0) { + std::cout << helpMessage << std::endl; + std::exit(0); + } else if (std::strcmp(argv[i], "-p") == 0) { + printMetadata = true; + //} else if ((std::strcmp(argv[i], "-w") == 0 || std::strcmp(argv[i], "-v") == 0) && + // eospacWarn == Verbosity::Quiet) { + // eospacWarn = Verbosity::Verbose; + //} else if ((std::strcmp(argv[i], "-d") == 0 || std::strcmp(argv[i], "-vv") == 0) && + // eospacWarn != Verbosity::Debug) { + // eospacWarn = Verbosity::Debug; + } else if (std::strcmp(argv[i], "-s") == 0) { + savename = argv[++i]; + } else { + filename = std::string(argv[i]); + } + } +} diff --git a/ipcress2spiner/parse_cli.hpp b/ipcress2spiner/parse_cli.hpp new file mode 100644 index 0000000..759960f --- /dev/null +++ b/ipcress2spiner/parse_cli.hpp @@ -0,0 +1,71 @@ +//====================================================================== +// sesame2spiner tool for converting eospac to spiner +// Author: Jonah Miller (jonahm@lanl.gov) +// © 2021-2023. 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 _SESAME2SPINER_PARSER_HPP_ +#define _SESAME2SPINER_PARSER_HPP_ + +//#include "io_eospac.hpp" +#include +#include + +const std::string DEFAULT_SAVENAME = "opacity.sp5"; +const std::string EXAMPLESTRING = R"( +# air.dat +# These are comments. +# The "#" character must be at the beginning of a line. +# only matid is required. All others override defaults. +matid = 5030 +name = air +# rho is in g/cm^3 +rhomin = 1e-2 +rhomax = 10 +numrho = 64 +# T is in Kelvin +Tmin = 252 +Tmax = 1e4 +numT = 32 +# sie is in erg/g +siemin = 1e12 +siemax = 1e16 +numsie = 32 + + +# titanium.dat +matid = 2961 +name = titanium +# These set the number of grid poitns per decade +# for each variable. The default is 50 points. +numrho/decade = 30 +numT/decade = 25 +numSie/decade = 15 + + +# steel.dat +matid=4272 +rhomin = 1e-2 +Tmin = 1 +# These shrink lograithm of bounds +# by a fraction of the total interval <= 1 +shrinklRhoBounds = 0.15 +shrinklTBounds = 0.15 +shrinkleBounds = 0.5 +)"; + +void parseCLI(int argc, char *argv[], std::string &savename, + std::string &filename, bool &printMetadata, + std::string &helpMessage); + +#endif // _SESAME2SPINER_PARSER_HPP_ diff --git a/ipcress2spiner/parser.cpp b/ipcress2spiner/parser.cpp new file mode 100644 index 0000000..7e8b056 --- /dev/null +++ b/ipcress2spiner/parser.cpp @@ -0,0 +1,161 @@ +//------------------------------------------------------------------------------ +// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "parse_cli.hpp" +#include "parser.hpp" + +Params::Params(const std::string &input_file) { + std::ifstream config_file(input_file); + if (config_file.is_open()) { + Parse(config_file); + } else { + throw std::runtime_error("Couldn't open config file " + input_file + "\n"); + } +} + +void Params::Parse(std::istream &s) { + std::string line; + while (getline(s, line)) { + line.erase(std::remove_if(line.begin(), line.end(), isspace), line.end()); + if (line[0] == '#' || line.empty()) continue; + auto delimiter_pos = line.find("="); + auto name = line.substr(0, delimiter_pos); + auto value = line.substr(delimiter_pos + 1); + // make sure there's no trailing comment + auto comment_pos = value.find("#"); + if (comment_pos != std::string::npos) { + value.erase(comment_pos); + } + params_[name] = value; + } +} + +void Params::Print(std::ostream &s) const { + s << "\nParams for " << params_.at("name") << "\n"; + for (const auto &pair : params_) { + s << pair.first << " = " << pair.second << "\n"; + } +} + +template <> +std::string Params::Get(const std::string &key) const { + return params_.at(key); +} + +template <> +int Params::Get(const std::string &key) const { + return std::stoi(params_.at(key)); +} + +template <> +double Params::Get(const std::string &key) const { + return std::stod(params_.at(key)); +} + +template <> +float Params::Get(const std::string &key) const { + return std::stof(params_.at(key)); +} + +template <> +bool Params::Get(const std::string &key) const { + std::string val = params_.at(key); + if ((val == "true") || (val == "1")) { + return true; + } else if ((val == "false") || (val == "0")) { + return false; + } else { + throw std::runtime_error("The value of " + key + " is not a boolean.\n"); + } +} + +void add_param(const Params &p, std::vector ¶ms, std::vector &matids, + const std::string &filename) { + + if (!p.Contains("matid")) { + if (p.Contains("name")) { + const auto &name = p.Get("name"); + std::cerr << "Material " << name << " in file " << filename << "is missing matid.\n" + << "Example input files:\n" + << EXAMPLESTRING << std::endl; + std::exit(1); + } else { + std::cerr << "A Material in file " << filename << "has no name.\n" + << "Example input files:\n" + << EXAMPLESTRING << std::endl; + std::exit(1); + } + } + matids.push_back(p.Get("matid")); + params.push_back(p); +} + +void parse_file(std::vector ¶ms, std::vector &matids, std::istream &s, + const std::string &filename) { + std::string line; + Params p; + size_t line_num = 1; + while (getline(s, line)) { + line.erase(std::remove_if(line.begin(), line.end(), isspace), line.end()); + if (line[0] == '#' || line.empty()) { + line_num++; + continue; + } + if (line[0] == '<') { + auto term = line.find('>'); + if (term == std::string::npos) { + std::cerr << "Missing closing > on line " << line_num << " of file " << filename + << "\n"; + std::exit(1); + } + + auto name = line.substr(1, term - 1); + if (!p.Empty()) add_param(p, params, matids, filename); + + p.Clear(); + p.Set("name", name); + } else { + auto delimiter_pos = line.find("="); + auto name = line.substr(0, delimiter_pos); + auto value = line.substr(delimiter_pos + 1); + // make sure there's no trailing comment + auto comment_pos = value.find("#"); + if (comment_pos != std::string::npos) { + value.erase(comment_pos); + } + p.Set(name, value); + } + line_num++; + } + if (!p.Empty()) add_param(p, params, matids, filename); +} + +void AddMaterials(std::vector ¶ms, std::vector &matids, + const std::string &input_file) { + std::ifstream config_file(input_file); + if (config_file.is_open()) { + parse_file(params, matids, config_file, input_file); + } else { + throw std::runtime_error("Couldn't open config file " + input_file + "\n"); + } +} \ No newline at end of file diff --git a/ipcress2spiner/parser.hpp b/ipcress2spiner/parser.hpp new file mode 100644 index 0000000..10b9a1e --- /dev/null +++ b/ipcress2spiner/parser.hpp @@ -0,0 +1,52 @@ +//------------------------------------------------------------------------------ +// © 2021-2023. 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 SESAME2SPINER_PARSER_HPP_ +#define SESAME2SPINER_PARSER_HPP_ + +#include +#include +#include +#include +#include + +// Parse a simple parameter file with +// "#" denoting comments. +class Params { + public: + Params() = default; + Params(const std::string &input_file); + Params(std::stringstream &input) { Parse(input); } + + bool Contains(const std::string &key) const { return params_.count(key); } + template + T Get(const std::string &key) const; + template + T Get(const std::string &key, T default_value) const { + return Contains(key) ? Get(key) : default_value; + } + void Print(std::ostream &s) const; + void Clear() { params_.clear(); } + bool Empty() const { return params_.size() == 0; } + void Set(const std::string &key, const std::string &val) { params_[key] = val; } + + private: + void Parse(std::istream &s); + std::unordered_map params_; +}; + +void AddMaterials(std::vector ¶ms, std::vector &matids, + const std::string &input_file); + +#endif // SESAME2SPINER_PARSER_HPP_ From aa28f6d59ddba43c314c0607f4d30a9bef33bbe9 Mon Sep 17 00:00:00 2001 From: "Alex R. Long" Date: Thu, 9 Jul 2026 11:22:08 -0600 Subject: [PATCH 2/5] Remove unused files --- ipcress2spiner/CMakeLists.txt | 6 - ipcress2spiner/generate_files.cpp | 457 ------------------------------ ipcress2spiner/generate_files.hpp | 72 ----- ipcress2spiner/io_eospac.cpp | 424 --------------------------- ipcress2spiner/io_eospac.hpp | 97 ------- ipcress2spiner/parser.cpp | 161 ----------- ipcress2spiner/parser.hpp | 52 ---- 7 files changed, 1269 deletions(-) delete mode 100644 ipcress2spiner/generate_files.cpp delete mode 100644 ipcress2spiner/generate_files.hpp delete mode 100644 ipcress2spiner/io_eospac.cpp delete mode 100644 ipcress2spiner/io_eospac.hpp delete mode 100644 ipcress2spiner/parser.cpp delete mode 100644 ipcress2spiner/parser.hpp diff --git a/ipcress2spiner/CMakeLists.txt b/ipcress2spiner/CMakeLists.txt index 0329aac..e6bff75 100644 --- a/ipcress2spiner/CMakeLists.txt +++ b/ipcress2spiner/CMakeLists.txt @@ -14,14 +14,8 @@ project(ipcress2spiner LANGUAGES CXX Fortran) add_executable(ipcress2spiner - #io_eospac.cpp - #io_eospac.hpp - #generate_files.cpp - #generate_files.hpp parse_cli.cpp parse_cli.hpp - #parser.cpp - #parser.hpp main.cpp ) diff --git a/ipcress2spiner/generate_files.cpp b/ipcress2spiner/generate_files.cpp deleted file mode 100644 index e5453f9..0000000 --- a/ipcress2spiner/generate_files.cpp +++ /dev/null @@ -1,457 +0,0 @@ -//====================================================================== -// sesame2spiner tool for converting eospac to spiner -// Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2025. Triad National Security, LLC. All rights reserved. This -// program was produced under U.S. Government contract 89233218CNA000001 -// for Los Alamos National Laboratory (LANL), which is operated by Triad -// National Security, LLC for the U.S. Department of Energy/National -// Nuclear Security Administration. All rights in the program are -// reserved by Triad National Security, LLC, and the U.S. Department of -// Energy/National Nuclear Security Administration. The Government is -// granted for itself and others acting on its behalf a nonexclusive, -// paid-up, irrevocable worldwide license in this material to reproduce, -// prepare derivative works, distribute copies to the public, perform -// publicly and display publicly, and to permit others to do so. -//====================================================================== - -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -#ifndef SPINER_USE_HDF -#error "HDF5 must be enabled" -#endif // SPINER_USE_HDF - -#include -#include -#include -#include -#include -#include -#include - -#include "generate_files.hpp" -#include "io_eospac.hpp" -#include "parse_cli.hpp" -#include "parser.hpp" - -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, const bool addSubtables, - Verbosity eospacWarn) { - - const int matid = metadata.matid; - std::string sMatid = std::to_string(matid); - - herr_t status = 0; - 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); - - // Dependent variables metadata - status += H5LTset_attribute_string(loc, sMatid.c_str(), SP5::Offsets::messageName, - SP5::Offsets::message); - status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Offsets::rho, - &lRhoBounds.offset, 1); - status += - H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Offsets::T, &lTBounds.offset, 1); - status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Offsets::sie, - &leBounds.offset, 1); - - // Material metadata - status += - H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Material::exchangeCoefficient, - &metadata.exchangeCoefficient, 1); - status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Material::meanAtomicMass, - &metadata.meanAtomicMass, 1); - status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Material::meanAtomicNumber, - &metadata.meanAtomicNumber, 1); - status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Material::solidBulkModulus, - &metadata.solidBulkModulus, 1); - status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Material::normalDensity, - &metadata.normalDensity, 1); - status += H5LTset_attribute_string(loc, sMatid.c_str(), SP5::Material::comments, - metadata.comments.c_str()); - status += H5LTset_attribute_int(loc, sMatid.c_str(), SP5::Material::matid, - &metadata.matid, 1); - status += H5LTset_attribute_string(loc, sMatid.c_str(), SP5::Material::name, - metadata.name.c_str()); - - lTGroup = H5Gcreate(matGroup, SP5::Depends::logRhoLogT, H5P_DEFAULT, H5P_DEFAULT, - H5P_DEFAULT); - - leGroup = H5Gcreate(matGroup, SP5::Depends::logRhoLogSie, H5P_DEFAULT, H5P_DEFAULT, - H5P_DEFAULT); - coldGroup = - H5Gcreate(matGroup, SP5::Depends::coldCurve, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - - 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); - // 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); - // currently unused - // status += mask.saveHDF(coldGroup, SP5::Fields::mask); - // 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}; - 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); - status += H5Gclose(lTGroup); - status += H5Gclose(coldGroup); - status += H5Gclose(matGroup); - - return status; -} - -herr_t saveAllMaterials(const std::string &savename, - const std::vector &filenames, bool printMetadata, - Verbosity eospacWarn) { - std::vector params; - std::vector matids; - std::unordered_map used_names; - std::unordered_set used_matids; - SesameMetadata metadata; - hid_t file; - herr_t status = H5_SUCCESS; - - for (auto const &filename : filenames) { - AddMaterials(params, matids, filename); - } - - std::cout << "Saving to file " << savename << std::endl; - file = H5Fcreate(savename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - - // singularity version - H5LTset_attribute_string(file, "/", "singularity_version", SINGULARITY_VERSION); - // log type. 0 for true, 1 for NQT1, 2 for NQT2, -1 for single precision true - int log_type = singularity::FastMath::Settings::log_type; - H5LTset_attribute_int(file, "/", SP5::logType, &log_type, 1); - - std::cout << "Processing " << matids.size() << " materials..." << std::endl; - - for (size_t i = 0; i < matids.size(); i++) { - int matid = matids[i]; - if (used_matids.count(matid) > 0) { - std::cerr << "...Duplicate matid " << matid << " detected. Skipping." << std::endl; - continue; - } - used_matids.insert(matid); - - std::cout << "..." << matid << std::endl; - - eosGetMetadata(matid, metadata, Verbosity::Debug); - if (printMetadata) std::cout << metadata << std::endl; - - std::string name = params[i].Get("name", metadata.name); - if (name == "-1" || name == "") { - std::string new_name = "material_" + std::to_string(i); - std::cerr << "...WARNING: no reasonable name found. " - << "Using a default name: " << new_name << std::endl; - name = new_name; - } - if (used_names.count(name) > 0) { - used_names[name] += 1; - std::string new_name = name + "_" + std::to_string(used_names[name]); - std::cerr << "...WARNING: Name " << name << " already used. " - << "Using name: " << new_name << std::endl; - name = new_name; - } else { - used_names[name] = 1; - } - - Bounds lRhoBounds, lTBounds, leBounds; - getMatBounds(i, matid, metadata, params[i], lRhoBounds, lTBounds, leBounds); - - if (eospacWarn == Verbosity::Debug) { - std::cout << "bounds for log(rho), log(T), log(sie) are:\n" - << lRhoBounds << lTBounds << leBounds << std::endl; - } - - 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; - } - } - - std::cout << "Cleaning up." << std::endl; - status += H5Fclose(file); - if (status != H5_SUCCESS) { - std::cerr << "WARNING: problem with HDf5" << std::endl; - } - 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) { - - // The "epsilon" shifts here are required to avoid eospac - // extrapolation errors at table bounds - constexpr Real TINY = std::numeric_limits::epsilon(); - auto TinyShift = [=](Real val, int sign) { - Real shift = std::abs(std::min(10 * val * TINY, TINY)); - return val + sign * shift; - }; - - Real rhoMin = params.Get("rhomin", TinyShift(metadata.rhoMin, 1)); - Real rhoMax = params.Get("rhomax", metadata.rhoMax); - Real TMin = params.Get("Tmin", TinyShift(metadata.TMin, 1)); - Real TMax = params.Get("Tmax", metadata.TMax); - Real sieMin = params.Get("siemin", TinyShift(metadata.sieMin, 1)); - Real sieMax = params.Get("siemax", metadata.sieMax); - - checkValInMatBounds(matid, "rhoMin", rhoMin, metadata.rhoMin, metadata.rhoMax); - checkValInMatBounds(matid, "rhoMax", rhoMax, metadata.rhoMin, metadata.rhoMax); - checkValInMatBounds(matid, "TMin", TMin, metadata.TMin, metadata.TMax); - checkValInMatBounds(matid, "TMax", TMax, metadata.TMin, metadata.TMax); - checkValInMatBounds(matid, "sieMin", sieMin, metadata.sieMin, metadata.sieMax); - checkValInMatBounds(matid, "sieMax", sieMax, metadata.sieMin, metadata.sieMax); - - Real shrinklRhoBounds = params.Get("shrinklRhoBounds", 0.); - Real shrinklTBounds = params.Get("shrinklTBounds", 0.); - Real shrinkleBounds = params.Get("shrinkleBounds", 0.); - - shrinklRhoBounds = std::min(1., std::max(shrinklRhoBounds, 0.)); - shrinklTBounds = std::min(1., std::max(shrinklTBounds, 0.)); - shrinkleBounds = std::min(1., std::max(shrinkleBounds, 0.)); - - if (shrinklRhoBounds > 0 && (params.Contains("rhomin") || params.Contains("rhomax"))) { - std::cerr << "WARNING [" << matid << "]: " - << "shrinklRhoBounds > 0 and rhomin or rhomax set" << std::endl; - } - if (shrinklTBounds > 0 && (params.Contains("Tmin") || params.Contains("Tmax"))) { - std::cerr << "WARNING [" << matid << "]: " - << "shrinklTBounds > 0 and Tmin or Tmax set" << std::endl; - } - if (shrinkleBounds > 0 && (params.Contains("siemin") || params.Contains("siemax"))) { - std::cerr << "WARNING [" << matid << "]: " - << "shrinkleBounds > 0 and siemin or siemax set" << std::endl; - } - - int ppdRho = params.Get("numrho/decade", PPD_DEFAULT_RHO); - int numRhoDefault = Bounds::getNumPointsFromPPD(rhoMin, rhoMax, ppdRho); - - int ppdT = params.Get("numT/decade", PPD_DEFAULT_T); - int numTDefault = Bounds::getNumPointsFromPPD(TMin, TMax, ppdT); - - int ppdSie = params.Get("numSie/decade", PPD_DEFAULT_T); - int numSieDefault = Bounds::getNumPointsFromPPD(sieMin, sieMax, ppdSie); - - int numRho = params.Get("numrho", numRhoDefault); - int numT = params.Get("numT", numTDefault); - int numSie = params.Get("numsie", numSieDefault); - - constexpr Real TAnchor = 298.15; - Real rhoAnchor = params.Get("rho_fine_center", metadata.normalDensity); - if (std::isnan(rhoAnchor) || rhoAnchor <= 0 || rhoAnchor > 1e8) { - std::cerr << "WARNING [" << matid << "] " - << "normal density ill defined. Setting it to a sensible default." - << std::endl; - rhoAnchor = 1; - } - - // Piecewise grids stuff - const bool piecewiseRho = params.Get("piecewiseRho", true); - const bool piecewiseT = params.Get("piecewiseT", true); - const bool piecewiseSie = params.Get("piecewiseSie", true); - - const Real ppd_factor_rho_lo = - params.Get("rhoCoarseFactorLo", COARSE_FACTOR_DEFAULT_RHO_LO); - const Real ppd_factor_rho_hi = - params.Get("rhoCoarseFactorHi", COARSE_FACTOR_DEFAULT_RHO_HI); - const Real ppd_factor_T = params.Get("TCoarseFactor", COARSE_FACTOR_DEFAULT_T); - const Real ppd_factor_sie = params.Get("sieCoarseFactor", COARSE_FACTOR_DEFAULT_T); - const Real rho_fine_diameter = - params.Get("rhoFineDiameterDecades", RHO_FINE_DIAMETER_DEFAULT); - const Real TSplitPoint = params.Get("TSplitPoint", T_SPLIT_POINT_DEFAULT); - const Real rho_fine_center = rhoAnchor; - - // These override the rho center/diameter settings - Real rho_fine_min = params.Get("rhoFineMin", -1); - Real rho_fine_max = params.Get("rhoFineMax", -1); - if (rho_fine_min * rho_fine_max < 0) { - std::cerr << "WARNING [" << matid << "]: " - << "Either rhoFineMin or rhoFineMax is set while the other is still unset." - << " Both must be set to be sensible. Ignoring." << std::endl; - rho_fine_min = rho_fine_max = -1; - } - - // Forces density and temperature to be in a region where an offset - // is not needed. This improves resolution at low densities and - // temperatures. - - // Extrapolation and other resolution tricks will be explored in the - // future. - if (rhoMin < STRICTLY_POS_MIN_RHO) rhoMin = STRICTLY_POS_MIN_RHO; - if (TMin < STRICTLY_POS_MIN_T) TMin = STRICTLY_POS_MIN_T; - - if (piecewiseRho) { - if (rho_fine_min > 0) { - lRhoBounds = Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, - rho_fine_min, rho_fine_max, ppdRho, ppd_factor_rho_lo, - ppd_factor_rho_hi, true, shrinklRhoBounds); - } else { - lRhoBounds = - Bounds(Bounds::ThreeGrids(), rhoMin, rhoMax, rho_fine_center, rho_fine_diameter, - ppdRho, ppd_factor_rho_lo, ppd_factor_rho_hi, true, shrinklRhoBounds); - } - } else { - lRhoBounds = Bounds(rhoMin, rhoMax, numRho, true, shrinklRhoBounds, rhoAnchor); - } - if (piecewiseT) { - lTBounds = Bounds(Bounds::TwoGrids(), TMin, TMax, TAnchor, TSplitPoint, ppdT, - ppd_factor_T, true, shrinklTBounds); - } else { - lTBounds = Bounds(TMin, TMax, numT, true, shrinklTBounds, TAnchor); - } - if (piecewiseSie) { - // compute temperature as a reasonable anchor point - constexpr int NT = 1; - constexpr EOS_INTEGER nXYPairs = 2; - EOS_INTEGER tableHandle[NT]; - EOS_INTEGER tableType[NT] = {EOS_Ut_DT}; - EOS_REAL rho[2], T[2], sie[2], dx[2], dy[2]; - { - eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Ut_DT"}, Verbosity::Quiet); - EOS_INTEGER eospacEofRT = tableHandle[0]; - rho[0] = rho[1] = densityToSesame(rhoAnchor); - T[0] = temperatureToSesame(TAnchor); - T[1] = temperatureToSesame(TSplitPoint); - eosSafeInterpolate(&eospacEofRT, nXYPairs, rho, T, sie, dx, dy, "EofRT", - Verbosity::Quiet); - eosSafeDestroy(NT, tableHandle, Verbosity::Quiet); - } - const Real sieAnchor = sie[0]; - const Real sieSplitPoint = sie[1]; - leBounds = Bounds(Bounds::TwoGrids(), sieMin, sieMax, sieAnchor, sieSplitPoint, - ppdSie, ppd_factor_sie, true, shrinkleBounds); - } else { - leBounds = Bounds(sieMin, sieMax, numSie, true, shrinkleBounds); - } - - std::cout << "lRho bounds are\n" - << lRhoBounds << "lT bounds are\n" - << lTBounds << "lSie bounds are \n" - << leBounds << std::endl; - - return; -} - -bool checkValInMatBounds(int matid, const std::string &name, Real val, Real vmin, - Real vmax) { - if (val < vmin || val > vmax) { - std::cerr << "WARNING [" << matid << "]: " << name - << " out of sesame table bounds. Consider changing this.\n" - << "\t" << name << ", [bounds] = " << val << ", [" << vmin << ", " << vmax - << "]" << std::endl; - return false; - } - return true; -} diff --git a/ipcress2spiner/generate_files.hpp b/ipcress2spiner/generate_files.hpp deleted file mode 100644 index 543a60f..0000000 --- a/ipcress2spiner/generate_files.hpp +++ /dev/null @@ -1,72 +0,0 @@ -//====================================================================== -// sesame2spiner tool for converting eospac to spiner -// Author: Jonah Miller (jonahm@lanl.gov) -// © 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 _SESAME2SPINER_GENERATE_FILES_HPP_ -#define _SESAME2SPINER_GENERATE_FILES_HPP_ - -#include -#include - -#include -#include - -#include - -#include "io_eospac.hpp" -#include "parser.hpp" - -using namespace EospacWrapper; - -constexpr int PPD_DEFAULT_RHO = 350; -constexpr int PPD_DEFAULT_T = 100; -constexpr Real STRICTLY_POS_MIN_RHO = 1e-8; -constexpr Real STRICTLY_POS_MIN_T = 1e-2; -constexpr Real COARSE_FACTOR_DEFAULT_RHO_LO = 3; -constexpr Real COARSE_FACTOR_DEFAULT_RHO_HI = 5; -constexpr Real COARSE_FACTOR_DEFAULT_T = 1.5; -constexpr Real RHO_FINE_DIAMETER_DEFAULT = 1.5; -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, 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); - -bool checkValInMatBounds(int matid, const std::string &name, Real val, Real vmin, - Real vmax); - -int getNumPointsFromPPD(Real min, Real max, int ppd); - -#endif // _SESAME2SPINER_GENERATE_FILES_HPP_ diff --git a/ipcress2spiner/io_eospac.cpp b/ipcress2spiner/io_eospac.cpp deleted file mode 100644 index 2329072..0000000 --- a/ipcress2spiner/io_eospac.cpp +++ /dev/null @@ -1,424 +0,0 @@ -//====================================================================== -// sesame2spiner tool for converting eospac to spiner -// Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2025. Triad National Security, LLC. All rights reserved. This -// program was produced under U.S. Government contract 89233218CNA000001 -// for Los Alamos National Laboratory (LANL), which is operated by Triad -// National Security, LLC for the U.S. Department of Energy/National -// Nuclear Security Administration. All rights in the program are -// reserved by Triad National Security, LLC, and the U.S. Department of -// Energy/National Nuclear Security Administration. The Government is -// granted for itself and others acting on its behalf a nonexclusive, -// paid-up, irrevocable worldwide license in this material to reproduce, -// prepare derivative works, distribute copies to the public, perform -// publicly and display publicly, and to permit others to do so. -//====================================================================== - -#include -#include -#include -#include - -#include - -#include "io_eospac.hpp" - -// TODO: more error checking of bounds? -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; - EOS_INTEGER tableHandle[NT]; - EOS_INTEGER eospacPofRT, eospacTofRE, eospacEofRT; - 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 - 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]; - - // Make sie faster moving index, for easy inversion - // TODO: is this right? - Ps.resize(rhos.size(), sies.size()); - Ps.setRange(0, leBounds.grid); - Ps.setRange(1, lRhoBounds.grid); - Ts.copyMetadata(Ps); - bMods.copyMetadata(Ps); - dPdRho.copyMetadata(Ps); - dPde.copyMetadata(Ps); - dTdRho.copyMetadata(Ps); - dTde.copyMetadata(Ps); - 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 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_pack[iflat]); - Ps(j, i) = pressureFromSesame(P_pack[iflat]); - bMods(j, i) = bulkModulusFromSesame(std::max(bMod, 0.0)); - 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 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; - EOS_INTEGER tableHandle[NT]; - EOS_INTEGER eospacPofRT, eospacTofRE, eospacEofRT; - 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; - makeInterpPoints(rhos, lRhoBounds); - makeInterpPoints(Ts, lTBounds); - - // Load tables - 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]; - - // Make T faster moving index, for easy inversion - Ps.resize(rhos.size(), Ts.size()); - Ps.setRange(0, lTBounds.grid); - Ps.setRange(1, lRhoBounds.grid); - sies.copyMetadata(Ps); - bMods.copyMetadata(Ps); - dPdRho.copyMetadata(Ps); - dPdE.copyMetadata(Ps); - dTdRho.copyMetadata(Ps); - dTde.copyMetadata(Ps); - dEdRho.copyMetadata(Ps); - dEdT.copyMetadata(Ps); - mask.copyMetadata(Ps); - - // 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 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[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); -} - -void eosColdCurves(int matid, const Bounds &lRhoBounds, DataBox &Ps, DataBox &sies, - DataBox &dPdRho, DataBox &dEdRho, DataBox &bMods, DataBox &mask, - Verbosity eospacWarn) { - using namespace EospacWrapper; - - constexpr int NT = 2; - EOS_INTEGER tableHandle[NT]; - EOS_INTEGER eospacPColdCurve, eospacSieColdCurve; - EOS_INTEGER tableType[NT] = {EOS_Pc_D, EOS_Uc_D}; - - // indep vars - 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); - eospacPColdCurve = tableHandle[0]; - eospacSieColdCurve = tableHandle[1]; - - // vars - Ps.resize(rhos.size()); - Ps.setRange(0, lRhoBounds.grid); - sies.copyMetadata(Ps); - dPdRho.copyMetadata(Ps); - dEdRho.copyMetadata(Ps); - bMods.copyMetadata(Ps); - mask.copyMetadata(Ps); - - 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_pack[i]); - sies(i) = sieFromSesame(sie_pack[i]); - bMods(i) = bulkModulusFromSesame(std::max(bMod, 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. - } -} - -void eosColdCurveMask(int matid, const Bounds &lRhoBounds, const int numSie, - const DataBox &sieColdCurve, DataBox &mask, Verbosity eospacWarn) { - using namespace EospacWrapper; - - constexpr int NT = 1; - constexpr EOS_INTEGER nXYPairs = 1; - EOS_INTEGER tableHandle[NT]; - EOS_INTEGER eospacTofRE; - EOS_INTEGER tableType[NT] = {EOS_T_DUt}; - - // Load tables - eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_T_DUt"}, eospacWarn); - eospacTofRE = tableHandle[0]; - - // Interpolatable vars - EOS_REAL var[1], dx[1], dy[1]; - - // rho is as expected. sie will be filled in the loop. - std::vector rhos; - makeInterpPoints(rhos, lRhoBounds); - - // vars - mask.resize(rhos.size(), numSie); - mask.setRange(1, lRhoBounds.grid); - - // sie - Bounds coldBounds(-1, 1, numSie, false); - mask.setRange(0, coldBounds.grid); - - // 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); - Real sieColdAbs = std::abs(sieCold); - Bounds sieBounds(sieCold - sieColdAbs, sieCold + sieColdAbs, numSie, false); - for (int i = 0; i < numSie; i++) { - Real sie = sieToSesame(sieBounds.grid.x(i)); - // T - bool no_errors = eosSafeInterpolate(&eospacTofRE, nXYPairs, &rho, &sie, var, dx, dy, - "TofRE", eospacWarn); - mask(j, i) = no_errors ? 1.0 : 0.0; - } - } - 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; - - 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++) { - 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/ipcress2spiner/io_eospac.hpp b/ipcress2spiner/io_eospac.hpp deleted file mode 100644 index 2a2c166..0000000 --- a/ipcress2spiner/io_eospac.hpp +++ /dev/null @@ -1,97 +0,0 @@ -//====================================================================== -// sesame2spiner tool for converting eospac to spiner -// Author: Jonah Miller (jonahm@lanl.gov) -// © 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 _SESAME2SPINER_IO_EOSPAC_HPP_ -#define _SESAME2SPINER_IO_EOSPAC_HPP_ - -#include -#include - -#include // eospac API - -#ifndef SINGULARITY_USE_SPINER_WITH_HDF5 -#error "HDF5 must be enabled" -#endif - -#include -#include -#include -#include -#include - -#include - -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 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 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, - Verbosity eospacWarn = Verbosity::Quiet); - -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); - -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/ipcress2spiner/parser.cpp b/ipcress2spiner/parser.cpp deleted file mode 100644 index 7e8b056..0000000 --- a/ipcress2spiner/parser.cpp +++ /dev/null @@ -1,161 +0,0 @@ -//------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This -// program was produced under U.S. Government contract 89233218CNA000001 -// for Los Alamos National Laboratory (LANL), which is operated by Triad -// National Security, LLC for the U.S. Department of Energy/National -// Nuclear Security Administration. All rights in the program are -// reserved by Triad National Security, LLC, and the U.S. Department of -// Energy/National Nuclear Security Administration. The Government is -// granted for itself and others acting on its behalf a nonexclusive, -// paid-up, irrevocable worldwide license in this material to reproduce, -// prepare derivative works, distribute copies to the public, perform -// publicly and display publicly, and to permit others to do so. -//------------------------------------------------------------------------------ - -#include -#include -#include -#include -#include -#include -#include -#include - -#include "parse_cli.hpp" -#include "parser.hpp" - -Params::Params(const std::string &input_file) { - std::ifstream config_file(input_file); - if (config_file.is_open()) { - Parse(config_file); - } else { - throw std::runtime_error("Couldn't open config file " + input_file + "\n"); - } -} - -void Params::Parse(std::istream &s) { - std::string line; - while (getline(s, line)) { - line.erase(std::remove_if(line.begin(), line.end(), isspace), line.end()); - if (line[0] == '#' || line.empty()) continue; - auto delimiter_pos = line.find("="); - auto name = line.substr(0, delimiter_pos); - auto value = line.substr(delimiter_pos + 1); - // make sure there's no trailing comment - auto comment_pos = value.find("#"); - if (comment_pos != std::string::npos) { - value.erase(comment_pos); - } - params_[name] = value; - } -} - -void Params::Print(std::ostream &s) const { - s << "\nParams for " << params_.at("name") << "\n"; - for (const auto &pair : params_) { - s << pair.first << " = " << pair.second << "\n"; - } -} - -template <> -std::string Params::Get(const std::string &key) const { - return params_.at(key); -} - -template <> -int Params::Get(const std::string &key) const { - return std::stoi(params_.at(key)); -} - -template <> -double Params::Get(const std::string &key) const { - return std::stod(params_.at(key)); -} - -template <> -float Params::Get(const std::string &key) const { - return std::stof(params_.at(key)); -} - -template <> -bool Params::Get(const std::string &key) const { - std::string val = params_.at(key); - if ((val == "true") || (val == "1")) { - return true; - } else if ((val == "false") || (val == "0")) { - return false; - } else { - throw std::runtime_error("The value of " + key + " is not a boolean.\n"); - } -} - -void add_param(const Params &p, std::vector ¶ms, std::vector &matids, - const std::string &filename) { - - if (!p.Contains("matid")) { - if (p.Contains("name")) { - const auto &name = p.Get("name"); - std::cerr << "Material " << name << " in file " << filename << "is missing matid.\n" - << "Example input files:\n" - << EXAMPLESTRING << std::endl; - std::exit(1); - } else { - std::cerr << "A Material in file " << filename << "has no name.\n" - << "Example input files:\n" - << EXAMPLESTRING << std::endl; - std::exit(1); - } - } - matids.push_back(p.Get("matid")); - params.push_back(p); -} - -void parse_file(std::vector ¶ms, std::vector &matids, std::istream &s, - const std::string &filename) { - std::string line; - Params p; - size_t line_num = 1; - while (getline(s, line)) { - line.erase(std::remove_if(line.begin(), line.end(), isspace), line.end()); - if (line[0] == '#' || line.empty()) { - line_num++; - continue; - } - if (line[0] == '<') { - auto term = line.find('>'); - if (term == std::string::npos) { - std::cerr << "Missing closing > on line " << line_num << " of file " << filename - << "\n"; - std::exit(1); - } - - auto name = line.substr(1, term - 1); - if (!p.Empty()) add_param(p, params, matids, filename); - - p.Clear(); - p.Set("name", name); - } else { - auto delimiter_pos = line.find("="); - auto name = line.substr(0, delimiter_pos); - auto value = line.substr(delimiter_pos + 1); - // make sure there's no trailing comment - auto comment_pos = value.find("#"); - if (comment_pos != std::string::npos) { - value.erase(comment_pos); - } - p.Set(name, value); - } - line_num++; - } - if (!p.Empty()) add_param(p, params, matids, filename); -} - -void AddMaterials(std::vector ¶ms, std::vector &matids, - const std::string &input_file) { - std::ifstream config_file(input_file); - if (config_file.is_open()) { - parse_file(params, matids, config_file, input_file); - } else { - throw std::runtime_error("Couldn't open config file " + input_file + "\n"); - } -} \ No newline at end of file diff --git a/ipcress2spiner/parser.hpp b/ipcress2spiner/parser.hpp deleted file mode 100644 index 10b9a1e..0000000 --- a/ipcress2spiner/parser.hpp +++ /dev/null @@ -1,52 +0,0 @@ -//------------------------------------------------------------------------------ -// © 2021-2023. 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 SESAME2SPINER_PARSER_HPP_ -#define SESAME2SPINER_PARSER_HPP_ - -#include -#include -#include -#include -#include - -// Parse a simple parameter file with -// "#" denoting comments. -class Params { - public: - Params() = default; - Params(const std::string &input_file); - Params(std::stringstream &input) { Parse(input); } - - bool Contains(const std::string &key) const { return params_.count(key); } - template - T Get(const std::string &key) const; - template - T Get(const std::string &key, T default_value) const { - return Contains(key) ? Get(key) : default_value; - } - void Print(std::ostream &s) const; - void Clear() { params_.clear(); } - bool Empty() const { return params_.size() == 0; } - void Set(const std::string &key, const std::string &val) { params_[key] = val; } - - private: - void Parse(std::istream &s); - std::unordered_map params_; -}; - -void AddMaterials(std::vector ¶ms, std::vector &matids, - const std::string &input_file); - -#endif // SESAME2SPINER_PARSER_HPP_ From cba114f0a14c269f31987fa2f80220d5af63789e Mon Sep 17 00:00:00 2001 From: "Alex R. Long" Date: Thu, 30 Jul 2026 22:39:35 -0600 Subject: [PATCH 3/5] Add HDF5 writing routines + Mkae Spiner databoxes for each opacity type + Call HDF5 to make converted file from IPCRESS data + Update to-do list --- ipcress2spiner/CMakeLists.txt | 10 +- ipcress2spiner/main.cpp | 195 ++++++++++++++-------------------- ipcress2spiner/parse_cli.hpp | 15 ++- 3 files changed, 96 insertions(+), 124 deletions(-) diff --git a/ipcress2spiner/CMakeLists.txt b/ipcress2spiner/CMakeLists.txt index e6bff75..bd015be 100644 --- a/ipcress2spiner/CMakeLists.txt +++ b/ipcress2spiner/CMakeLists.txt @@ -11,11 +11,13 @@ # prepare derivative works, distribute copies to the public, perform # publicly and display publicly, and to permit others to do so. #------------------------------------------------------------------------------ -project(ipcress2spiner LANGUAGES CXX Fortran) +project(ipcress2spiner VERSION 0.8.0.0 LANGUAGES CXX Fortran) add_executable(ipcress2spiner parse_cli.cpp parse_cli.hpp + make_spiner_databox.hpp + interpolate_opacity.hpp main.cpp ) @@ -38,8 +40,12 @@ target_link_libraries(ipcress2spiner ${GANDOLF_LIBRARY} PortsofCall::PortsofCall ) -set_property(TARGET ipcress2spiner PROPERTY LINKER_LANGUAGE Fortran) +set_property(TARGET ipcress2spiner PROPERTY LINKER_LANGUAGE Fortran) +# Define the full version as a string macro +target_compile_definitions(ipcress2spiner PRIVATE + IPCRESS2SPINER_VERSION=\"${PROJECT_VERSION}\" +) install(TARGETS ipcress2spiner DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/ipcress2spiner/main.cpp b/ipcress2spiner/main.cpp index 12c4390..4d0844e 100644 --- a/ipcress2spiner/main.cpp +++ b/ipcress2spiner/main.cpp @@ -1,7 +1,7 @@ //====================================================================== -// ipcress2spiner tool for converting eospac to spiner +// ipcress2spiner tool for converting ipcresse to spiner // Author: Alex Long (along@lanl.gov) -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 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 @@ -37,19 +37,27 @@ #include #include -//#include "generate_files.hpp" -//#include "io_eospac.hpp" +#include "generate_files.hpp" +#include "interpolate_opacity.hpp" #include "parse_cli.hpp" - -/* - * Function prototypes for Gandolf C function equivalents - */ - +#include "make_spiner_databox.hpp" + +// ARL: My to-do list +// \TODO Thread log interpolation points and log interpolation mode options through command line +// \TODO Add gray opacties +// \TODO Add ionization state? +// \TODO clang-format +// \TODO Add tests +// \TODO Compare interpolation to draco's routines +// \TODO Delete char fields from new[] + +// Function prototypes for Gandolf C function equivalents, prepend "C" to prevent name mangling +// gandolf.h does not expose the interface to these functions, which is why they are declared here extern "C" void c_gmatids( char*, int*, int*, int*, int*); extern "C" void c_gkeys( char*, int*,char*[], int*, int*, int*); -extern void c_gsizeof( char*, int*, char*, int*, int*); +extern "C" void c_gsizeof( char*, int*, char*, int*, int*); extern "C" void c_gchgrids( char*, int*, int*, int*, int*, int*, int*, int*); -extern void c_ggetgray( char*, int*, char*,double*, int*, int*, +extern "C" void c_ggetgray( char*, int*, char*,double*, int*, int*, double*, int*, int*,double*, int*, int*, int*); extern "C" void c_ggetmg( char*, int*, char*,double*, int*, int*, double*, int*, int*, @@ -57,10 +65,8 @@ extern "C" void c_ggetmg( char*, int*, char*,double*, int*, int*, extern void c_ggetdata( char*, int*, char*,double*, int*, int*, int*); extern void c_ggetchar( char*, int*, char*,double*, int*, int*, int*); -/* - * Function prototypes for Gandolf Fortran subroutines - */ - +// Function prototypes for Gandolf Fortran subroutines +// gandolf.h does not expose the interface to these functions, which is why they are declared here extern void f_gmatids( char*, int*, int*, int*, int*); extern void f_gkeys( char*, int*, char*, int*, int*, int*); extern void f_gsizeof( char*, int*, char*, int*, int*); @@ -72,8 +78,6 @@ extern void f_ggetmg( char*, int*, char*,double*, int*, int*,double*, extern void f_ggetdata( char*, int*, char*,double*, int*, int*, int*); extern void f_ggetchar( char*, int*, char*,double*, int*, int*, int*); - - //constexpr int H5_SUCCESS = 0; // older HDF5's declare this, just declare it here void print_gchrids_output(const std::string &filename, const int mat_ID, const int nt, const int nrho, const int nhnu, const int ngray, const int nmg) { @@ -109,98 +113,10 @@ void print_ggetmg_output(const int nmg, const int post_ramg_nmg, } } -double interpolate_mg_opacity_data(const std::vector &T_grid, const std::vector &rho_grid, const std::vector &group_bounds, const std::vector &op_data, - double target_T, double target_rho, double target_hnu) { - - size_t n_T = T_grid.size(); - size_t n_rho = rho_grid.size(); - size_t n_groups = group_bounds.size()-1; - - auto T_L = std::distance(T_grid.begin(), std::lower_bound(T_grid.begin(), T_grid.end(), target_T)); - T_L = (target_T < T_grid[T_L]) ? T_L-1 : T_L; - T_L = (T_L == (T_grid.size()-1)) ? T_L-1: T_L; - auto T_R = T_L + 1; - auto log_T_L = log(T_grid[T_L]); - auto log_T_R = log(T_grid[T_R]); - auto log_T_target = log(target_T); - // fraction of T_2 to use - auto T_frac = (log_T_target - log_T_L) / (log_T_R - log_T_L); - - auto rho_L = std::distance(rho_grid.begin(), std::lower_bound(rho_grid.begin(), rho_grid.end(), target_rho)); - rho_L = (target_rho < rho_grid[rho_L]) ? rho_L-1 : rho_L; - rho_L = (rho_L == (rho_grid.size()-1)) ? rho_L-1: rho_L; - auto rho_R = rho_L + 1; - auto log_rho_L = log(rho_grid[rho_L]); - auto log_rho_R = log(rho_grid[rho_R]); - auto log_rho_target = log(target_rho); - // fraction of rho_2 to use - auto rho_frac = (log_rho_target - log_rho_L) / (log_rho_R - log_rho_L); - - auto group = std::distance(group_bounds.begin(), std::lower_bound(group_bounds.begin(), group_bounds.end(), target_hnu)); - group = (target_hnu < group_bounds[group]) ? group-1 : group; - group = (group == (group_bounds.size()-1)) ? group-1 : group; - - //std::cout<<"Targets--T: "< build_opacity_spiner_databox(const std::vector &temperature_points, const std::vector &density_points, const std::vector &group_bounds, const std::vector &ramg) { - - Spiner::RegularGrid1D new_temperature(temperature_points.front(), temperature_points.back(), temperature_points.size()); - Spiner::RegularGrid1D new_density(density_points.front(), density_points.back(), density_points.size()); - Spiner::RegularGrid1D new_group_bounds(group_bounds.front(), group_bounds.back(), 2*group_bounds.size()); - Spiner::DataBox opacity_databox(temperature_points.size(), density_points.size(), 2*group_bounds.size()); - opacity_databox.setRange(0, new_temperature); - opacity_databox.setRange(1, new_density); - opacity_databox.setRange(2, new_group_bounds); - - std::cout<<"T rho hnu opac(sq_cm/g)"< keys; - for(int ikey = 0;ikey<32;++ikey) { - keys.push_back(new char(16)); + for (size_t ikey =0;ikey(keys.size()); int nkeys = -1; c_gkeys(filename_char, &mat_ID, keys.data(), &kkeys, &nkeys, &ier); - std::cout<<"nkeys: "< total_keys) { + std::cout<<"POTENTIAL ERROR: nkeys: "< nice_keys; for (int j=0;j mg_opac_keywords{ramg_keyword, rsmg_keyword, pmg_keyword}; + std::array mg_fields{ SP5::Fields::ramg, SP5::Fields::rsmg, SP5::Fields::pmg}; - for (auto key : mg_opac_keywords) { + // TODO Thread these variables through the input + bool log_T_rho_hnu = true; // form a new grid from max and min evenly spaced in log10 + bool log_interpolation = true; // interpolate logarithmically - std::cout<<"Interpolating and building databox for "< opacity_data(nmg, 0.0); c_ggetmg(filename_char, &mat_ID, key.data(), temperature_points.data(), &nt, &post_nt, @@ -280,9 +227,29 @@ int main(int argc, char *argv[]) { group_bounds.data(), &nhnu, &post_nhnu, opacity_data.data(), &nmg, &post_nmg, &ier); - auto spiner_opacity_databox = build_opacity_spiner_databox(temperature_points, density_points, group_bounds, opacity_data); - // save here with keyword? + auto [spiner_opacity_databox, new_group_bounds] = build_opacity_spiner_databox( + temperature_points, density_points, group_bounds, opacity_data, log_T_rho_hnu, + log_interpolation); + + // only save the group bounds once for this material + if (i==0) { + const hsize_t dimensions[] = {static_cast(new_group_bounds.size())}; + std::string dataset_name("group bounds"); + hid_t dataspace_id = H5Screate_simple(1, dimensions, nullptr); + hid_t dataset_id = H5Dcreate2(matGroup, dataset_name.c_str(), H5T_IEEE_F64LE, dataspace_id, + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + const herr_t status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, + new_group_bounds.data()); + } + std::cout<<"Saving databox for "< #include -const std::string DEFAULT_SAVENAME = "opacity.sp5"; +const std::string DEFAULT_SAVENAME = "converted_ipcress.h5";; const std::string EXAMPLESTRING = R"( # air.dat # These are comments. @@ -68,4 +67,4 @@ void parseCLI(int argc, char *argv[], std::string &savename, std::string &filename, bool &printMetadata, std::string &helpMessage); -#endif // _SESAME2SPINER_PARSER_HPP_ +#endif // _IPCRESS2SPINER_PARSER_HPP_ From ce63b036826b5300a30f4dff2a1dc4da133870a5 Mon Sep 17 00:00:00 2001 From: "Alex R. Long" Date: Mon, 3 Aug 2026 08:57:08 -0600 Subject: [PATCH 4/5] + Add supporting header files for ipcress2spiner + Remove commented out code and unused variables from sesame version --- cmake/Findgandolf.cmake | 47 ++++++++++++++ ipcress2spiner/CMakeLists.txt | 2 +- ipcress2spiner/generate_files.hpp | 48 ++++++++++++++ ipcress2spiner/interpolate_opacity.hpp | 90 ++++++++++++++++++++++++++ ipcress2spiner/make_spiner_databox.hpp | 73 +++++++++++++++++++++ ipcress2spiner/parse_cli.cpp | 10 --- ipcress2spiner/parse_cli.hpp | 41 ------------ 7 files changed, 259 insertions(+), 52 deletions(-) create mode 100644 cmake/Findgandolf.cmake create mode 100644 ipcress2spiner/generate_files.hpp create mode 100644 ipcress2spiner/interpolate_opacity.hpp create mode 100644 ipcress2spiner/make_spiner_databox.hpp diff --git a/cmake/Findgandolf.cmake b/cmake/Findgandolf.cmake new file mode 100644 index 0000000..efd46b5 --- /dev/null +++ b/cmake/Findgandolf.cmake @@ -0,0 +1,47 @@ +#------------------------------------------------------------------------------# +# © 2021-2023. 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. +#------------------------------------------------------------------------------# + +# Find the native EOSPAC headers and libraries. +# +# EOSPAC_INCLUDE_DIRS - where to find eos_Interface.h, etc. +# EOSPAC_LIBRARIES - List of libraries when using eospac6. +# EOSPAC_FOUND - True if eospac found. + +#TODO: Add EOSPAC_MODULES (possibly part of include dirs) and clarify which interface we need + +# Look for the header file. +FIND_PATH(GANDOLF_INC_DIR NAMES gandolf.h) + +# Look for the library. +FIND_LIBRARY(GANDOLF_LIB_DIR NAMES libgandolf.a) +#message("Gandolf lib dir; ${GANDOLF_LIB_DIR}") + +# handle the QUIETLY and REQUIRED arguments and set EOSPAC_FOUND to TRUE if +# all listed variables are TRUE +INCLUDE(FindPackageHandleStandardArgs) +FIND_PACKAGE_HANDLE_STANDARD_ARGS(gandolf DEFAULT_MSG GANDOLF_LIB_DIR GANDOLF_INC_DIR) + +# Copy the results to the output variables. +SET(GANDOLF_LIBRARY ${GANDOLF_LIB_DIR}) +SET(GANDOLF_INCLUDE_DIRECTORY ${GANDOLF_INC_DIR}) + +MARK_AS_ADVANCED(GANDOLF_INCLUDE_DIRECTORY GANDOLF_LIBRARY) + +if(GANDOLF_INCLUDE_DIRECTORY AND GANDOLF_LIBRARY) + add_library(gandolf::gandolf STATIC IMPORTED) + set_target_properties(gandolf::gandolf PROPERTIES + IMPORTED_LOCATION "${GANDOLF_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${GANDOLF_INCLUDE_DIRECTORY}") +endif() + diff --git a/ipcress2spiner/CMakeLists.txt b/ipcress2spiner/CMakeLists.txt index bd015be..1163614 100644 --- a/ipcress2spiner/CMakeLists.txt +++ b/ipcress2spiner/CMakeLists.txt @@ -1,5 +1,5 @@ #------------------------------------------------------------------------------ -# © 2021-2023. Triad National Security, LLC. All rights reserved. This +# © 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 diff --git a/ipcress2spiner/generate_files.hpp b/ipcress2spiner/generate_files.hpp new file mode 100644 index 0000000..db0f826 --- /dev/null +++ b/ipcress2spiner/generate_files.hpp @@ -0,0 +1,48 @@ +//====================================================================== +// generate_files function that write spiner file +// Author: Alex R. Long (along@lanl.gov) +// © 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. +//====================================================================== + +#ifndef _IPCRESS2SPINER_GENERATE_FILES_HPP_ +#define _IPCRESS2SPINER_GENERATE_FILES_HPP_ + +#include +#include +#include +#include +#include +#include + +herr_t saveMaterial(hid_t loc, hid_t matGroup, const int matid, + const std::string &sMatid, + const std::string &sp5_field_name, + const std::vector &group_bounds, + Spiner::DataBox &opacity) { + + + double zero_offset =0.0; + herr_t status = 0; + // Dependent variables metadata + status += H5LTset_attribute_string(loc, sMatid.c_str(), SP5::Offsets::messageName, + SP5::Offsets::message); + status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Offsets::rho, + &zero_offset, 1); + status += H5LTset_attribute_double(loc, sMatid.c_str(), SP5::Offsets::T, &zero_offset, 1); + + status += opacity.saveHDF(matGroup, sp5_field_name); + + return status; +} + +#endif // _IPCRESS2SPINER_GENERATE_FILES_HPP_ diff --git a/ipcress2spiner/interpolate_opacity.hpp b/ipcress2spiner/interpolate_opacity.hpp new file mode 100644 index 0000000..4f0e9b2 --- /dev/null +++ b/ipcress2spiner/interpolate_opacity.hpp @@ -0,0 +1,90 @@ +//====================================================================== +// ipcress2spiner tool for converting ipcress to spiner +// Author: Alex R. Long (along@lanl.gov) +// © 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. +//====================================================================== + +#ifndef _IPCRESS2SPINER_INTERPOLATE_HPP_ +#define _IPCRESS2SPINER_INTERPOLATE_HPP_ + +inline double interpolate_mg_opacity_data(const std::vector &T_grid, const std::vector &rho_grid, const std::vector &group_bounds, const std::vector &op_data, + double target_T, double target_rho, double target_hnu, const double log_interpolation) { + + size_t n_T = T_grid.size(); + size_t n_rho = rho_grid.size(); + size_t n_groups = group_bounds.size()-1; + + auto T_L = std::distance(T_grid.begin(), std::lower_bound(T_grid.begin(), T_grid.end(), target_T)); + T_L = (target_T < T_grid[T_L]) ? T_L-1 : T_L; + T_L = (T_L == (T_grid.size()-1)) ? T_L-1: T_L; + auto T_R = T_L + 1; + // fraction of T_2 to use + double T_frac = -1.0; // defaulted to invalid value + if (log_interpolation) { + auto log_T_L = log(T_grid[T_L]); + auto log_T_R = log(T_grid[T_R]); + auto log_T_target = log(target_T); + T_frac = (log_T_target - log_T_L) / (log_T_R - log_T_L); + } + else { + T_frac = (target_T - T_grid[T_L]) / (T_grid[T_R] - T_grid[T_L]); + } + + auto rho_L = std::distance(rho_grid.begin(), std::lower_bound(rho_grid.begin(), rho_grid.end(), target_rho)); + rho_L = (target_rho < rho_grid[rho_L]) ? rho_L-1 : rho_L; + rho_L = (rho_L == (rho_grid.size()-1)) ? rho_L-1: rho_L; + auto rho_R = rho_L + 1; + // fraction of rho_2 to use + double rho_frac = -1.0; // defaulted to invalid value + if (log_interpolation) { + auto log_rho_L = log(rho_grid[rho_L]); + auto log_rho_R = log(rho_grid[rho_R]); + auto log_rho_target = log(target_rho); + rho_frac = (log_rho_target - log_rho_L) / (log_rho_R - log_rho_L); + } + else { + rho_frac = (target_rho - rho_grid[rho_L]) / (rho_grid[rho_R] - rho_grid[rho_L]); + } + + auto group = std::distance(group_bounds.begin(), std::lower_bound(group_bounds.begin(), group_bounds.end(), target_hnu)); + group = (target_hnu < group_bounds[group]) ? group-1 : group; + group = (group == (group_bounds.size()-1)) ? group-1 : group; + + // get the adjacent rows of the opacity index, looks like op_[T]_[rho]_[hnu] with L indicating + // left or index n and "R" indicating right index (n+1) + auto op_L_L = op_data[ T_L * (n_rho*n_groups) + rho_L*(n_groups) + group]; + auto op_L_R = op_data[ T_L * (n_rho*n_groups) + rho_R*(n_groups) + group]; + auto op_R_L = op_data[ T_R * (n_rho*n_groups) + rho_L*(n_groups) + group]; + auto op_R_R = op_data[ T_R * (n_rho*n_groups) + rho_R*(n_groups) + group]; + + // reduce in temperature dimension + auto op_L =(1.0-T_frac) * op_L_L + T_frac*op_R_L; + auto op_R =(1.0-T_frac) * op_L_R + T_frac*op_R_R; + + // reduce in density dimension and return + double interp_opac = (1.0-rho_frac) * op_L + rho_frac*op_R; + if (false) { + std::cout<<"Targets--T: "<, std::vector> build_opacity_spiner_databox(const std::vector &temperature_points, const std::vector &density_points, const std::vector &group_bounds, const std::vector &ramg, const bool log_T_rho_hnu, const bool log_interpolation, const bool verbose_mode = false) { + + const auto n_groups = group_bounds.size() -1; + const double temperature_start = (log_T_rho_hnu) ? log10(temperature_points.front()) : temperature_points.front(); + const double temperature_end = (log_T_rho_hnu) ? log10(temperature_points.back()) : temperature_points.back(); + const double density_start = (log_T_rho_hnu) ? log10(density_points.front()) : density_points.front(); + const double density_end = (log_T_rho_hnu) ? log10(density_points.back()) : density_points.back(); + const double hnu_start = (log_T_rho_hnu) ? log10(group_bounds.front()) : group_bounds.front(); + const double hnu_end = (log_T_rho_hnu) ? log10(group_bounds.back()) : group_bounds.back(); + const double hnu_midpoints_start = (log_T_rho_hnu) ? log10( 0.5*(group_bounds[0]+group_bounds[1])) : 0.5*(group_bounds[0] + group_bounds[1]); + const double hnu_midpoints_end = (log_T_rho_hnu) ? log10(0.5*(group_bounds[n_groups]+group_bounds[n_groups-1])) : 0.5*(group_bounds[n_groups]+group_bounds[n_groups-1]); + + Spiner::RegularGrid1D new_temperature(temperature_start, temperature_end, temperature_points.size()); + Spiner::RegularGrid1D new_density(density_start, density_end, density_points.size()); + Spiner::RegularGrid1D new_group_bounds(hnu_start, hnu_end, group_bounds.size()); + Spiner::RegularGrid1D new_group_centers(hnu_midpoints_start, hnu_midpoints_end, group_bounds.size()-1); + Spiner::DataBox opacity_databox(temperature_points.size(), density_points.size(), n_groups); + opacity_databox.setRange(0, new_temperature); + opacity_databox.setRange(1, new_density); + opacity_databox.setRange(2, new_group_centers); + + if(verbose_mode) { + std::cout<<"T rho hnu opac(sq_cm/g)"< vector_new_group_bounds; + for (size_t i=0;i #include -//#include "io_eospac.hpp" #include "parse_cli.hpp" void parseCLI(int argc, char *argv[], std::string &savename, @@ -43,9 +42,6 @@ void parseCLI(int argc, char *argv[], std::string &savename, << "\t-d: same as -vv\n" << "\t-h: print this message\n" << "\n" - // ARL: Don't have input filenames - //<< "Several example input files:\n" - //<< EXAMPLESTRING << "\n" << std::endl; helpMessage = helpStream.str(); @@ -65,12 +61,6 @@ void parseCLI(int argc, char *argv[], std::string &savename, std::exit(0); } else if (std::strcmp(argv[i], "-p") == 0) { printMetadata = true; - //} else if ((std::strcmp(argv[i], "-w") == 0 || std::strcmp(argv[i], "-v") == 0) && - // eospacWarn == Verbosity::Quiet) { - // eospacWarn = Verbosity::Verbose; - //} else if ((std::strcmp(argv[i], "-d") == 0 || std::strcmp(argv[i], "-vv") == 0) && - // eospacWarn != Verbosity::Debug) { - // eospacWarn = Verbosity::Debug; } else if (std::strcmp(argv[i], "-s") == 0) { savename = argv[++i]; } else { diff --git a/ipcress2spiner/parse_cli.hpp b/ipcress2spiner/parse_cli.hpp index b0f0f4d..174372a 100644 --- a/ipcress2spiner/parse_cli.hpp +++ b/ipcress2spiner/parse_cli.hpp @@ -21,47 +21,6 @@ #include const std::string DEFAULT_SAVENAME = "converted_ipcress.h5";; -const std::string EXAMPLESTRING = R"( -# air.dat -# These are comments. -# The "#" character must be at the beginning of a line. -# only matid is required. All others override defaults. -matid = 5030 -name = air -# rho is in g/cm^3 -rhomin = 1e-2 -rhomax = 10 -numrho = 64 -# T is in Kelvin -Tmin = 252 -Tmax = 1e4 -numT = 32 -# sie is in erg/g -siemin = 1e12 -siemax = 1e16 -numsie = 32 - - -# titanium.dat -matid = 2961 -name = titanium -# These set the number of grid poitns per decade -# for each variable. The default is 50 points. -numrho/decade = 30 -numT/decade = 25 -numSie/decade = 15 - - -# steel.dat -matid=4272 -rhomin = 1e-2 -Tmin = 1 -# These shrink lograithm of bounds -# by a fraction of the total interval <= 1 -shrinklRhoBounds = 0.15 -shrinklTBounds = 0.15 -shrinkleBounds = 0.5 -)"; void parseCLI(int argc, char *argv[], std::string &savename, std::string &filename, bool &printMetadata, From 2a5ff4bebd45d39638977cff097c9b9aeeeba4c8 Mon Sep 17 00:00:00 2001 From: "Alex R. Long" Date: Tue, 4 Aug 2026 14:32:43 -0600 Subject: [PATCH 5/5] Address comments from Wollaeger + Use hints to find Gandolf in Findgandolf.cmake + Guard against Gandolf not found cmake errors + Update copyright dates --- cmake/Findgandolf.cmake | 36 ++++++++++++++++++-------- ipcress2spiner/CMakeLists.txt | 13 +++++----- ipcress2spiner/main.cpp | 1 + ipcress2spiner/make_spiner_databox.hpp | 4 +-- ipcress2spiner/parse_cli.cpp | 2 +- 5 files changed, 35 insertions(+), 21 deletions(-) diff --git a/cmake/Findgandolf.cmake b/cmake/Findgandolf.cmake index efd46b5..0351ced 100644 --- a/cmake/Findgandolf.cmake +++ b/cmake/Findgandolf.cmake @@ -1,5 +1,5 @@ #------------------------------------------------------------------------------# -# © 2021-2023. Triad National Security, LLC. All rights reserved. This +# © 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 @@ -12,22 +12,37 @@ # publicly and display publicly, and to permit others to do so. #------------------------------------------------------------------------------# -# Find the native EOSPAC headers and libraries. +# ARL NOTE: Based off FindEOSPAC.cmake from singularity-eos + +# Find the native GANDOLF headers and libraries. +# +# GANDOLF_INCLUDE_DIRECTORY - Where to find gandolf.h, etc. +# GANDOLF_LIBRARY - Gandolf library +# GANDOLF_FOUND - True if gandolf found. # -# EOSPAC_INCLUDE_DIRS - where to find eos_Interface.h, etc. -# EOSPAC_LIBRARIES - List of libraries when using eospac6. -# EOSPAC_FOUND - True if eospac found. +# Current modules set GANDOLF_INC_DIR and GANDOLF_LIB_DIR and prepend the top +# level of the Gandolf build to CMAKE_PREFIX_PATH + + +# if environment variables are set, use them as hints to FIND calls +set(GANDOLF_INC_DIR_HINTS "") +if(DEFINED ENV{GANDOLF_INC_DIR}) + list(APPEND GANDOLF_INC_DIR_HINTS "$ENV{GANDOLF_INC_DIR}") +endif() + +set(GANDOLF_LIB_DIR_HINTS "") +if(DEFINED ENV{GANDOLF_LIB_DIR}) + list(APPEND GANDOLF_LIB_DIR_HINTS "$ENV{GANDOLF_LIB_DIR}") +endif() -#TODO: Add EOSPAC_MODULES (possibly part of include dirs) and clarify which interface we need # Look for the header file. -FIND_PATH(GANDOLF_INC_DIR NAMES gandolf.h) +FIND_PATH(GANDOLF_INC_DIR NAMES gandolf.h HINTS ${GANDOLF_INC_DIR_HINTS}) # Look for the library. -FIND_LIBRARY(GANDOLF_LIB_DIR NAMES libgandolf.a) -#message("Gandolf lib dir; ${GANDOLF_LIB_DIR}") +FIND_LIBRARY(GANDOLF_LIB_DIR NAMES libgandolf.a HINTS ${GANDOLF_LIB_DIR_HINTS}) -# handle the QUIETLY and REQUIRED arguments and set EOSPAC_FOUND to TRUE if +# handle the QUIETLY and REQUIRED arguments and set GANDOLF_FOUND to TRUE if # all listed variables are TRUE INCLUDE(FindPackageHandleStandardArgs) FIND_PACKAGE_HANDLE_STANDARD_ARGS(gandolf DEFAULT_MSG GANDOLF_LIB_DIR GANDOLF_INC_DIR) @@ -44,4 +59,3 @@ if(GANDOLF_INCLUDE_DIRECTORY AND GANDOLF_LIBRARY) IMPORTED_LOCATION "${GANDOLF_LIBRARY}" INTERFACE_INCLUDE_DIRECTORIES "${GANDOLF_INCLUDE_DIRECTORY}") endif() - diff --git a/ipcress2spiner/CMakeLists.txt b/ipcress2spiner/CMakeLists.txt index 1163614..fc4ca8a 100644 --- a/ipcress2spiner/CMakeLists.txt +++ b/ipcress2spiner/CMakeLists.txt @@ -13,6 +13,9 @@ #------------------------------------------------------------------------------ project(ipcress2spiner VERSION 0.8.0.0 LANGUAGES CXX Fortran) +find_package(gandolf) + +if(GANDOLF_FOUND) add_executable(ipcress2spiner parse_cli.cpp parse_cli.hpp @@ -21,10 +24,6 @@ add_executable(ipcress2spiner main.cpp ) -# Turn on verbosity right before your package search -#set(CMAKE_FIND_DEBUG_MODE ON) -find_package(gandolf) -#set(CMAKE_FIND_DEBUG_MODE OFF) target_include_directories(ipcress2spiner PUBLIC @@ -47,6 +46,6 @@ target_compile_definitions(ipcress2spiner PRIVATE IPCRESS2SPINER_VERSION=\"${PROJECT_VERSION}\" ) install(TARGETS ipcress2spiner DESTINATION ${CMAKE_INSTALL_BINDIR}) - - -# TODO: Add tests for sesame2spiner here. +else() +message("GANDOLF was not found, not building ipcress2spiner") +endif() diff --git a/ipcress2spiner/main.cpp b/ipcress2spiner/main.cpp index 4d0844e..f691a45 100644 --- a/ipcress2spiner/main.cpp +++ b/ipcress2spiner/main.cpp @@ -44,6 +44,7 @@ // ARL: My to-do list // \TODO Thread log interpolation points and log interpolation mode options through command line +// \TODO Thread verbosity option through command-line options to print more info // \TODO Add gray opacties // \TODO Add ionization state? // \TODO clang-format diff --git a/ipcress2spiner/make_spiner_databox.hpp b/ipcress2spiner/make_spiner_databox.hpp index 4984f73..a94e2ab 100644 --- a/ipcress2spiner/make_spiner_databox.hpp +++ b/ipcress2spiner/make_spiner_databox.hpp @@ -61,9 +61,9 @@ inline std::pair, std::vector> build_opacity_sp } // return new group bounds to write separately to hdf5 file - std::vector vector_new_group_bounds; + std::vector vector_new_group_bounds(group_bounds.size()); for (size_t i=0;i