diff --git a/CMakeLists.txt b/CMakeLists.txt index f65349da751..aaecf898839 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,5 +1,5 @@ # ------------------------------------------------------------------------------# -# © 2021-2025. Triad National Security, LLC. All rights reserved. This program +# © 2021-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. @@ -11,6 +11,8 @@ # publicly and display publicly, and to permit others to do so. # ------------------------------------------------------------------------------# +# This file was made in part with generative AI + cmake_minimum_required(VERSION 3.19) # Disable "in-source" builds @@ -174,6 +176,7 @@ endif() # singularity-eos Library # ------------------------------------------------------------------------------# +# Load CMake modules for finding dependencies include(singularity-eos/Eigen3) include(singularity-eos/eospac) include(singularity-eos/hdf5) @@ -190,12 +193,31 @@ add_library(singularity-eos::singularity-eos_Common ALIAS singularity-eos_Common add_library(singularity-eos_Interface INTERFACE) add_library(singularity-eos::singularity-eos_Interface ALIAS singularity-eos_Interface) -target_link_libraries(singularity-eos_Interface INTERFACE singularity-eos_Common) +target_link_libraries(singularity-eos_Interface INTERFACE singularity-eos_Common singularity-utils::singularity-utils) target_link_libraries(singularity-eos INTERFACE singularity-eos_Interface) # ld has problems with Clang's libomp.so, use ldd instead target_link_options(singularity-eos_Interface INTERFACE $<$:-fuse-ld=lld>) +# ------------------------------------------------------------------------------# +# Determine if we're in submodule mode +# ------------------------------------------------------------------------------# +# This must happen before finding dependencies + +# checks if this is our build, or we've been imported via `add_subdirectory` +if(NOT CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME) + message( + STATUS + "Detected that `singularity-eos` is a subproject, will configure build in submodule mode" + ) + set(SINGULARITY_SUBMODULE_MODE ON) +endif() + +if(SINGULARITY_FORCE_SUBMODULE_MODE) + message(STATUS "Building as though project was a submodule.") + set(SINGULARITY_SUBMODULE_MODE ON) +endif() + # ------------------------------------------------------------------------------# # Compiler & language setup # ------------------------------------------------------------------------------# @@ -225,31 +247,6 @@ if(SINGULARITY_BUILD_PYTHON) set(CMAKE_POSITION_INDEPENDENT_CODE ON) endif() -# require at least C++17, but allow newer versions to become a client requirement if -# explicitly set at build time (needed for building with newer Kokkos) -if(CMAKE_CXX_STANDARD) - target_compile_features(singularity-eos_Interface INTERFACE cxx_std_${CMAKE_CXX_STANDARD}) -else() - target_compile_features(singularity-eos_Interface INTERFACE cxx_std_17) -endif() - -# checks if this is our build, or we've been imported via `add_subdirectory` NB: -# this should make the `option(SINGULARITY_SUBMODULE_MODE ...)` unnecessary -if(CMAKE_PROJECT_NAME STREQUAL PROJECT_NAME) - set(CMAKE_CXX_EXTENSIONS OFF) -else() - message( - STATUS - "Detected that `singularity-eos` is a subproject, will configure build in submodule mode" - ) - set(SINGULARITY_SUBMODULE_MODE ON) -endif() - -if(SINGULARITY_FORCE_SUBMODULE_MODE) - message(STATUS "Building as though project was a submodule.") - set(SINGULARITY_SUBMODULE_MODE ON) -endif() - # Use to determine if Eigen is used or not set(SINGULARITY_USE_EIGEN OFF @@ -266,6 +263,83 @@ if(SINGULARITY_BUILD_CLOSURE AND NOT SINGULARITY_USE_KOKKOSKERNELS) CACHE BOOL "" FORCE) endif() +# ------------------------------------------------------------------------------# +# Import/Find External Dependencies +# ------------------------------------------------------------------------------# +# Must happen before building singularity-utils (which needs ports-of-call and spiner) +# Dependency order: Kokkos -> ports-of-call -> spiner -> singularity-utils + +if(SINGULARITY_SUBMODULE_MODE) + # add all submodules + message(STATUS "singularity-eos configuring in submodule mode.") + + # Kokkos first (ports-of-call may depend on it if PORTABILITY_STRATEGY_KOKKOS is set) + if(SINGULARITY_USE_KOKKOS) + singularity_import_kokkos() + if(SINGULARITY_USE_KOKKOSKERNELS) + singularity_import_kokkoskernels() + endif() + endif() + + # Then ports-of-call + singularity_import_ports_of_call() + + # Then spiner (depends on ports-of-call) + if(SINGULARITY_USE_SPINER) + singularity_import_spiner() + endif() + + # Eigen (if needed) + if(SINGULARITY_USE_EIGEN) + singularity_import_eigen() + endif() +else() + # use system packages + + # Kokkos first + if(SINGULARITY_USE_KOKKOS) + singularity_find_kokkos() + if(SINGULARITY_USE_KOKKOSKERNELS) + singularity_find_kokkoskernels() + endif() + endif() + + # Then ports-of-call + singularity_find_ports_of_call() + + # Then spiner + if(SINGULARITY_USE_SPINER) + singularity_find_spiner() + endif() + + # Eigen (if needed) + if(SINGULARITY_USE_EIGEN) + singularity_find_eigen() + endif() +endif() + +# Build singularity-utils first - it's a dependency for both +# singularity-eos and sesame2spiner +add_subdirectory(singularity-utils) + +# Enable HDF5 and EOSPAC now that languages and dependencies are set up +if(SINGULARITY_USE_SPINER_WITH_HDF5) + singularity_enable_hdf5(singularity-eos_Common) +endif() + +if(SINGULARITY_USE_EOSPAC) + # NB This will add the `eospac-wrapper` directory. + singularity_enable_eospac(singularity-eos_Common) +endif() + +# require at least C++17, but allow newer versions to become a client requirement if +# explicitly set at build time (needed for building with newer Kokkos) +if(CMAKE_CXX_STANDARD) + target_compile_features(singularity-eos_Interface INTERFACE cxx_std_${CMAKE_CXX_STANDARD}) +else() + target_compile_features(singularity-eos_Interface INTERFACE cxx_std_17) +endif() + # ------------------------------------------------------------------------------# # De-thaw some options # ------------------------------------------------------------------------------# @@ -300,13 +374,28 @@ if(SINGULARITY_BUILD_EXAMPLES) add_subdirectory(example) endif() -# add subdirs +# Enable dependencies for singularity-eos_Interface +singularity_enable_ports_of_call(singularity-eos_Interface) + +if(SINGULARITY_USE_SPINER) + singularity_enable_spiner(singularity-eos_Interface) + # if(SINGULARITY_USE_SPINER_WITH_HDF5) + # singularity_enable_hdf5(singularity-eos_Interface) endif() +endif() + +# ------------------------------------------------------------------------------# +# Add subdirectories +# ------------------------------------------------------------------------------# + if(SINGULARITY_BUILD_PYTHON) add_subdirectory(python) endif() if(SINGULARITY_BUILD_SESAME2SPINER) add_subdirectory(sesame2spiner) + # Link sesame2spiner library into singularity-eos + target_link_libraries(singularity-eos_Interface INTERFACE sesame2spiner-lib) + target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_USE_SESAME2SPINER) endif() if(SINGULARITY_BUILD_STELLARCOLLAPSE2SPINER) @@ -379,59 +468,6 @@ endif() # building are not located outside of the source directory, except for explicit # cases. -if(SINGULARITY_USE_SPINER_WITH_HDF5) - singularity_enable_hdf5(singularity-eos_Common) -endif() - -if(SINGULARITY_USE_EOSPAC) - # NB This will add the `eospac-wrapper` directory. - singularity_enable_eospac(singularity-eos_Common) -endif() - -if(SINGULARITY_SUBMODULE_MODE) - # add all submodules - message(STATUS "singularity-eos configuring in submodule mode.") - singularity_import_ports_of_call() - if(SINGULARITY_USE_SPINER) - singularity_import_spiner() - endif() - if(SINGULARITY_USE_KOKKOS) - singularity_import_kokkos() - if(SINGULARITY_USE_KOKKOSKERNELS) - singularity_import_kokkoskernels() - endif() - endif() - - if(SINGULARITY_USE_EIGEN) - singularity_import_eigen() - endif() -else() - # use system packages - singularity_find_ports_of_call() - if(SINGULARITY_USE_SPINER) - singularity_find_spiner() - endif() - - if(SINGULARITY_USE_KOKKOS) - singularity_find_kokkos() - if(SINGULARITY_USE_KOKKOSKERNELS) - singularity_find_kokkoskernels() - endif() - endif() - if(SINGULARITY_USE_EIGEN) - singularity_find_eigen() - endif() - -endif() - -singularity_enable_ports_of_call(singularity-eos_Interface) - -if(SINGULARITY_USE_SPINER) - singularity_enable_spiner(singularity-eos_Interface) - # if(SINGULARITY_USE_SPINER_WITH_HDF5) - # singularity_enable_hdf5(singularity-eos_Interface) endif() -endif() - # Both the interface (headers) and the library (compiled) need to link to Kokkos # see get_sg_eos.cpp if(SINGULARITY_USE_KOKKOS) @@ -589,8 +625,7 @@ target_compile_options( # `-Wclass-memaccess now default with -Wall but we explicitly # manage this ourselves in our serialization routines. $<$,$>:-Wno-class-memaccess> - # Intel compilers enable -ffast-math even in Debug builds. This disables finite-math-only in - # Debug builds to allow NaN/Inf checks without pages of warnings + # Disable finite-math-only in Debug builds to allow NaN/Inf checks (Intel compilers) $<$,$,$>>:-fno-finite-math-only> ) if (SINGULARITY_STRICT_WARNINGS) diff --git a/cmake/install.cmake b/cmake/install.cmake index dcc19a4948f..668b1eb809d 100644 --- a/cmake/install.cmake +++ b/cmake/install.cmake @@ -65,6 +65,22 @@ install( NAMESPACE "singularity-eos::" DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/singularity-eos) +# Install singularity-utils (required by singularity-eos_Interface) +install( + TARGETS singularity-utils + EXPORT singularity-eos_Interface + DESTINATION ${CMAKE_INSTALL_LIBDIR}) + +# Install sesame2spiner-lib if it was built (required by singularity-eos_Interface when enabled) +if(TARGET sesame2spiner-lib) + install( + TARGETS sesame2spiner-lib + EXPORT singularity-eos_Interface + LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR} + ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR} + INCLUDES DESTINATION include + ) +endif() install( TARGETS singularity-eos_Interface EXPORT singularity-eos_Interface @@ -123,6 +139,12 @@ export( FILE ${CMAKE_CURRENT_BINARY_DIR}/singularity-eos_Common.cmake NAMESPACE singularity-eos::) +# Export singularity-utils with singularity-eos_Interface +export( + TARGETS singularity-utils + FILE ${CMAKE_CURRENT_BINARY_DIR}/singularity-utils.cmake + NAMESPACE singularity-utils::) + export( EXPORT singularity-eos_Interface FILE ${CMAKE_CURRENT_BINARY_DIR}/singularity-eos_Interface.cmake diff --git a/cmake/singularity-eos/hdf5.cmake b/cmake/singularity-eos/hdf5.cmake index 323536f9f99..9d4ca0411ab 100644 --- a/cmake/singularity-eos/hdf5.cmake +++ b/cmake/singularity-eos/hdf5.cmake @@ -40,7 +40,9 @@ macro(singularity_enable_hdf5 target) # If we're doing an in-tree cray build, skip all this stuff. Fixes # in-tree builds on Roci + Chicoma This is if (NOT (IN_TREE AND # CRAY_THING)), but no composite expressions in cmake - if(NOT HDF5_C_COMPILER_EXECUTABLE_NO_INTERROGATE OR NOT SINGULARITY_FORCE_SUBMODULE_MODE) + # Check both HDF5_C_COMPILER_NO_INTERROGATE and HDF5_C_COMPILER_EXECUTABLE_NO_INTERROGATE + # as the variable name may differ across CMake/HDF5 versions + if(NOT ((HDF5_C_COMPILER_NO_INTERROGATE OR HDF5_C_COMPILER_EXECUTABLE_NO_INTERROGATE) AND SINGULARITY_FORCE_SUBMODULE_MODE)) # I'm bailing out here if these arn't filled in. I don't know if this is # correct for every downstream use, but right now we need to enforce some kind diff --git a/eospac-wrapper/CMakeLists.txt b/eospac-wrapper/CMakeLists.txt index f2e0b1bf9cc..8285191c567 100644 --- a/eospac-wrapper/CMakeLists.txt +++ b/eospac-wrapper/CMakeLists.txt @@ -29,7 +29,7 @@ endif() target_include_directories(eospac-wrapper PUBLIC - $ + $ $ ) diff --git a/example/benchmark_spiner_cpu.cpp b/example/benchmark_spiner_cpu.cpp index 35783e84591..195435fd378 100644 --- a/example/benchmark_spiner_cpu.cpp +++ b/example/benchmark_spiner_cpu.cpp @@ -48,7 +48,7 @@ // Spiner headers #include -#include +#include #include #include #include @@ -448,4 +448,4 @@ int main(int argc, char *argv[]) { std::cout << "Benchmark complete for material " << std::to_string(matid) << "\n"; } return 0; -} \ No newline at end of file +} diff --git a/example/benchmark_spiner_kokkos.cpp b/example/benchmark_spiner_kokkos.cpp index 9f9b2cc102f..4f7b49faa0b 100644 --- a/example/benchmark_spiner_kokkos.cpp +++ b/example/benchmark_spiner_kokkos.cpp @@ -37,7 +37,7 @@ // Spiner headers #include -#include +#include #include #include #include diff --git a/example/eos_grid.cpp b/example/eos_grid.cpp index e89ab6a1176..271df371992 100644 --- a/example/eos_grid.cpp +++ b/example/eos_grid.cpp @@ -28,7 +28,7 @@ #include // This contains useful tools for preventing things like divide by zero -#include +#include // Needed to import the eos models #include diff --git a/example/get_sound_speed_press.cpp b/example/get_sound_speed_press.cpp index 7114f1ad0f0..ab793268cf0 100644 --- a/example/get_sound_speed_press.cpp +++ b/example/get_sound_speed_press.cpp @@ -17,7 +17,7 @@ #include // This contains useful tools for preventing things like divide by zero -#include +#include // Needed to import the eos models #include // One way of initializing models with modifiers diff --git a/example/map_pt_space.cpp b/example/map_pt_space.cpp index 7fbeefd715d..79f67a31d45 100644 --- a/example/map_pt_space.cpp +++ b/example/map_pt_space.cpp @@ -41,9 +41,9 @@ #include // This contains useful tools for preventing things like divide by zero -#include +#include // 1D root finding -#include +#include // Needed to import the eos models #include diff --git a/example/plugin/dust/dust.hpp b/example/plugin/dust/dust.hpp index 03c737d4ec2..123711756ea 100644 --- a/example/plugin/dust/dust.hpp +++ b/example/plugin/dust/dust.hpp @@ -21,8 +21,8 @@ // Base stuff #include #include -#include #include +#include namespace singularity { diff --git a/example/plugin/dust/dust_variant.hpp b/example/plugin/dust/dust_variant.hpp index cba779a5374..88d257f2760 100644 --- a/example/plugin/dust/dust_variant.hpp +++ b/example/plugin/dust/dust_variant.hpp @@ -30,7 +30,7 @@ // Base stuff #include #include -#include +#include // EOS models #include diff --git a/example/pte_2mat.cpp b/example/pte_2mat.cpp index d51af00c548..8ae66778710 100644 --- a/example/pte_2mat.cpp +++ b/example/pte_2mat.cpp @@ -62,9 +62,9 @@ #include // This contains logic for indexers -#include +#include // This contains useful tools for preventing things like divide by zero -#include +#include // The PTE closures #include // Needed to import the eos models diff --git a/python/module.hpp b/python/module.hpp index 00309d763ea..834292f13e9 100644 --- a/python/module.hpp +++ b/python/module.hpp @@ -25,7 +25,7 @@ #include #include -#include +#include #include namespace py = pybind11; diff --git a/sesame2spiner/CMakeLists.txt b/sesame2spiner/CMakeLists.txt index 5d622d92be5..29c9ace5154 100644 --- a/sesame2spiner/CMakeLists.txt +++ b/sesame2spiner/CMakeLists.txt @@ -1,5 +1,5 @@ #------------------------------------------------------------------------------ -# © 2021-2023. Triad National Security, LLC. All rights reserved. This +# © 2021-2025. Triad National Security, LLC. All rights reserved. This # program was produced under U.S. Government contract 89233218CNA000001 # for Los Alamos National Laboratory (LANL), which is operated by Triad # National Security, LLC for the U.S. Department of Energy/National @@ -12,30 +12,75 @@ # publicly and display publicly, and to permit others to do so. #------------------------------------------------------------------------------ + +# ------------------------------------------------------------------------------ +# sesame2spiner Library +# ------------------------------------------------------------------------------ + +add_library(sesame2spiner-lib + sesame2spiner/io_eospac.cpp + sesame2spiner/io_eospac.hpp + sesame2spiner/generate_files.cpp + sesame2spiner/generate_files.hpp + sesame2spiner/parser.cpp + sesame2spiner/parser.hpp +) + +# Set library output name (without -lib suffix) +set_target_properties(sesame2spiner-lib PROPERTIES OUTPUT_NAME sesame2spiner) + +target_include_directories(sesame2spiner-lib + PUBLIC + $ + $ + $ +) + +# sesame2spiner library depends on singularity-utils (not full singularity-eos) +target_link_libraries(sesame2spiner-lib + PUBLIC + singularity-utils::singularity-utils + eospac-wrapper + singularity-eos::singularity-eos_Common # For HDF5 +) + +# Version information for library +target_compile_definitions(sesame2spiner-lib + PUBLIC + SESAME2SPINER_VERSION=\"${PROJECT_VERSION}\" +) + +# Create alias for consistency +add_library(sesame2spiner::sesame2spiner ALIAS sesame2spiner-lib) + +# ------------------------------------------------------------------------------ +# sesame2spiner Executable +# ------------------------------------------------------------------------------ + add_executable(sesame2spiner - io_eospac.cpp - io_eospac.hpp - generate_files.cpp - generate_files.hpp - parse_cli.cpp - parse_cli.hpp - parser.cpp - parser.hpp - main.cpp + src/parse_cli.cpp + src/parse_cli.hpp + src/main.cpp ) target_include_directories(sesame2spiner - PUBLIC - $ - $ + PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR}/src ) target_link_libraries(sesame2spiner - PUBLIC - singularity-eos::singularity-eos + PRIVATE + sesame2spiner-lib ) +# ------------------------------------------------------------------------------ +# Installation +# ------------------------------------------------------------------------------ + install(TARGETS sesame2spiner DESTINATION ${CMAKE_INSTALL_BINDIR}) -# TODO: Add tests for sesame2spiner here. +install(DIRECTORY sesame2spiner/ + DESTINATION include/sesame2spiner +) + diff --git a/sesame2spiner/generate_files.cpp b/sesame2spiner/sesame2spiner/generate_files.cpp similarity index 98% rename from sesame2spiner/generate_files.cpp rename to sesame2spiner/sesame2spiner/generate_files.cpp index 23bfa4406e3..191dcfc408b 100644 --- a/sesame2spiner/generate_files.cpp +++ b/sesame2spiner/sesame2spiner/generate_files.cpp @@ -22,24 +22,23 @@ #include #include +#ifdef SPINER_USE_HDF #include #include - -#ifndef SPINER_USE_HDF +#else #error "HDF5 must be enabled" #endif // SPINER_USE_HDF #include #include -#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; @@ -47,6 +46,8 @@ using singularity::spiner_table_builder::constructRhoBounds; using singularity::spiner_table_builder::constructTBounds; using singularity::spiner_table_builder::SpinerTableGridParams; +namespace sesame2spiner { + 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, @@ -187,7 +188,7 @@ herr_t saveAllMaterials(const std::string &savename, file = H5Fcreate(savename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); // singularity version - H5LTset_attribute_string(file, "/", "singularity_version", SINGULARITY_VERSION); + H5LTset_attribute_string(file, "/", "singularity_version", SESAME2SPINER_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); @@ -476,3 +477,4 @@ bool checkValInMatBounds(int matid, const std::string &name, Real val, Real vmin } return true; } +} // namespace sesame2spiner diff --git a/sesame2spiner/generate_files.hpp b/sesame2spiner/sesame2spiner/generate_files.hpp similarity index 98% rename from sesame2spiner/generate_files.hpp rename to sesame2spiner/sesame2spiner/generate_files.hpp index 7fc5e0d890b..545129d67cf 100644 --- a/sesame2spiner/generate_files.hpp +++ b/sesame2spiner/sesame2spiner/generate_files.hpp @@ -32,6 +32,8 @@ using namespace EospacWrapper; using singularity::spiner_table_builder::SpinerTableGridParams; +namespace sesame2spiner { + constexpr int PPD_DEFAULT_RHO = 350; constexpr int PPD_DEFAULT_T = 100; constexpr Real STRICTLY_POS_MIN_RHO = 1e-8; @@ -76,5 +78,6 @@ bool checkValInMatBounds(int matid, const std::string &name, Real val, Real vmin Real vmax); int getNumPointsFromPPD(Real min, Real max, int ppd); +} // namespace sesame2spiner #endif // _SESAME2SPINER_GENERATE_FILES_HPP_ diff --git a/sesame2spiner/io_eospac.cpp b/sesame2spiner/sesame2spiner/io_eospac.cpp similarity index 99% rename from sesame2spiner/io_eospac.cpp rename to sesame2spiner/sesame2spiner/io_eospac.cpp index 23290727f5b..7fcd15abef9 100644 --- a/sesame2spiner/io_eospac.cpp +++ b/sesame2spiner/sesame2spiner/io_eospac.cpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2025. Triad National Security, LLC. All rights reserved. This +// © 2021-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 @@ -23,6 +23,7 @@ #include "io_eospac.hpp" +namespace sesame2spiner { // 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, @@ -422,3 +423,4 @@ void modifyNames(TableSplit split, std::vector &names) { } } } // namespace impl +} // namespace sesame2spiner diff --git a/sesame2spiner/io_eospac.hpp b/sesame2spiner/sesame2spiner/io_eospac.hpp similarity index 94% rename from sesame2spiner/io_eospac.hpp rename to sesame2spiner/sesame2spiner/io_eospac.hpp index 2a2c166da2d..e93022d00a0 100644 --- a/sesame2spiner/io_eospac.hpp +++ b/sesame2spiner/sesame2spiner/io_eospac.hpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2025. Triad National Security, LLC. All rights reserved. This +// © 2021-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 @@ -27,13 +27,14 @@ #endif #include -#include -#include -#include +#include +#include +#include #include #include +namespace sesame2spiner { using EospacWrapper::Verbosity; constexpr int NGRIDS = 3; using singularity::TableSplit; @@ -93,5 +94,6 @@ T select(TableSplit split, T a, T b, T c) { } void modifyNames(TableSplit split, std::vector &names); } // namespace impl +} // namespace sesame2spiner #endif // _SESAME2SPINER_IO_EOSPAC_HPP_ diff --git a/sesame2spiner/parser.cpp b/sesame2spiner/sesame2spiner/parser.cpp similarity index 97% rename from sesame2spiner/parser.cpp rename to sesame2spiner/sesame2spiner/parser.cpp index 7e8b0568b24..23e6cd63fe1 100644 --- a/sesame2spiner/parser.cpp +++ b/sesame2spiner/sesame2spiner/parser.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-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 @@ -14,6 +14,7 @@ #include #include +#include #include #include #include @@ -21,9 +22,10 @@ #include #include -#include "parse_cli.hpp" #include "parser.hpp" +namespace sesame2spiner { + Params::Params(const std::string &input_file) { std::ifstream config_file(input_file); if (config_file.is_open()) { @@ -158,4 +160,5 @@ void AddMaterials(std::vector ¶ms, std::vector &matids, } else { throw std::runtime_error("Couldn't open config file " + input_file + "\n"); } -} \ No newline at end of file +} +} // namespace sesame2spiner diff --git a/sesame2spiner/parser.hpp b/sesame2spiner/sesame2spiner/parser.hpp similarity index 69% rename from sesame2spiner/parser.hpp rename to sesame2spiner/sesame2spiner/parser.hpp index 10b9a1e1d2f..9b79ca8b13b 100644 --- a/sesame2spiner/parser.hpp +++ b/sesame2spiner/sesame2spiner/parser.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-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 @@ -21,6 +21,49 @@ #include #include +namespace sesame2spiner { +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 +)"; + // Parse a simple parameter file with // "#" denoting comments. class Params { @@ -49,4 +92,5 @@ class Params { void AddMaterials(std::vector ¶ms, std::vector &matids, const std::string &input_file); +} // namespace sesame2spiner #endif // SESAME2SPINER_PARSER_HPP_ diff --git a/sesame2spiner/main.cpp b/sesame2spiner/src/main.cpp similarity index 85% rename from sesame2spiner/main.cpp rename to sesame2spiner/src/main.cpp index a7308b0edc2..f9eadd5d684 100644 --- a/sesame2spiner/main.cpp +++ b/sesame2spiner/src/main.cpp @@ -1,7 +1,7 @@ //====================================================================== // sesame2spiner tool for converting eospac to spiner // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-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 @@ -30,13 +30,14 @@ #endif #include -#include +#include #include -#include "generate_files.hpp" -#include "io_eospac.hpp" #include "parse_cli.hpp" +#include +#include +using EospacWrapper::Verbosity; int main(int argc, char *argv[]) { std::vector filenames; std::string savename, helpMessage; @@ -52,7 +53,8 @@ int main(int argc, char *argv[]) { << "-----------------------------------------\n" << std::endl; - status = saveAllMaterials(savename, filenames, printMetadata, eospacWarn); + status = + sesame2spiner::saveAllMaterials(savename, filenames, printMetadata, eospacWarn); std::cout << "Done." << std::endl; diff --git a/sesame2spiner/parse_cli.cpp b/sesame2spiner/src/parse_cli.cpp similarity index 96% rename from sesame2spiner/parse_cli.cpp rename to sesame2spiner/src/parse_cli.cpp index bf4f9171e82..a9104bbfc07 100644 --- a/sesame2spiner/parse_cli.cpp +++ b/sesame2spiner/src/parse_cli.cpp @@ -22,8 +22,9 @@ #include #include -#include "io_eospac.hpp" #include "parse_cli.hpp" +#include +#include void parseCLI(int argc, char *argv[], std::string &savename, std::vector &filenames, bool &printMetadata, @@ -46,7 +47,7 @@ void parseCLI(int argc, char *argv[], std::string &savename, << "\t-h: print this message\n" << "\n" << "Several example input files:\n" - << EXAMPLESTRING << "\n" + << sesame2spiner::EXAMPLESTRING << "\n" << std::endl; helpMessage = helpStream.str(); diff --git a/sesame2spiner/parse_cli.hpp b/sesame2spiner/src/parse_cli.hpp similarity index 64% rename from sesame2spiner/parse_cli.hpp rename to sesame2spiner/src/parse_cli.hpp index 7c4c35da60b..582b5759917 100644 --- a/sesame2spiner/parse_cli.hpp +++ b/sesame2spiner/src/parse_cli.hpp @@ -17,53 +17,13 @@ #ifndef _SESAME2SPINER_PARSER_HPP_ #define _SESAME2SPINER_PARSER_HPP_ -#include "io_eospac.hpp" +#include #include #include const std::string DEFAULT_SAVENAME = "materials.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 -)"; +using EospacWrapper::Verbosity; void parseCLI(int argc, char *argv[], std::string &savename, std::vector &filenames, bool &printMetadata, Verbosity &eospacWarn, std::string &helpMessage); diff --git a/sesame2spiner/test.cpp b/sesame2spiner/test/test.cpp similarity index 100% rename from sesame2spiner/test.cpp rename to sesame2spiner/test/test.cpp diff --git a/singularity-eos/base/constants.hpp b/singularity-eos/base/constants.hpp index 2f0af656c55..b7edc9b3d9a 100644 --- a/singularity-eos/base/constants.hpp +++ b/singularity-eos/base/constants.hpp @@ -16,6 +16,7 @@ #define SINGULARITY_EOS_BASE_CONSTANTS_HPP_ #include +#include // TableSplit namespace singularity { @@ -37,7 +38,6 @@ constexpr unsigned long all_values = (1 << 10) - 1; constexpr size_t MAX_NUM_LAMBDAS = 3; enum class DataStatus { Deallocated = 0, OnDevice = 1, OnHost = 2, UnManaged = 3 }; enum class TableStatus { OnTable = 0, OffBottom = 1, OffTop = 2 }; -enum class TableSplit { Total = 0, ElectronOnly = 1, IonCold = 2 }; constexpr Real ROOM_TEMPERATURE = 293; // K constexpr Real ATMOSPHERIC_PRESSURE = 1e6; diff --git a/singularity-eos/base/error_utils.hpp b/singularity-eos/base/error_utils.hpp index 4e0527f60f3..3153ce21560 100644 --- a/singularity-eos/base/error_utils.hpp +++ b/singularity-eos/base/error_utils.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2025. 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 @@ -15,74 +15,6 @@ #ifndef SINGULARITY_EOS_BASE_ERROR_UTILS_HPP_ #define SINGULARITY_EOS_BASE_ERROR_UTILS_HPP_ -#include -#include - -#include - -namespace singularity { -namespace error_utils { - -using PortsOfCall::printf; - -constexpr double _NORMAL_FACTOR = 1.0e10; - -struct is_normal_or_zero { - template - constexpr bool PORTABLE_FORCEINLINE_FUNCTION operator()(valT value) const { - static_assert(std::is_floating_point::value); - return (value == valT{0}) || - (std::isnormal(_NORMAL_FACTOR * value) && std::isnormal(value)); - } -}; - -struct is_strictly_positive { - template - constexpr bool PORTABLE_FORCEINLINE_FUNCTION operator()(valT value) const { - static_assert(std::is_arithmetic::value); - return value > valT{0}; - } -}; - -struct is_non_negative { - template - constexpr bool PORTABLE_FORCEINLINE_FUNCTION operator()(valT value) const { - static_assert(std::is_arithmetic::value); - return value >= valT{0}; - } -}; - -// Checks whether a value obeys some sort of provided condition. If not, returns true and -// prints the provided error message, variable name, and value (but does not abort!) -template -PORTABLE_FORCEINLINE_FUNCTION bool violates_condition(valT &&value, condT &&condition, - nameT &&var_name) { - const bool good = condition(std::forward(value)); - if (!good) { - printf("### ERROR: Bad singularity-eos value\n Var: %s\n Value: %.15e\n", - var_name, value); - } - return !good; -} - -// Create specialized functions using the above conditions -template -PORTABLE_FORCEINLINE_FUNCTION bool bad_value(valT &&value, nameT &&var_name) { - return violates_condition(std::forward(value), is_normal_or_zero{}, - std::forward(var_name)); -} -template -PORTABLE_FORCEINLINE_FUNCTION bool non_positive_value(valT &&value, nameT &&var_name) { - return violates_condition(std::forward(value), is_strictly_positive{}, - std::forward(var_name)); -} -template -PORTABLE_FORCEINLINE_FUNCTION bool negative_value(valT &&value, nameT &&var_name) { - return violates_condition(std::forward(value), is_strictly_positive{}, - std::forward(var_name)); -} - -} // namespace error_utils -} // namespace singularity +#include #endif // #ifndef SINGULARITY_EOS_BASE_ERROR_UTILS_HPP_ diff --git a/singularity-eos/base/fast-math/logs.hpp b/singularity-eos/base/fast-math/logs.hpp index a986d8a1c25..b11e6107103 100644 --- a/singularity-eos/base/fast-math/logs.hpp +++ b/singularity-eos/base/fast-math/logs.hpp @@ -1,5 +1,5 @@ //====================================================================== -// © 2021-2025. 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,286 +12,9 @@ // publicly and display publicly, and to permit others to do so. //====================================================================== -#ifndef _SINGULARITY_EOS_UTILS_FAST_MATH_LOGS_ -#define _SINGULARITY_EOS_UTILS_FAST_MATH_LOGS_ -#include -#include -#include +#ifndef _SINGULARITY_EOS_BASE_FAST_MATH_LOGS_ +#define _SINGULARITY_EOS_BASE_FAST_MATH_LOGS_ -#include -#include -/* - * These functions are for use when moving to and from the gridding space - * for Spiner grids. - * - * The core idea here is to move onto a grid that's roughly - * logarithmically spaced. It doesn't actually matter if that grid is - * EXACTLY logarithmic. Just that it is approximately so. - * There are however some constraints that do matter: - * - * 1. The function we use to translate to "grid space" must be - * precise, meaning you get the same answer each time. - * - * 2. The function must be invertible, meaning you can go back to - * "linear space." Ideally it is also continuous. - * - * 3. The function and its inverse must be fast. - * - * To meet these constraints, we split a real number into its mantissa - * + exponent in base 2. The mantissa is a real number on the interval - * [0.5, 1) and the exponent is an integer such that - * - * mantissa * 2^exponent = number - * - * The log in base 2 is then - * - * lg(number) = lg(mantissa) + exponent - * - * so our goal is to approximate lg(mantissa) satisfying the above - * criteria. If we do, then we will achieve the desired goals. We have - * two approaches. The first is a linear approximation of - * lg(mantissa). The second is a quadratic one. The latter allows for - * continuity of the derivative. - * - * See ArXiv:2206.08957 - */ +#include -// TODO(JMM): All the math here is all assuming inputs are -// doubles. Should add single-precision overloads. - -namespace singularity { -namespace FastMath { - -// TODO(JMM): switch from preprocessor macros to CMake generated C++ -enum LogType { SINGLE = -1, DOUBLE = 0, NQT1 = 1, NQT2 = 2 }; -namespace Settings { -#ifdef SINGULARITY_USE_TRUE_LOG_GRIDDING -constexpr bool TRUE_LOGS = true; -#else -constexpr bool TRUE_LOGS = false; -#endif // SINGULARITY_USE_TRUE_LOG_GRIDDING -constexpr bool NQT_LOGS = !TRUE_LOGS; - -#ifdef SINGULARITY_USE_SINGLE_LOGS -constexpr bool FP32 = true; -#else -constexpr bool FP32 = false; -#endif // SINGULARITY_USE_SINGLE_LOGS -constexpr bool FP64 = !FP32; - -#ifdef SINGULARITY_NQT_PORTABLE -constexpr bool NQT_PORTABLE = true; -#else -constexpr bool NQT_PORTABLE = false; -#endif // SINGULARITY_NQT_PORTABLE - -#ifdef SINGULARITY_NQT_ORDER_1 -constexpr bool NQT_O1 = true; -#else -constexpr bool NQT_O1 = false; -#endif // SINGULARITY_NQT_O1 -constexpr LogType log_type = (TRUE_LOGS && FP32) ? LogType::SINGLE - : (TRUE_LOGS && FP64) ? LogType::DOUBLE - : NQT_O1 ? LogType::NQT1 - : LogType::NQT2; -} // namespace Settings - -template -PORTABLE_FORCEINLINE_FUNCTION auto sgn(const T &x) { - return (T(0) < x) - (x < T(0)); -} - -// TODO(JMM): Add templating and concepts to enforce bit width -// equivalence in input/output types. -PORTABLE_FORCEINLINE_FUNCTION -auto as_int(double f) { return *reinterpret_cast(&f); } -PORTABLE_FORCEINLINE_FUNCTION -auto as_double(std::int64_t i) { return *reinterpret_cast(&i); } - -// Magic numbers constexpr because C++ doesn't constexpr reinterpret casts -// These numbers will be different for different precisions and endianness. -namespace FP64LE { // 64 bit, little endian -constexpr std::int64_t one = 1; -// as_int(1.0) == 2^62 - 2^52 -constexpr std::int64_t one_as_int = (one << 62) - (one << 52); -// 1./static_cast(as_int(2.0) - as_int(1.0)) == 2^-52 -constexpr double scale_down = 2.22044604925031e-16; -// as_int(2.0) - as_int(1.0) = 2^52, but note the type -constexpr double scale_up = (one << 52); -// 2^52 - 1 -constexpr std::int64_t mantissa_mask = (one << 52) - one; -// 2^26 - 1 -constexpr std::int64_t low_mask = (one << 26) - 1; -} // namespace FP64LE - -// First order interpolation based NQTs -// ---------------------------------------------------------------------- -// Reference implementations, however the integer cast implementation -// below is probably faster. -PORTABLE_FORCEINLINE_FUNCTION -double lg_o1_portable(const double x) { - int e; - PORTABLE_REQUIRE(x > 0, "log divergent for x <= 0"); - const double m = frexp(x, &e); - return 2 * (m - 1) + e; -} - -PORTABLE_FORCEINLINE_FUNCTION -double pow2_o1_portable(const double x) { - const int flr = std::floor(x); - const double remainder = x - flr; - const double mantissa = 0.5 * (remainder + 1); - const double exponent = flr + 1; - return ldexp(mantissa, exponent); -} - -// Integer aliased versions -PORTABLE_FORCEINLINE_FUNCTION -double lg_o1_aliased(const double x) { - using namespace FP64LE; - PORTABLE_REQUIRE(x > 0, "log divergent for x <= 0"); - return static_cast(as_int(x) - one_as_int) * scale_down; -} - -PORTABLE_FORCEINLINE_FUNCTION -double pow2_o1_aliased(const double x) { - using namespace FP64LE; - return as_double(static_cast(x * scale_up) + one_as_int); -} -// ---------------------------------------------------------------------- - -// Second-order interpolation based NQTs -// These implementations are due to Peter Hammond -// ---------------------------------------------------------------------- -// Portable versions that use frexp/ldexp rather than integer aliasing -PORTABLE_FORCEINLINE_FUNCTION -double lg_o2_portable(const double x) { - PORTABLE_REQUIRE(x > 0, "log divergent for x <= 0"); - constexpr double four_thirds = 4. / 3.; - int e; - const double m = frexp(x, &e); - return e - four_thirds * (m - 2) * (m - 1); -} - -// This version uses the exact formula -PORTABLE_FORCEINLINE_FUNCTION -double pow2_o2_portable(const double x) { - // log2(mantissa). should go between -1 and 0 - const int flr = std::floor(x); - const double lm = x - flr - 1; - const double mantissa = 0.5 * (3 - std::sqrt(1 - 3 * lm)); - const double exponent = flr + 1; - return ldexp(mantissa, exponent); -} - -// Integer aliased/bithacked versions -PORTABLE_FORCEINLINE_FUNCTION -double lg_o2_aliased(const double x) { - using namespace FP64LE; - PORTABLE_REQUIRE(x > 0, "log divergent for x <= 0"); - const std::int64_t x_as_int = as_int(x) - one_as_int; - const std::int64_t frac_as_int = x_as_int & mantissa_mask; - const std::int64_t frac_high = frac_as_int >> 26; - const std::int64_t frac_low = frac_as_int & low_mask; - const std::int64_t frac_squared = - frac_high * frac_high + ((frac_high * frac_low) >> 25); - - return static_cast(x_as_int + ((frac_as_int - frac_squared) / 3)) * scale_down; -} - -PORTABLE_FORCEINLINE_FUNCTION -double pow2_o2_aliased(const double x) { - using namespace FP64LE; - constexpr std::int64_t a = 9007199254740992; // 2 * 2^52 - constexpr double b = 67108864; // 2^26 - constexpr std::int64_t c = 18014398509481984; // 4 * 2^52 - const std::int64_t x_as_int = static_cast(x * scale_up); - const std::int64_t frac_as_int = x_as_int & mantissa_mask; - const std::int64_t frac_sqrt = - static_cast(b * std::sqrt(static_cast(c - 3 * frac_as_int))); - - return as_double(x_as_int + a - frac_sqrt - frac_as_int + one_as_int); -} -// ---------------------------------------------------------------------- - -PORTABLE_FORCEINLINE_FUNCTION -double fastlg(const double x) { -#ifdef SINGULARITY_NQT_ORDER_1 -#ifdef SINGULARITY_NQT_PORTABLE - return lg_o1_portable(x); -#else - return lg_o1_aliased(x); -#endif // PORTABLE -#else -#ifdef SINGULARITY_NQT_PORTABLE - return lg_o2_portable(x); -#else - return lg_o2_aliased(x); -#endif // PORTABLE -#endif // NQT_ORDER -} - -PORTABLE_FORCEINLINE_FUNCTION -double fastpow2(const double x) { -#ifdef SINGULARITY_NQT_ORDER_1 -#ifdef SINGULARITY_NQT_PORTABLE - return pow2_o1_portable(x); -#else - return pow2_o1_aliased(x); -#endif // PORTABLE -#else -#ifdef SINGULARITY_NQT_PORTABLE - return pow2_o2_portable(x); -#else - return pow2_o2_aliased(x); -#endif // PORTABLE -#endif // NQT_ORDER -} - -PORTABLE_FORCEINLINE_FUNCTION -double lg(const double x) { - PORTABLE_REQUIRE(x > 0, "log divergent for x <= 0"); - if constexpr (Settings::NQT_LOGS) { - // Default expression - return fastlg(x); - } else { - if constexpr (Settings::FP64) { - // double precision is default - return std::log2(x); - } else { - return std::log2f(x); - } - } -} - -PORTABLE_FORCEINLINE_FUNCTION -double pow2(const double x) { - if constexpr (Settings::NQT_LOGS) { - // Default expression - return fastpow2(x); - } else { - if constexpr (Settings::FP64) { - // double precision is default - return std::exp2(x); - } else { - return std::exp2f(x); - } - } -} - -PORTABLE_FORCEINLINE_FUNCTION -double log10(const double x) { - constexpr double LOG2OLOG10 = 0.301029995663981195; - return LOG2OLOG10 * lg(x); -} - -PORTABLE_FORCEINLINE_FUNCTION -double pow10(const double x) { - constexpr double LOG10OLOG2 = 3.321928094887362626; - return pow2(LOG10OLOG2 * x); -} - -} // namespace FastMath -} // namespace singularity - -#endif // _SINGULARITY_EOS_UTILS_FAST_MATH_LOGS_ +#endif // _SINGULARITY_EOS_BASE_FAST_MATH_LOGS_ diff --git a/singularity-eos/base/generic_indexer.hpp b/singularity-eos/base/generic_indexer.hpp index 70c35a2ea79..82dbe3d8ca7 100644 --- a/singularity-eos/base/generic_indexer.hpp +++ b/singularity-eos/base/generic_indexer.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2025. 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 @@ -15,48 +15,6 @@ #ifndef SINGULARITY_EOS_BASE_GENERIC_INDEXER_HPP_ #define SINGULARITY_EOS_BASE_GENERIC_INDEXER_HPP_ -#include - -namespace singularity { - -template -struct GenericIndexer { - Arr arr_; - Map map_; - - template - constexpr GenericIndexer(ArrT_ &&arr_in, MapT_ &&map_in) - : arr_(std::forward(arr_in)), map_(std::forward(map_in)) {} - - // & : non-const lvalue - template - constexpr decltype(auto) operator[](I i) & { - return arr_[map_[i]]; - } - - // const& : const lvalue - template - constexpr decltype(auto) operator[](I i) const & { - return arr_[map_[i]]; - } - - // && : rvalue (indexer is a temporary) — forward arr_’s value category - template - constexpr decltype(auto) operator[](I i) && { - return std::forward(arr_)[map_[i]]; // move only if Arr is a value type - } - - // const rvalue indexer - template - constexpr decltype(auto) operator[](I i) const && { - return std::forward(arr_)[map_[i]]; // preserves const - } -}; - -// CTAD: preserve references for lvalues, decay for rvalues -template -GenericIndexer(ArrT_ &&, MapT_ &&) -> GenericIndexer; - -} // namespace singularity +#include #endif // #ifndef SINGULARITY_EOS_BASE_GENERIC_INDEXER_HPP_ diff --git a/singularity-eos/base/hermite.hpp b/singularity-eos/base/hermite.hpp index ebee68258fa..5727a516daf 100644 --- a/singularity-eos/base/hermite.hpp +++ b/singularity-eos/base/hermite.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 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 @@ -15,85 +15,6 @@ #ifndef SINGULARITY_EOS_BASE_HERMITE_HPP_ #define SINGULARITY_EOS_BASE_HERMITE_HPP_ -#include -#include +#include -namespace singularity { -namespace hermite { - -// psi0 and its derivatives -PORTABLE_INLINE_FUNCTION Real psi0(Real z) { - return z * z * z * (z * (-6.0 * z + 15.0) - 10.0) + 1.0; -} -PORTABLE_INLINE_FUNCTION Real dpsi0(Real z) { - return z * z * (z * (-30.0 * z + 60.0) - 30.0); -} -PORTABLE_INLINE_FUNCTION Real ddpsi0(Real z) { - return z * (z * (-120.0 * z + 180.0) - 60.0); -} - -// psi1 and its derivatives -PORTABLE_INLINE_FUNCTION Real psi1(Real z) { - return z * (z * z * (z * (-3.0 * z + 8.0) - 6.0) + 1.0); -} -PORTABLE_INLINE_FUNCTION Real dpsi1(Real z) { - return z * z * (z * (-15.0 * z + 32.0) - 18.0) + 1.0; -} -PORTABLE_INLINE_FUNCTION Real ddpsi1(Real z) { - return z * (z * (-60.0 * z + 96.0) - 36.0); -} - -// psi2 and its derivatives -PORTABLE_INLINE_FUNCTION Real psi2(Real z) { - return 0.5 * z * z * (z * (z * (-z + 3.0) - 3.0) + 1.0); -} -PORTABLE_INLINE_FUNCTION Real dpsi2(Real z) { - return 0.5 * z * (z * (z * (-5.0 * z + 12.0) - 9.0) + 2.0); -} -PORTABLE_INLINE_FUNCTION Real ddpsi2(Real z) { - return 0.5 * (z * (z * (-20.0 * z + 36.0) - 18.0) + 2.0); -} - -// biquintic hermite polynomial -PORTABLE_INLINE_FUNCTION Real h5(Real fi[36], Real w0t, Real w1t, Real w2t, Real w0mt, - Real w1mt, Real w2mt, Real w0d, Real w1d, Real w2d, - Real w0md, Real w1md, Real w2md) { - return fi[0] * w0d * w0t + fi[1] * w0md * w0t + fi[2] * w0d * w0mt + - fi[3] * w0md * w0mt + fi[4] * w0d * w1t + fi[5] * w0md * w1t + - fi[6] * w0d * w1mt + fi[7] * w0md * w1mt + fi[8] * w0d * w2t + - fi[9] * w0md * w2t + fi[10] * w0d * w2mt + fi[11] * w0md * w2mt + - fi[12] * w1d * w0t + fi[13] * w1md * w0t + fi[14] * w1d * w0mt + - fi[15] * w1md * w0mt + fi[16] * w2d * w0t + fi[17] * w2md * w0t + - fi[18] * w2d * w0mt + fi[19] * w2md * w0mt + fi[20] * w1d * w1t + - fi[21] * w1md * w1t + fi[22] * w1d * w1mt + fi[23] * w1md * w1mt + - fi[24] * w2d * w1t + fi[25] * w2md * w1t + fi[26] * w2d * w1mt + - fi[27] * w2md * w1mt + fi[28] * w1d * w2t + fi[29] * w1md * w2t + - fi[30] * w1d * w2mt + fi[31] * w1md * w2mt + fi[32] * w2d * w2t + - fi[33] * w2md * w2t + fi[34] * w2d * w2mt + fi[35] * w2md * w2mt; -} - -// cubic hermite polynomial -// psi0 and its derivatives -PORTABLE_INLINE_FUNCTION Real xpsi0(Real z) { return z * z * (2.0 * z - 3.0) + 1.0; } -PORTABLE_INLINE_FUNCTION Real xdpsi0(Real z) { return z * (6.0 * z - 6.0); } - -// psi1 and its derivatives -PORTABLE_INLINE_FUNCTION Real xpsi1(Real z) { return z * (z * (z - 2.0) + 1.0); } - -PORTABLE_INLINE_FUNCTION Real xdpsi1(Real z) { return z * (3.0 * z - 4.0) + 1.0; } - -// bicubic hermite polynomial -PORTABLE_INLINE_FUNCTION Real h3(Real fi[16], Real w0t, Real w1t, Real w0mt, Real w1mt, - Real w0d, Real w1d, Real w0md, Real w1md) { - return fi[0] * w0d * w0t + fi[1] * w0md * w0t + fi[2] * w0d * w0mt + - fi[3] * w0md * w0mt + fi[4] * w0d * w1t + fi[5] * w0md * w1t + - fi[6] * w0d * w1mt + fi[7] * w0md * w1mt + fi[8] * w1d * w0t + - fi[9] * w1md * w0t + fi[10] * w1d * w0mt + fi[11] * w1md * w0mt + - fi[12] * w1d * w1t + fi[13] * w1md * w1t + fi[14] * w1d * w1mt + - fi[15] * w1md * w1mt; -} - -} // namespace hermite -} // namespace singularity - -#endif // SINGULARITY_EOS_BASE_MATH_UTILS_HPP_ +#endif // SINGULARITY_EOS_BASE_HERMITE_HPP_ diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp index c950832f797..c62ea08e9f0 100644 --- a/singularity-eos/base/indexable_types.hpp +++ b/singularity-eos/base/indexable_types.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2025. 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 @@ -15,341 +15,6 @@ #ifndef SINGULARITY_EOS_BASE_INDEXABLE_TYPES_ #define SINGULARITY_EOS_BASE_INDEXABLE_TYPES_ -#include -#include -#include +#include -#include -#include -#include - -namespace singularity { -namespace IndexerUtils { - -// Disable unreachable code warnings for SafeSet and SafeGet -#if defined(__clang__) -#pragma clang diagnostic push -#pragma clang diagnostic ignored "-Wunreachable-code" -#elif defined(__GNUC__) -#pragma GCC diagnostic push -#pragma GCC diagnostic ignored "-Wunreachable-code" -#elif defined(_MSC_VER) -#pragma warning(push) -#pragma warning(disable : 4702) // unreachable code -#endif - -// Identifies an indexer as a type-based indexer -template -struct is_type_indexer : std::false_type {}; -template -struct is_type_indexer::is_type_indexable)>> - : std::bool_constant::is_type_indexable> {}; -template -constexpr bool is_type_indexer_v = is_type_indexer::value; - -namespace impl { - -// Simple way to switch between pure type indexing or also allowing intergers -enum class AllowedIndexing { Numeric, TypeOnly }; - -// The "safe" version of Get(). This function will ONLY return a value IF that -// type-based index is present in the Indexer OR if the Indexer doesn't support -// type-based indexing. -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, const T &t, - std::size_t const idx, Real &out) { - // If null then nothing happens - if (variadic_utils::is_nullptr(lambda)) { - return false; - } - - // Return value if type index is available - if constexpr (variadic_utils::is_indexable_v) { - out = lambda[t]; - return true; - } - - // Do nothing if lambda has type indexing BUT doesn't have this type index - if constexpr (is_type_indexer_v) { - return false; - } - - // Fall back to numeric indexing if allowed - if constexpr (AI == AllowedIndexing::Numeric) { - if constexpr (variadic_utils::has_int_index_v) { - out = lambda[idx]; - return true; - } - } - - // Something else... - return false; -} -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t const idx, - Real &out) { - return SafeGet(lambda, T{}, idx, out); -} - -// Break out "Set" functionality from "Get". The original "Get()" did both, but -// the "safe" version needs to separate that functionality for setting the -// values in a lambda -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, const T &t, - std::size_t const idx, Real const in) { - // If null then nothing happens - if (variadic_utils::is_nullptr(lambda)) { - return false; - } - - // Return value if type index is available - if constexpr (variadic_utils::is_indexable_v) { - lambda[t] = in; - return true; - } - - // Do nothing if lambda has type indexing BUT doesn't have this type index - if constexpr (is_type_indexer_v) { - return false; - } - - // Fall back to numeric indexing if allowed - if constexpr (AI == AllowedIndexing::Numeric) { - if constexpr (variadic_utils::has_int_index_v) { - lambda[idx] = in; - return true; - } - } - - // Something else... - return false; -} -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const idx, - Real const in) { - return SafeSet(lambda, T{}, idx, in); -} - -// Warnings are back after this point -#if defined(__clang__) -#pragma clang diagnostic pop -#elif defined(__GNUC__) -#pragma GCC diagnostic pop -#elif defined(_MSC_VER) -#pragma warning(pop) -#endif - -// Same as above but causes an error condition (static or dynamic) if the value -// can't be obtained. Note that the `decltype(auto)` is intended to preserve the -// value category of the square bracket operator of the `Indexer_t` type. This -// allows references to be returned since there is also no possibility of the -// call doing nothing (i.e. like SafeGet and SafeSet), and thus it can be used -// for either setting or getting values. This should also allow for `const` -// correctness downstream in the wrappers where `lambda` is `const` -template -PORTABLE_FORCEINLINE_FUNCTION decltype(auto) -SafeMustGetSet(Indexer_t &&lambda, const T &t, std::size_t const idx) { - // Error on null pointer - PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda), - "Indexer can't be nullptr"); - - if constexpr (is_type_indexer_v) { - // Return type-based index. Static assert that type MUST exist in indexer - static_assert(variadic_utils::is_indexable_v); - // Use std::forward to maintain value category for lambda, and use - // parentheses to do the same for the output of the lambda[] operation - return (std::forward(lambda)[t]); - } else if constexpr (AI == AllowedIndexing::Numeric) { - // Fall back to numerical indexing if allowed - static_assert(variadic_utils::has_int_index_v); - // Use std::forward to maintain value category for lambda, and use - // parentheses to do the same for the output of the lambda[] operation - return (std::forward(lambda)[idx]); - } else { - // Something else that can't be compiled... - static_assert(variadic_utils::dependent_false_v, - "Indexer must either be designated as type-based through a " - "`is_type_indexable` boolean data member or SafeGet/SafeSet function " - "must be called with a numerical index"); - } -} -template -PORTABLE_FORCEINLINE_FUNCTION decltype(auto) SafeMustGetSet(Indexer_t &&lambda, - std::size_t const idx) { - return SafeMustGetSet(lambda, T{}, idx); -} -} // namespace impl - -// Overload when numerical index is provided -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t const idx, - Real &out) { - return impl::SafeGet(lambda, idx, out); -} -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, const T &t, - std::size_t const idx, Real &out) { - return impl::SafeGet(lambda, t, idx, out); -} - -// Overload when numerical index isn't provided -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, Real &out) { - std::size_t idx = 0; - return impl::SafeGet(lambda, idx, out); -} -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, const T &t, - Real &out) { - std::size_t idx = 0; - return impl::SafeGet(lambda, t, idx, out); -} - -// Overload when numerical index is provided -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const idx, - Real const in) { - return impl::SafeSet(lambda, idx, in); -} -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, const T &t, - std::size_t const idx, Real const in) { - return impl::SafeSet(lambda, t, idx, in); -} - -// Overload when numerical index isn't provided -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, Real const in) { - std::size_t idx = 0; - return impl::SafeSet(lambda, idx, in); -} -template -PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, const T &t, Real const in) { - std::size_t idx = 0; - return impl::SafeSet(lambda, t, idx, in); -} - -// Overload when numerical index is provided -template -PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda, - std::size_t const idx) { - return impl::SafeMustGetSet(lambda, idx); -} -template -PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda, const T &t, - std::size_t const idx) { - return impl::SafeMustGetSet(lambda, t, idx); -} - -// Overload when numerical index isn't provided -template -PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda) { - std::size_t idx = 0; - return impl::SafeMustGetSet(lambda, idx); -} -template -PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda, const T &t) { - std::size_t idx = 0; - return impl::SafeMustGetSet(lambda, t, idx); -} - -// Overload when numerical index is provided -template -PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t const idx, - Real const in) { - impl::SafeMustGetSet(lambda, idx) = in; -} -template -PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, const T &t, - std::size_t const idx, Real const in) { - impl::SafeMustGetSet(lambda, t, idx) = in; -} - -// Overload when numerical index isn't provided -template -PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, Real const in) { - std::size_t idx = 0; - impl::SafeMustGetSet(lambda, idx) = in; -} -template -PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, const T &t, - Real const in) { - std::size_t idx = 0; - impl::SafeMustGetSet(lambda, t, idx) = in; -} - -// This is a convenience struct to easily build a small indexer with -// a set of indexable types. -template -class VariadicIndexerBase { - public: - // Any class that wants to be recognized as indexable (so that we don't - // accidentally fall back to integer indexing when we don't want to) needs to - // include this. - constexpr static bool is_type_indexable = true; - - // JHP: another option for the `is_type_indexable` flag is to take the ADL - // route. Essentially this would involve defining a friend function that - // could be defined in an appropriate namesapce so that theoretically a host - // code could use a TPL container with type-based indexing and allow that - // container to be flagged in our code as acceptable. This seems like a bit - // of a heavy hammer for what we need here though. We can easily change this - // if a TPL provides a type that is being used for this purpose. - - VariadicIndexerBase() = default; - - PORTABLE_FORCEINLINE_FUNCTION - VariadicIndexerBase(const Data_t &data) : data_(data) {} - - template ::value>> - PORTABLE_FORCEINLINE_FUNCTION Real &operator[](const T &t) { - constexpr std::size_t idx = variadic_utils::GetIndexInTL(); - return data_[idx]; - } - - PORTABLE_FORCEINLINE_FUNCTION - Real &operator[](const std::size_t idx) { return data_[idx]; } - - template ::value>> - PORTABLE_FORCEINLINE_FUNCTION const Real &operator[](const T &t) const { - constexpr std::size_t idx = variadic_utils::GetIndexInTL(); - return data_[idx]; - } - - PORTABLE_FORCEINLINE_FUNCTION - const Real &operator[](const std::size_t idx) const { return data_[idx]; } - - static inline constexpr std::size_t size() { return sizeof...(Ts); } - - private: - Data_t data_; -}; -// uses an array -template -using VariadicIndexer = VariadicIndexerBase, Ts...>; -// uses a Real* -template -using VariadicPointerIndexer = VariadicIndexerBase; - -} // namespace IndexerUtils - -namespace IndexableTypes { -struct MeanIonizationState {}; -struct LogDensity {}; -struct LogTemperature {}; -struct MeanAtomicMass {}; -struct MeanAtomicNumber {}; -struct ElectronFraction {}; -struct RootStatus {}; -struct TableStatus {}; -struct MassFractions { - std::size_t n; - PORTABLE_FORCEINLINE_FUNCTION - MassFractions(const std::size_t index_) : n(index_) {} -}; -} // namespace IndexableTypes -} // namespace singularity #endif // SINGULARITY_EOS_BASE_INDEXABLE_TYPES_ diff --git a/singularity-eos/base/math_utils.hpp b/singularity-eos/base/math_utils.hpp index f24fe4dbe81..ddd61e14d55 100644 --- a/singularity-eos/base/math_utils.hpp +++ b/singularity-eos/base/math_utils.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2022. 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 @@ -15,50 +15,6 @@ #ifndef SINGULARITY_EOS_BASE_MATH_UTILS_HPP_ #define SINGULARITY_EOS_BASE_MATH_UTILS_HPP_ -#include - -namespace singularity { -namespace math_utils { - -/* - * Replaces "SQUARE, CUBE, etc - * Use ternary operator instead of if constexpr to enable compile-time tail-recursion - * for various edge cases. Following accelerants are used: - * x^{2*i} = x^i x^i - * x^i = i x^{i - 1} - * x^{-i} = 1/x^i - */ -template -PORTABLE_FORCEINLINE_FUNCTION constexpr Real pow(const Real x) { - // defining Ppos is important. Otherwise compiler is not smart - // enough to terminate recursion. - constexpr int Ppos = (P < 0) ? -P : P; - return (P >= 0) - ? ((Ppos % 2) ? x * pow(x) : pow(x) * pow(x)) - : 1 / pow(x); -} - -template <> -PORTABLE_FORCEINLINE_FUNCTION constexpr Real pow<0>(const Real x) { - return 1.0; -} - -template <> -PORTABLE_FORCEINLINE_FUNCTION constexpr Real pow<1>(const Real x) { - return x; -} - -template -PORTABLE_FORCEINLINE_FUNCTION constexpr auto ipow10() { - return pow

