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/cmake/Findgandolf.cmake b/cmake/Findgandolf.cmake new file mode 100644 index 0000000..0351ced --- /dev/null +++ b/cmake/Findgandolf.cmake @@ -0,0 +1,61 @@ +#------------------------------------------------------------------------------# +# © 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. +#------------------------------------------------------------------------------# + +# 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. +# +# 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() + + +# Look for the header file. +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 HINTS ${GANDOLF_LIB_DIR_HINTS}) + +# 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) + +# 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 new file mode 100644 index 0000000..fc4ca8a --- /dev/null +++ b/ipcress2spiner/CMakeLists.txt @@ -0,0 +1,51 @@ +#------------------------------------------------------------------------------ +# © 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. +#------------------------------------------------------------------------------ +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 + make_spiner_databox.hpp + interpolate_opacity.hpp + main.cpp +) + + +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) +# Define the full version as a string macro +target_compile_definitions(ipcress2spiner PRIVATE + IPCRESS2SPINER_VERSION=\"${PROJECT_VERSION}\" +) +install(TARGETS ipcress2spiner DESTINATION ${CMAKE_INSTALL_BINDIR}) +else() +message("GANDOLF was not found, not building ipcress2spiner") +endif() 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: "< +#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 "interpolate_opacity.hpp" +#include "parse_cli.hpp" +#include "make_spiner_databox.hpp" + +// 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 +// \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 "C" void c_gsizeof( char*, int*, char*, int*, int*); +extern "C" void c_gchgrids( char*, int*, int*, int*, int*, int*, 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*, + 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 +// 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*); +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< 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 (size_t ikey =0;ikey(keys.size()); + int nkeys = -1; + c_gkeys(filename_char, &mat_ID, keys.data(), &kkeys, &nkeys, &ier); + if(nkeys > total_keys) { + std::cout<<"POTENTIAL ERROR: 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}; + std::array mg_fields{ SP5::Fields::ramg, SP5::Fields::rsmg, SP5::Fields::pmg}; + + // 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 + + for (size_t i=0;i 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, 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 "<, 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(group_bounds.size()); + for (size_t i=0;i +#include +#include +#include +#include +#include +#include + +#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" + << 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], "-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..174372a --- /dev/null +++ b/ipcress2spiner/parse_cli.hpp @@ -0,0 +1,29 @@ +//====================================================================== +// ipcress2spiner tool for converting eospac 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_PARSER_HPP_ +#define _IPCRESS2SPINER_PARSER_HPP_ + +#include +#include + +const std::string DEFAULT_SAVENAME = "converted_ipcress.h5";; + +void parseCLI(int argc, char *argv[], std::string &savename, + std::string &filename, bool &printMetadata, + std::string &helpMessage); + +#endif // _IPCRESS2SPINER_PARSER_HPP_