(10); -} - -PORTABLE_FORCEINLINE_FUNCTION auto pow10(const Real x) { - constexpr Real ln10 = 2.30258509299405e+00; - return std::exp(ln10 * x); -} - -} // namespace math_utils -} // namespace singularity +#include #endif // SINGULARITY_EOS_BASE_MATH_UTILS_HPP_ diff --git a/singularity-eos/base/robust_utils.hpp b/singularity-eos/base/robust_utils.hpp index 6c51d89abbd..c33776e5fc9 100644 --- a/singularity-eos/base/robust_utils.hpp +++ b/singularity-eos/base/robust_utils.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2026. 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 @@ -15,66 +15,6 @@ #ifndef SINGULARITY_EOS_BASE_ROBUST_UTILS_HPP_ #define SINGULARITY_EOS_BASE_ROBUST_UTILS_HPP_ -#include -#include -#include -#include - -namespace singularity { -namespace robust { - -template -PORTABLE_FORCEINLINE_FUNCTION constexpr auto SMALL() { - return 10 * std::numeric_limits::min(); -} - -template -PORTABLE_FORCEINLINE_FUNCTION constexpr auto EPS() { - return 10 * std::numeric_limits::epsilon(); -} - -template -PORTABLE_FORCEINLINE_FUNCTION constexpr T min_exp_arg() { - return (std::numeric_limits::min_exponent - 1) * M_LN2; -} - -template -PORTABLE_FORCEINLINE_FUNCTION constexpr T max_exp_arg() { - return std::numeric_limits::max_exponent * M_LN2; -} - -template -PORTABLE_FORCEINLINE_FUNCTION auto make_positive(const T val) { - return std::max(val, EPS()); -} - -PORTABLE_FORCEINLINE_FUNCTION -Real make_bounded(const Real val, const Real vmin, const Real vmax) { - return std::min(std::max(val, vmin + EPS()), vmax * (1.0 - EPS())); -} - -template -PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val) { - if constexpr (std::is_unsigned_v) { - return 1; - } else { - return (T(0) <= val) - (val < T(0)); - } -} - -template -PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) { - return a / (b + sgn(b) * SMALL()); -} - -template -PORTABLE_FORCEINLINE_FUNCTION T safe_arg_exp(const T &x) { - return x < min_exp_arg() ? 0.0 - : x > max_exp_arg() ? std::numeric_limits::infinity() - : std::exp(x); -} - -} // namespace robust -} // namespace singularity +#include #endif // SINGULARITY_EOS_BASE_ROBUST_UTILS_HPP_ diff --git a/singularity-eos/base/root-finding-1d/root_finding.hpp b/singularity-eos/base/root-finding-1d/root_finding.hpp index ceb289e5af5..f0469bca7ed 100644 --- a/singularity-eos/base/root-finding-1d/root_finding.hpp +++ b/singularity-eos/base/root-finding-1d/root_finding.hpp @@ -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,600 +12,9 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ -#ifndef _SINGULARITY_EOS_UTILS_ROOT_FINDING_HPP_ -#define _SINGULARITY_EOS_UTILS_ROOT_FINDING_HPP_ +#ifndef _SINGULARITY_EOS_BASE_ROOT_FINDING_1D_ROOT_FINDING_HPP_ +#define _SINGULARITY_EOS_BASE_ROOT_FINDING_1D_ROOT_FINDING_HPP_ -// Implementation based on gsl root finder API -// Code originally taken from nubhlight, LA-UR-19-20336 -// Miller, Ryan, Dolence, (2019). ApJ Supplements, 241(2), 30. -// TODO: this code still contains many "C"-isms. -// TOOD: if-logic could be removed for speed. Careful, though! -// ~JMM +#include -// #include // debug -// #include -#include -#include -#include -#include -#include - -#define SINGULARITY_ROOT_DEBUG (0) -#define SINGULARITY_ROOT_VERBOSE (0) - -#define SINGULARITY_MY_SIGN(x) (x > 0) - (x < 0) - -namespace RootFinding1D { -constexpr const int SECANT_NITER_MAX{1000}; -constexpr const int BISECT_NITER_MAX{1000}; -constexpr const int BISECT_REG_MAX{1000}; -constexpr const int NEWTON_RAPHSON_NITER_MAX{100}; -enum class Status { SUCCESS = 0, FAIL = 1 }; - -/* - // TODO: Something like this would be nice - friend std::ostream& operator<< (std::ostream& os, const RootCounts& c) { - Real tot = c.total(); - std::vector percents(c.nbins_,0); - for (int i = 0; i < c.nbins_; i++) { - percents[i] = (100.*c.counts_[i] / tot); - } - os << "\n********** ROOT FINDING *********\n" - << " ITERATIONS PERCENTAGE\n"; - for (int i = 0; i < c.more_; i++) { - os << " " - << std::right << std::setw(4) - << i - << " " - << std::right << std::setw(3) - << percents[i] << "\n"; - } - os << " more " - << std::right << std::setw(3) - << percents[c.more_] << "\n" - << "*********************************\n" - << std::endl; - return os; - } -*/ -class RootCounts { - private: - static constexpr std::size_t nbins_{15}; - mutable Real counts_[nbins_]; - - public: - PORTABLE_INLINE_FUNCTION - RootCounts() { - for (std::size_t i{0}; i < nbins_; ++i) - counts_[i] = 0; - } - PORTABLE_INLINE_FUNCTION void reset() { - for (std::size_t i{0}; i < nbins_; ++i) - counts_[i] = 0; - } - PORTABLE_INLINE_FUNCTION void increment(std::size_t i) const { - PORTABLE_REQUIRE(i < nbins_, "Index in bounds"); -#ifdef PORTABILITY_STRATEGY_NONE - counts_[i] += 1; -#endif // PORTABILITY_STRATEGY_NONE - } - PORTABLE_INLINE_FUNCTION Real total() const { - Real tot{1.e-20}; - for (std::size_t i{0}; i < nbins_; ++i) - tot += counts_[i]; - return tot; - } - PORTABLE_INLINE_FUNCTION const Real &operator[](const std::size_t i) const { - PORTABLE_REQUIRE(i < nbins_, "Index in bounds"); - return counts_[i]; - } - PORTABLE_INLINE_FUNCTION Real &operator[](const std::size_t i) { - PORTABLE_REQUIRE(i < nbins_, "Index in bounds"); - return counts_[i]; - } - PORTABLE_INLINE_FUNCTION void print_counts() const { - for (std::size_t i{0}; i < nbins_; ++i) - printf("%e\n", counts_[i]); - } - PORTABLE_INLINE_FUNCTION std::size_t nBins() const { return nbins_; } - PORTABLE_INLINE_FUNCTION std::size_t more() const { return nbins_ - 1; } -}; - -PORTABLE_INLINE_FUNCTION bool check_bracket(const Real ya, const Real yb) { - return (ya * yb <= 0.0); -} - -template -PORTABLE_INLINE_FUNCTION bool set_bracket(const T &f, Real &a, const Real guess, Real &b, - Real &ya, const Real yg, Real &yb, - const bool &verbose = false) { - constexpr std::size_t max_search_depth = 6; - Real dx = b - a; - for (std::size_t level = 0; level < max_search_depth; level++) { - const std::size_t nlev = (1 << level); - for (std::size_t i = 0; i < nlev; i++) { - const Real x = a + (i + 0.5) * dx; - const Real yx = f(x); - if (check_bracket(yx, yg)) { - if (x < guess) { - a = x; - ya = yx; - b = guess; - yb = yg; - } else { - a = guess; - ya = yg; - b = x; - yb = yx; - } - return true; - } - } - } - // if we get here then we failed to bracket a root - if (verbose) { - printf("set_bracket failed to bound a root! %.14e %.14e %.14e %.14e %.14e %.14e\n", a, - guess, b, ya, yg, yb); - } - return false; -} - -// solves for f(x,params) - ytarget = 0 -template -PORTABLE_INLINE_FUNCTION Status regula_falsi(const T &f, const Real ytarget, - const Real guess, Real a, Real b, - const Real xtol, const Real ytol, - Real &xroot, - const RootCounts *counts = nullptr, - const bool &verbose = false) { - constexpr std::size_t max_iter = SECANT_NITER_MAX; - auto func = [&](const Real x) { return f(x) - ytarget; }; - Real ya = func(a); - Real yg = func(guess); - Real yb; - - if (check_bracket(ya, yg)) { - b = guess; - yb = yg; - } else { - yb = func(b); - if (check_bracket(yg, yb)) { - a = guess; - ya = yg; - } else { - // ya, yg, and yb have the same sign - if (!set_bracket(func, a, guess, b, ya, yg, yb, verbose)) { - if (verbose) { - printf("regula_falsi failed! %.14e %.14e %.14e %.14e\n", ytarget, guess, a, b); - } - return Status::FAIL; - } - } - } - - Real sign = (ya < 0 ? 1.0 : -1.0); - ya *= sign; - yb *= sign; - - std::size_t b1 = 0; - std::size_t b2 = 0; - std::size_t iteration_count = 0; - while (b - a > 2.0 * xtol && (std::abs(ya) > ytol || std::abs(yb) > ytol) && - iteration_count < max_iter) { - Real c = (a * yb - b * ya) / (yb - ya); - // guard against roundoff because ya or yb is sufficiently close to zero - if (c == a) { - b = a; - continue; - } else if (c == b) { - a = b; - continue; - } - Real yc = sign * func(c); - if (yc > 0.0) { - b = c; - yb = yc; - b1++; - ya *= (b1 > 1 ? 0.5 : 1.0); - b2 = 0; - } else if (yc < 0.0) { - a = c; - ya = yc; - b2++; - yb *= (b2 > 1 ? 0.5 : 1.0); - b1 = 0; - } else { - a = c; - b = c; - } - iteration_count++; - } - auto status = Status::SUCCESS; - if (iteration_count == max_iter) { - if (verbose) { - printf("root finding reached the maximum number of iterations. likely not " - "converged\n"); - } - status = Status::FAIL; - } - if (counts != nullptr) { - if (iteration_count < counts->nBins()) { - counts->increment(iteration_count); - } else { - counts->increment(counts->more()); - } - } - xroot = 0.5 * (a + b); - return status; -} - -// solves for f(x,params) - ytarget = 0 -// WARNING: this root finding expects a different callable f than the other -// root finding methods. f should return a tuple of (f(x), f'(x)) where f'(x) -// is the derivative of f with respect to x. -template -PORTABLE_INLINE_FUNCTION Status newton_raphson(const T &f, const Real ytarget, - const Real guess, const Real a, - const Real b, const Real ytol, Real &xroot, - const RootCounts *counts = nullptr, - const bool &verbose = false, - const bool &fail_on_bound_root = true) { - - constexpr std::size_t max_iter = NEWTON_RAPHSON_NITER_MAX; - Real _x = guess; - Real _xold = 0.0; - auto status = Status::SUCCESS; - - Real yg; - Real dfunc; - - std::size_t iter; - - for (iter = 0; iter < max_iter; iter++) { - std::tie(yg, dfunc) = f(_x); // C++11 tuple unpacking - - // check if we are converged already - if (std::abs(yg - ytarget) < std::abs(ytol * ytarget)) break; - - // not converged; compute the next step - _xold = _x; - _x = _x - (yg - ytarget) / dfunc; - - // check if we are out of bounds - // CAUTION: we do not set the root to the boundary value in this case - // because one might want to handle this on a case by case basis - // (e.g. if the boundary is a physical boundary, then one might want to - // set the root to the boundary value). - // Per default, we fail if the root is out of bounds controlled by - // fail_on_bound_root. - if ((_x <= a && _xold <= a) || (_x >= b && _xold >= b)) { - if (verbose) { - printf("newton_raphson out of bounds! %.14e %.14e %.14e %.14e\n", ytarget, guess, - a, b); - } - if (fail_on_bound_root) { - status = Status::FAIL; - } - break; - } - _x = std::max(std::min(_x, b), a); - } - if (iter >= max_iter) { - if (verbose) { - printf("root finding reached the maximum number of iterations. likely not " - "converged\n"); - } - status = Status::FAIL; - } - - if (counts != nullptr) { - if (iter < counts->nBins()) { - counts->increment(iter); - } else { - counts->increment(counts->more()); - } - } - xroot = _x; - return status; -} - -template -PORTABLE_INLINE_FUNCTION Status findRoot(const T &f, const Real ytarget, Real xguess, - const Real xmin, const Real xmax, - const Real xtol, const Real ytol, Real &xroot, - const RootCounts *counts = nullptr) { - Status status; - - // first check if we're at the max or min values - const Real fmax = f(xmax); - const Real errmax = fabs(fmax - ytarget) / (fabs(fmax) + ytol); - if (errmax < ytol) { - xroot = xmax; - if (counts != nullptr) { - counts->increment(0); - } - return Status::SUCCESS; - } - const Real fmin = f(xmin); - const Real errmin = fabs(fmin - ytarget) / (fabs(fmin) + ytol); - if (errmin < ytol) { - xroot = xmin; - if (counts != nullptr) { - counts->increment(0); - } - return Status::SUCCESS; - } - - if (xguess >= xmax) xguess = xmax - xtol; - if (xguess <= xmin) xguess = xmin + xtol; - - // Next try Secant - status = secant(f, ytarget, xguess, xmin, xmax, xtol, ytol, xroot, counts); - if (status == Status::SUCCESS) return status; - -#if SINGULARITY_ROOT_DEBUG - if (isnan(xroot)) { - fprintf(stderr, "xroot is nan after secant\n"); - } -#endif - -// Secant failed. Try bisection. -#if SINGULARITY_ROOT_VERBOSE - fprintf(stderr, - "\n\nRoot finding. Secant failed. Trying bisection.\n" - "\txguess = %.10g\n" - "\tytarget = %.10g\n" - "\txmin = %.10g\n" - "\txmax = %.10g\n", - xguess, ytarget, xmin, xmax); -#endif - status = bisect(f, ytarget, xguess, xmin, xmax, xtol, ytol, xroot, counts); - - // Check for something horrible happening - if (std::isnan(xroot) || std::isinf(xroot)) { -#if SINGULARITY_ROOT_DEBUG - fprintf(stderr, "xroot is nan after bisection\n"); -#endif - return Status::FAIL; - } - if (xroot < xmin) return Status::FAIL; - if (xroot > xmax) return Status::FAIL; - - return status; -} - -template -PORTABLE_INLINE_FUNCTION Status secant(const T &f, const Real ytarget, const Real xguess, - const Real xmin, const Real xmax, const Real xtol, - const Real ytol, Real &xroot, - const RootCounts *counts = nullptr) { - Real dx; - Real x_last, y, yp, ym, dyNum, dyDen, dy; - - Real x = xguess; - std::size_t iter{0}; - for (iter = 0; iter < SECANT_NITER_MAX; ++iter) { - x_last = x; - dx = fabs(1.e-7 * x) + xtol; - y = f(x) - ytarget; - yp = f(x + dx); - ym = f(x - dx); - dyNum = yp - ym; - dyDen = (2. * dx); - dy = dyNum / dyDen; - x -= y / dy; - if (x < xmin) x = xmin; - if (x > xmax) x = xmax; - if (std::isnan(x) || std::isinf(x)) { - // can't recover from this -#if SINGULARITY_ROOT_DEBUG - fprintf(stderr, - "\n\n[secant]: NAN or out-of-bounds detected!\n" - "\txguess = %.10e\n" - "\tytarget = %.10e\n" - "\tx = %.10e\n" - "\tx_last = %.10e\n" - "\tx_min = %.10e\n" - "\tx_max = %.10e\n" - "\ty = %.10e\n" - "\tdx = %.10e\n" - "\typ = %.10e\n" - "\tym = %.10e\n" - "\tdyNum = %.10e\n" - "\tdyDen = %.10e\n" - "\tdy = %.10e\n" - "\titer = %d\n" - "\tsign x = %d\n", - xguess, ytarget, x, x_last, xmin, xmax, y, dx, yp, ym, dyNum, dyDen, dy, - iter, (int)SINGULARITY_MY_SIGN(x)); -#endif - if (counts != nullptr) { - counts->increment(counts->more()); - } - return Status::FAIL; - } - if (fabs(x - x_last) / (fabs(x) + xtol) < xtol) break; - } - // increment iter - ++iter; - if (counts != nullptr) { - if (iter < counts->nBins()) { - counts->increment(iter); - } else { - counts->increment(counts->more()); - } - } - xroot = x; - - y = f(x); - const Real frac_error = fabs(y - ytarget) / (fabs(y) + ytol); -#if SINGULARITY_ROOT_DEBUG - if (frac_error > ytol) { - fprintf(stderr, - "\n\n[secant]: Failed via too large yerror.\n" - "\tfractional error = %.10e\n" - "\tx = %.10e\n" - "\ty = %.10e\n" - "\tytarget = %.10e\n" - "\typ = %.10e\n" - "\tym = %.10e\n" - "\tdy = %.10e\n" - "\tdx = %.10e\n" - "\titer = %d\n", - frac_error, x, y, ytarget, yp, ym, dy, dx, iter); - } - if (fabs(x - x_last) > xtol) { - fprintf(stderr, - "\n\n[secant]: failed via dx too big.\n" - "\tfractional error = %.10e\n" - "\tx = %.10e\n" - "\tx_last = %.10e\n" - "\tdx = %.10e\n" - "\ty = %.10e\n" - "\tytarget = %.10e\n" - "\typ = %.10e\n" - "\tym = %.10e\n" - "\tdy = %.10e\n" - "\tdx = %.10e\n" - "\titer = %d\n", - frac_error, x, x_last, fabs(x - x_last), y, ytarget, yp, ym, dy, dx, iter); - } -#endif - - const int secant_failed = (std::abs(x - x_last) > xtol || std::abs(frac_error) > ytol || - std::isnan(x) || std::isinf(x)); - return secant_failed ? Status::FAIL : Status::SUCCESS; -} - -template -PORTABLE_INLINE_FUNCTION Status bisect(const T &f, const Real ytarget, const Real xguess, - const Real xmin, const Real xmax, const Real xtol, - const Real ytol, Real &xroot, - const RootCounts *counts = nullptr) { - Real xl, xr, fl, fr, dx; - - Real grow = 0.01; - Real x = xguess; - if (fabs(x) < xtol) { - x += 2. * xtol; - } - // do { // Try to find reasonable region for bisection - for (std::size_t i{0}; i < BISECT_REG_MAX; ++i) { - dx = fabs(grow * x); - xl = x - dx; - xr = x + dx; - fl = f(xl) - ytarget; - fr = f(xr) - ytarget; - grow *= 1.1; - if (fl * fr < 0.0 || xl < xmin || xr > xmax) break; - if (i > BISECT_REG_MAX - 2) { -#if SINGULARITY_ROOT_DEBUG - fprintf(stderr, - "\n\n[Bisect]: expanding region failed\n" - "\txl = %.10g\n" - "\txr = %.10g\n" - "\tfl = %.10g\n" - "\tfr = %.10g\n" - "\tdx = %.10g\n" - "\tgrow = %.10g\n", - xl, xr, fl, fr, dx, grow); -#endif - return Status::FAIL; - } - } - - // force back onto the bisection region - if (xr > xmax) { - xr = xmax; - fr = f(xr) - ytarget; - } - if (xl < xmin) { - xl = xmin; - fl = f(xl) - ytarget; - } - - // if they have the same sign, change that. - // if we can't fix it, fail. - if (fl * fr > 0) { - xl = xmin; - fl = f(xl) - ytarget; - if (fl * fr > 0) { - xr = xmax; - fr = f(xr) - ytarget; - if (fl * fr > 0) { -#if SINGULARITY_ROOT_DEBUG - Real il = f(xl); - Real ir = f(xr); - fprintf(stderr, - "\n\n[bisect]: fl*fr > 0!\n" - "\txguess = %.10e\n" - "\tytarget = %.10e\n" - "\txl = %.10e\n" - "\txr = %.10e\n" - "\tfl = %.10e\n" - "\tfr = %.10e\n" - "\til = %.10e\n" - "\tir = %.10e\n", - xguess, ytarget, xl, xr, fl, fr, il, ir); - std::size_t nx = 300; - Real dx = (xmax - xmin) / (nx - 1); - fprintf(stderr, "Area map:\nx\ty\n"); - for (std::size_t i = 0; i < nx; i++) { - fprintf(stderr, "%.4f\t%.4e\n", x + i * dx, f(x + i * dx)); - } -#endif - return Status::FAIL; - } - } - } - - for (std::size_t i{0}; i < BISECT_NITER_MAX; ++i) { - Real xm = 0.5 * (xl + xr); - Real fm = f(xm) - ytarget; - if (fl * fm <= 0) { - xr = xm; - fr = fm; - } else { - xl = xm; - fl = fm; - } - if (xr - xl < xtol) { - xroot = 0.5 * (xl + xr); - return Status::SUCCESS; - } - } - - xroot = 0.5 * (xl + xr); - - if (std::isnan(xroot)) { -#if SINGULARITY_ROOT_DEBUG - Real il = f(xl); - Real ir = f(xr); - fprintf(stderr, - "\n\n[bisect]: NAN DETECTED!\n" - "\txguess = %.10e\n" - "\tytarget = %.10e\n" - "\txl = %.10e\n" - "\txr = %.10e\n" - "\tdx = %.10e\n" - "\tgrow = %.10e\n" - "\txtol = %.10e\n" - "\tfl = %.10e\n" - "\tfr = %.10e\n" - "\til = %.10e\n" - "\tir = %.10e\n" - "\txmin = %.10e\n" - "\txmax = %.10e\n", - xguess, ytarget, xl, xr, dx, grow, xtol, fl, fr, il, ir, xmin, xmax); -#endif - } - return Status::FAIL; -} -// ---------------------------------------------------------------------- -} // namespace RootFinding1D - -#undef SINGULARITY_ROOT_DEBUG -#undef SINGULARITY_ROOT_VERBOSE -#undef SINGULARITY_MY_SIGN - -#endif // _SINGULARITY_EOS_UTILS_ROOT_FINDING_HPP_ +#endif // _SINGULARITY_EOS_BASE_ROOT_FINDING_1D_ROOT_FINDING_HPP_ diff --git a/singularity-eos/base/sp5/singularity_eos_sp5.hpp b/singularity-eos/base/sp5/singularity_eos_sp5.hpp index 67eedf64da9..c3f25d9393c 100644 --- a/singularity-eos/base/sp5/singularity_eos_sp5.hpp +++ b/singularity-eos/base/sp5/singularity_eos_sp5.hpp @@ -1,7 +1,7 @@ //------------------------------------------------------------------------------ // Extra definitions to the SP5 file format required for singularity-eos // Author: Jonah Miller (jonahm@lanl.gov) -// © 2021-2025. 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 @@ -14,64 +14,9 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ -#ifndef _SINGULARITY_EOS_UTILS_SP5_SINGULARITY_EOS_SP5_HPP_ -#define _SINGULARITY_EOS_UTILS_SP5_SINGULARITY_EOS_SP5_HPP_ +#ifndef _SINGULARITY_EOS_BASE_SP5_SINGULARITY_EOS_SP5_HPP_ +#define _SINGULARITY_EOS_BASE_SP5_SINGULARITY_EOS_SP5_HPP_ -namespace SP5 { +#include -constexpr char defaultSesFileName[] = "sesame_table.sp5"; -constexpr char logType[] = "log_type"; - -namespace Depends { -constexpr char logRhoLogSie[] = "dependsLogRhoLogSie"; -constexpr char logRhoLogT[] = "dependsLogRhoLogT"; -constexpr char coldCurve[] = "coldCurve"; -constexpr char massFrac[] = "massFrac"; -} // namespace Depends - -namespace SubTable { -constexpr char electronOnly[] = "electronOnly"; -constexpr char ionCold[] = "ionCold"; -} // namespace SubTable - -namespace Offsets { -constexpr char messageName[] = "interpretation"; -constexpr char message[] = - "All quantities are functions of log_10(X)\n" - "for X = density rho, temperature T, or specific internal energy sie\n" - "where conversion is X = 10^{Xlog} - Xoffset\n"; -constexpr char rho[] = "rhoOffset"; -constexpr char T[] = "TOffset"; -constexpr char sie[] = "sieOffset"; -} // namespace Offsets - -namespace Material { -constexpr char exchangeCoefficient[] = "exchangeCoefficient"; -constexpr char meanAtomicMass[] = "meanAtomicMass"; -constexpr char meanAtomicNumber[] = "meanAtomicNumber"; -constexpr char solidBulkModulus[] = "solidBulkModulus"; -constexpr char normalDensity[] = "normalDensity"; -constexpr char comments[] = "comments"; -constexpr char matid[] = "matid"; -constexpr char name[] = "name"; -} // namespace Material - -namespace Fields { -constexpr char P[] = "pressure"; -constexpr char sie[] = "specific internal energy"; -constexpr char T[] = "temperature"; -constexpr char bMod[] = "bulk modulus"; -constexpr char dPdRho[] = "dPdRho"; -constexpr char dPdE[] = "dPdE"; -constexpr char dTdRho[] = "dTdRho"; -constexpr char dTdE[] = "dTdE"; -constexpr char dEdRho[] = "dEdRho"; -constexpr char dEdT[] = "dEdT"; -constexpr char mask[] = "mask"; -constexpr char transitionMask[] = "transition mask"; -constexpr char massFrac[] = "mass fractions"; -} // namespace Fields - -} // namespace SP5 - -#endif // _SINGULARITY_EOS_UTILS_SP5_SINGULARITY_EOS_SP5_HPP_ +#endif // _SINGULARITY_EOS_BASE_SP5_SINGULARITY_EOS_SP5_HPP_ diff --git a/singularity-eos/base/spiner_table_utils.hpp b/singularity-eos/base/spiner_table_utils.hpp index 52f5a2c0209..dfeb6a63f44 100644 --- a/singularity-eos/base/spiner_table_utils.hpp +++ b/singularity-eos/base/spiner_table_utils.hpp @@ -24,264 +24,13 @@ #include #include -#include #include -#include +#include // Bounds #include -#include namespace singularity { namespace table_utils { -// For logarithmic interpolation, quantities may be negative. -// If they are, use offset to ensure negative values make sense. -template -class Bounds { - public: - using RegularGrid1D = Spiner::RegularGrid1D; - using Grid_t = Spiner::PiecewiseGrid1D; - // for tag dispatch in constructors - using OneGrid = std::integral_constant; - using TwoGrids = std::integral_constant; - using ThreeGrids = std::integral_constant; - - Bounds() : offset(0), piecewise(false), linmin_(0), linmax_(0) {} - - Bounds(const Real min, const Real max, const int N, const Real offset) - : grid(Grid_t(std::vector{RegularGrid1D(min, max, N)})), - offset(offset), piecewise(false), linmin_(min), linmax_(max) {} - Bounds(OneGrid, const Real min, const Real max, const int N, const Real offset) - : Bounds(min, max, N, offset) {} - - Bounds(Real min, Real max, int N, const bool convertToLog = false, - const Real shrinkRange = 0, - Real anchor_point = std::numeric_limits::signaling_NaN()) - : offset(0), piecewise(true), linmin_(min), linmax_(max) { - if (convertToLog) { - convertBoundsToLog_(min, max, shrinkRange); - if (!(std::isnan(anchor_point))) { - anchor_point = singularity::FastMath::log10(std::abs(anchor_point)); - } - } - if (!(std::isnan(anchor_point))) { - adjustForAnchor_(min, max, N, anchor_point); - } - grid = Grid_t(std::vector{RegularGrid1D(min, max, N)}); - } - Bounds(OneGrid, Real min, Real max, int N, const bool convertToLog = false, - const Real shrinkRange = 0, - Real anchor_point = std::numeric_limits::signaling_NaN()) - : Bounds(min, max, N, convertToLog, shrinkRange, anchor_point) {} - - Bounds(TwoGrids, Real global_min, Real global_max, Real anchor_point, Real splitPoint, - const Real ppd_fine, const Real ppd_factor, const bool convertToLog, - const Real shrinkRange = 0) - : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { - const Real ppd_coarse = (ppd_factor > 0) ? ppd_fine / ppd_factor : ppd_fine; - - if (convertToLog) { - convertBoundsToLog_(global_min, global_max, shrinkRange); - anchor_point = toLog_(anchor_point, offset); - splitPoint = toLog_(splitPoint, offset); - } - - checkInterval_(splitPoint, global_min, global_max, "Split point"); - checkInterval_(anchor_point, global_min, splitPoint, "Anchor point"); - - // add a point just to make sure we have enough points after adjusting for anchor - int N_fine = getNumPointsFromDensity(global_min, splitPoint, ppd_fine) + 1; - adjustForAnchor_(global_min, splitPoint, N_fine, anchor_point); - RegularGrid1D grid_lower(global_min, splitPoint, N_fine); - - const int N_upper = getNumPointsFromDensity(splitPoint, global_max, ppd_coarse); - RegularGrid1D grid_upper(splitPoint, global_max, N_upper); - - grid = Grid_t(std::vector{grid_lower, grid_upper}); - } - - Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, Real fine_min, - Real fine_max, const Real ppd_fine, const Real ppd_factor_lo, - const Real ppd_factor_hi, const bool convertToLog, const Real shrinkRange = 0) - : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { - - if (convertToLog) { - convertBoundsToLog_(global_min, global_max, shrinkRange); - anchor_point = toLog_(anchor_point, offset); - fine_min = toLog_(fine_min, offset); - fine_max = toLog_(fine_max, offset); - } - checkInterval_(anchor_point, global_min, global_max, "Anchor point"); - - grid = gridFromIntervals_(ThreeGrids(), global_min, global_max, anchor_point, - fine_min, fine_max, ppd_fine, ppd_factor_lo, ppd_factor_hi); - } - - Bounds(ThreeGrids, Real global_min, Real global_max, Real anchor_point, - Real log_fine_diameter, const Real ppd_fine, const Real ppd_factor_lo, - const Real ppd_factor_hi, const bool convertToLog, const Real shrinkRange = 0) - : offset(0), piecewise(true), linmin_(global_min), linmax_(global_max) { - - if (convertToLog) { - convertBoundsToLog_(global_min, global_max, shrinkRange); - anchor_point = toLog_(anchor_point, offset); - } - - checkInterval_(anchor_point, global_min, global_max, "Anchor point"); - Real mid_min = anchor_point - 0.5 * log_fine_diameter; - Real mid_max = anchor_point + 0.5 * log_fine_diameter; - - grid = gridFromIntervals_(ThreeGrids(), global_min, global_max, anchor_point, mid_min, - mid_max, ppd_fine, ppd_factor_lo, ppd_factor_hi); - } - - inline Real log2lin(const Real xl) const { - // JMM: Need to guard this with the linear bounds passed in. The - // reason is that our fast math routines, while completely - // invertible at the level of machine epsilon, do introduce error - // at machine epsilon, which can bring us out of the interpolation - // range of eospac. - return std::min(linmax_, - std::max(linmin_, singularity::FastMath::pow10(xl) - offset)); - } - inline Real i2lin(const int i) const { return log2lin(grid.x(i)); } - - friend std::ostream &operator<<(std::ostream &os, const Bounds &b) { - os << "Bounds: [" << b.grid.min() << ", " << b.grid.max() << "]" - << " + " << b.offset << "\n" - << "\tN = " << b.grid.nPoints() << "\n"; - for (int ig = 0; ig < b.grid.nGrids(); ++ig) { - os << "\t[ig,dx] = [" << ig << ", " << b.grid.dx(ig) << "]" - << "\n"; - } - return os; - } - - // This uses real logs - template - static int getNumPointsFromPPD(Real min, Real max, const T ppd) { - constexpr Real epsilon = std::numeric_limits::epsilon(); - const Real min_offset = 10 * std::abs(epsilon); - // 1.1 so that the y-intercept isn't identically zero - // when min < 0. - // min_offset to handle the case where min=0 - Real offset = 0; - if (min <= 0) offset = 1.1 * std::abs(min) + min_offset; - min += offset; - max += offset; - - Real lmin = std::log10(min); - Real lmax = std::log10(max); - int N = getNumPointsFromDensity(lmin, lmax, ppd); - return N; - } - template - static int getNumPointsFromDensity(const Real min, const Real max, const T density) { - Real delta = max - min; - int N = std::max(2, static_cast(std::ceil(density * delta))); - return N; - } - - private: - Real toLog_(Real val, const Real offset) { - val += offset; - return singularity::FastMath::log10(std::abs(val)); - } - - void convertBoundsToLog_(Real &min, Real &max, const Real shrinkRange = 0) { - // Log scales can't handle negative numbers or exactly zero. To - // deal with that, we offset. - constexpr Real epsilon = std::numeric_limits::epsilon(); - const Real min_offset = 10 * std::abs(epsilon); - // 1.1 so that the y-intercept isn't identically zero - // when min < 0. - // min_offset to handle the case where min=0 - if (min <= 0) offset = 1.1 * std::abs(min) + min_offset; - - min = toLog_(min, offset); - max = toLog_(max, offset); - - Real delta = max - min; - min += 0.5 * shrinkRange * delta; - max -= 0.5 * shrinkRange * delta; - } - - static void adjustForAnchor_(const Real min, Real &max, int &N, - const Real anchor_point) { - if (min < anchor_point && anchor_point < max) { - Real Nfrac = (anchor_point - min) / (max - min); - PORTABLE_REQUIRE((0 < Nfrac && Nfrac < 1), "anchor in bounds"); - - int Nanchor = static_cast(std::ceil(N * Nfrac)); - if ((Nanchor < 2) || (Nanchor >= N)) return; // not possible to shift this safely - - Real dx = (anchor_point - min) / static_cast(Nanchor - 1); - int Nmax_new = static_cast((max - min) / dx); - Real max_new = dx * (Nmax_new - 1) + min; - PORTABLE_REQUIRE(max_new <= max, "must not exceed table bounds"); - - N = Nmax_new; - max = max_new; - } else { - PORTABLE_ALWAYS_WARN("Anchor point out of bounds. Ignoring it."); - } - } - - static void checkInterval_(Real &p, const Real min, const Real max, - const std::string &name) { - if (p <= min) { - PORTABLE_ALWAYS_WARN(name + " less than minimum. Adjusting."); - Real eps = 0.1 * std::abs(min); - p = min + eps; - } - if (p >= max) { - PORTABLE_ALWAYS_WARN(name + " greater than maximum. Adjusting."); - Real eps = 0.1 * std::abs(max); - p = max - eps; - } - } - - static Grid_t gridFromIntervals_(ThreeGrids, Real global_min, Real global_max, - Real anchor_point, Real mid_min, Real mid_max, - const Real ppd_fine, const Real ppd_factor_lo, - const Real ppd_factor_hi) { - const Real ppd_lo = (ppd_factor_lo > 0) ? ppd_fine / ppd_factor_lo : ppd_fine; - const Real ppd_hi = (ppd_factor_hi > 0) ? ppd_fine / ppd_factor_hi : ppd_fine; - - if (mid_min <= global_min) { - PORTABLE_ALWAYS_WARN( - "Table bounds refined minimum lower than global minimum. Adjusting."); - Real delta = std::abs(anchor_point - global_min); - mid_min = anchor_point - 0.9 * delta; - } - if (mid_max >= global_max) { - PORTABLE_ALWAYS_WARN( - "Table bounds refined maximum greater than global maximum. Adjusting."); - Real delta = std::abs(global_max - anchor_point); - mid_max = anchor_point + 0.9 * delta; - } - - // add a point just to make sure we have enough points after adjusting for anchor - int N_fine = getNumPointsFromDensity(mid_min, mid_max, ppd_fine) + 1; - adjustForAnchor_(mid_min, mid_max, N_fine, anchor_point); - RegularGrid1D grid_middle(mid_min, mid_max, N_fine); - - const int N_lower = getNumPointsFromDensity(global_min, mid_min, ppd_lo); - RegularGrid1D grid_lower(global_min, mid_min, N_lower); - - const int N_upper = getNumPointsFromDensity(mid_max, global_max, ppd_hi); - RegularGrid1D grid_upper(mid_max, global_max, N_upper); - - return Grid_t(std::vector{grid_lower, grid_middle, grid_upper}); - } - - public: - Grid_t grid; - Real offset = 0; - bool piecewise = false; - - private: - Real linmin_, linmax_; -}; // JMM: Making this a struct with static methods, rather than a // namespace, saves a few "friend" declarations. Broadly these methods diff --git a/singularity-eos/base/variadic_utils.hpp b/singularity-eos/base/variadic_utils.hpp index 8b30c0cb5b2..b6c8698c971 100644 --- a/singularity-eos/base/variadic_utils.hpp +++ b/singularity-eos/base/variadic_utils.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2024. 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 @@ -15,250 +15,6 @@ #ifndef SINGULARITY_EOS_BASE_VARIADIC_UTILS_HPP_ #define SINGULARITY_EOS_BASE_VARIADIC_UTILS_HPP_ -#include -#include - -namespace singularity { -namespace variadic_utils { - -// Some generic variatic utilities -// ====================================================================== - -// Template parameter dependent boolean suitable for causing static_assert -// errors within `if constexpr` branches. Essentially the issue is that if the -// static_assert _always_ evaluates to false, then it will _always_ cause a -// compile time error even if that branch of the code will never be reached. -// Making the evaluation (superficially) dependent on the template deduction -// causes it to be evaluated after the `if constexpr` branching has already been -// determined. See https://en.cppreference.com/w/cpp/language/if.html#Constexpr_if -template -inline constexpr bool dependent_false_v = false; - -// Useful for generating nullptr of a specific pointer type -template -inline constexpr T *np() { - return nullptr; -} - -// C++14 implementation of std::remove_cvref (available since C++20) -// credit to CJ + Diego -template -struct remove_cvref { - typedef std::remove_cv_t> type; -}; - -// Helper types -template -using remove_cvref_t = typename remove_cvref::type; - -// SFINAE to check if a value is a null ptr -template >::value>::type> -constexpr inline bool is_nullptr(T &&t) { - return std::forward(t) == nullptr; -} -template >::value>::type * = nullptr> -constexpr inline bool is_nullptr(T &&) { - return false; -} - -// Backport of C++17 bool_constant. -// With C++17, can be replaced with -// using std::bool_constant -template -using bool_constant = std::integral_constant; - -// Implementation of std::conjunction/std::disjunction without C++17 -// With C++17, can be replaced with -// using std::disjunction -template -struct bool_pack {}; -template -using conjunction = std::is_same, bool_pack>; -template -struct disjunction : bool_constant::value> {}; - -// Checks if T is contained in the pack Ts -template -using contains = disjunction::value...>; - -template -constexpr bool contains_v() { - return contains::value; -} - -// variadic list -template -struct type_list {}; - -// variadic list of modifiers -template