diff --git a/CMakeLists.txt b/CMakeLists.txt index 87b8789d101..27a34074c16 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -352,6 +352,8 @@ list(APPEND libopenmc_SOURCES src/event.cpp src/file_utils.cpp src/finalize.cpp + src/gendf.cpp + src/gendf_parser.cpp src/geometry.cpp src/geometry_aux.cpp src/hdf5_interface.cpp diff --git a/bin/openmc b/bin/openmc new file mode 100755 index 00000000000..8d6dfab9426 Binary files /dev/null and b/bin/openmc differ diff --git a/geometry.xml b/geometry.xml new file mode 100644 index 00000000000..b5e29c3bbaa --- /dev/null +++ b/geometry.xml @@ -0,0 +1,65 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.0 1.0 +
0.0 0.0 0.0
+ + 9 + 20 10 +19 3 11 + 8 4 +18 2 12 + 7 5 +17 6 13 + 16 14 + 15 + 29 + 40 30 +39 23 31 + 28 24 +38 22 32 + 27 25 +37 26 33 + 36 34 + 35 + + diff --git a/include/openmc/capi.h b/include/openmc/capi.h index d8041ef4149..ac2df4c70c1 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -275,6 +275,79 @@ int openmc_properties_export(const char* filename); // \return Error code int openmc_properties_import(const char* filename); +//============================================================================== +// GENDF Library Functions +//============================================================================== + +//! Create GENDF library instance +//! \param[in] library_path Path to GENDF library directory +//! \param[in] n_energy_bounds Number of energy boundaries +//! \param[in] energy_bounds Energy group boundaries in eV +//! \param[in] energy_structure_name Energy structure name (e.g., "CCFE-709") +//! \param[out] lib_id Library instance ID +//! \return Error code +int openmc_gendf_library_create(const char* library_path, int n_energy_bounds, + const double* energy_bounds, const char* energy_structure_name, + int32_t* lib_id); + +//! Free GENDF library instance +//! \param[in] lib_id Library ID +//! \return Error code +int openmc_gendf_library_free(int32_t lib_id); + +//! Get number of energy groups in GENDF library +//! \param[in] lib_id Library ID +//! \param[out] n_groups Number of groups +//! \return Error code +int openmc_gendf_library_get_n_groups(int32_t lib_id, int* n_groups); + +//! Get energy group boundaries from GENDF library +//! \param[in] lib_id Library ID +//! \param[out] bounds Pointer to energy bounds array +//! \param[out] n Number of boundaries +//! \return Error code +int openmc_gendf_library_get_energy_bounds( + int32_t lib_id, const double** bounds, int* n); + +//! Check if nuclide is available in GENDF library +//! \param[in] lib_id Library ID +//! \param[in] nuclide Nuclide name +//! \param[out] has Whether nuclide is available +//! \return Error code +int openmc_gendf_library_has_nuclide( + int32_t lib_id, const char* nuclide, bool* has); + +//! Get list of available nuclides in GENDF library +//! \param[in] lib_id Library ID +//! \param[out] nuclides Array of nuclide name strings (caller must free) +//! \param[out] n Number of nuclides +//! \return Error code +int openmc_gendf_library_available_nuclides( + int32_t lib_id, char*** nuclides, int* n); + +//! Get cross-section data for nuclide and reaction +//! \param[in] lib_id Library ID +//! \param[in] nuclide Nuclide name +//! \param[in] mt ENDF MT reaction number +//! \param[in] n_energy_bounds Number of energy boundaries +//! \param[in] energy_bounds Energy boundaries for validation +//! \param[out] xs_data Cross-section array (caller must free with +//! openmc_gendf_free_xs) \param[out] n_groups Number of groups returned +//! \return Error code +int openmc_gendf_get_xs(int32_t lib_id, const char* nuclide, int32_t mt, + int n_energy_bounds, const double* energy_bounds, double** xs_data, + int* n_groups); + +//! Free cross-section array allocated by openmc_gendf_get_xs +//! \param[in] xs_data Array to free +void openmc_gendf_free_xs(double* xs_data); + +//! Free nuclide array allocated by openmc_gendf_library_available_nuclides +//! \param[in] nuclides Array of nuclide names to free (may be NULL) +//! \param[in] n Number of nuclides in array +//! \note Safe to call with NULL pointer or n=0 +void openmc_gendf_free_nuclides(char** nuclides, int n); + // Error codes extern int OPENMC_E_UNASSIGNED; extern int OPENMC_E_ALLOCATE; diff --git a/include/openmc/gendf.h b/include/openmc/gendf.h new file mode 100644 index 00000000000..2ee501277c3 --- /dev/null +++ b/include/openmc/gendf.h @@ -0,0 +1,256 @@ +//! \file gendf.h +//! \brief GENDF (Group-averaged ENDF) cross-section library for depletion + +#ifndef OPENMC_GENDF_H +#define OPENMC_GENDF_H + +#include +#include +#include + +#include "openmc/memory.h" // for unique_ptr +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +// GENDF tolerance constants (H1 fix - unified across C++ and Python) +//============================================================================== + +//! Relative tolerance for energy boundary matching +constexpr double GENDF_RTOL_MATCH = 1.0e-6; + +//! Relative tolerance for energy boundary mismatch warnings +constexpr double GENDF_RTOL_WARN = 1.0e-4; + +//! Absolute tolerance for energy matching (only used when values are near zero) +constexpr double GENDF_ATOL = 0.0; + +//============================================================================== +// Parser validation options and results (G9 fix) +//============================================================================== + +//! Validation options for GENDF parsing +struct GENDFParserOptions { + bool validate_za {true}; //!< Check ZA is physically valid + bool validate_xs_positive {true}; //!< Check XS values are non-negative + bool warn_short_lines {true}; //!< Warn about skipped short lines + bool require_mf1_header {true}; //!< Require MF=1, MT=451 header + int min_file_lines {10}; //!< Minimum expected lines in file +}; + +//! Result from parsing with diagnostics +struct GENDFParseResult { + int za {0}; //!< Z*1000 + A + int zam {0}; //!< Z*1000 + A + isomeric state + std::unordered_map> xs_data; //!< MT -> cross-sections + std::unordered_map> energy_data; //!< MT -> energy boundaries + int lines_read {0}; //!< Total lines read + int lines_skipped {0}; //!< Lines skipped (too short) + int negative_xs_count {0}; //!< Count of negative XS clamped + vector warnings; //!< Warning messages + bool success {false}; //!< Parsing succeeded + std::string error_message; //!< Error message if failed +}; + +//============================================================================== +//! Material data from a single GENDF file (MF=3 cross-sections only) +//============================================================================== + +class GENDFMaterial { +public: + // Constructors + GENDFMaterial() = default; + + //! Construct material from GENDF file + //! \param[in] filename Path to GENDF .asc file + explicit GENDFMaterial(const std::string& filename); + + // Methods + + //! Load material data from GENDF file + //! \param[in] filename Path to GENDF .asc file + void load_from_file(const std::string& filename); + + //! Get cross-section data for specific MT number with energy-aware alignment + //! \param[in] mt ENDF MT reaction number + //! \param[in] n_groups Number of energy groups + //! \param[in] library_bounds Library energy group boundaries for threshold alignment + //! \return Vector of cross-section values (one per group) + vector get_xs(int mt, int n_groups, const vector& library_bounds) const; + + //! Check if MT reaction exists in this material + //! \param[in] mt ENDF MT reaction number + //! \return True if reaction exists + bool has_mt(int mt) const; + + //! Get energy boundaries for specific MT number + //! \param[in] mt ENDF MT reaction number + //! \return Vector of energy boundaries for this reaction + const vector& get_energies(int mt) const; + + //! Check if MT has energy data + //! \param[in] mt ENDF MT reaction number + //! \return True if energy data exists + bool has_energies(int mt) const; + + // Accessors + const std::string& nuclide_name() const { return nuclide_name_; } + int za() const { return za_; } //!< Z*1000 + A + int zam() const { return zam_; } //!< Z*1000 + A + isomeric state + +private: + // Data members + std::string nuclide_name_; //!< Nuclide name (e.g., "U235") + int za_ {0}; //!< Z*1000 + A + int zam_ {0}; //!< Z*1000 + A + isomeric state + + //! Cross-section data: map MT -> vector (one value per group) + std::unordered_map> xs_data_; + + //! Energy data: map MT -> vector (energy boundaries for threshold alignment) + std::unordered_map> energy_data_; + + // Parsing methods + + //! Parse GENDF file (MF=3 only for speed) + //! \param[in] filename Path to GENDF .asc file + void parse_mf3_only(const std::string& filename); +}; + +//============================================================================== +//! GENDF library manager with caching +//============================================================================== + +class GENDFLibrary { +public: + // Constructors + + //! Construct GENDF library + //! \param[in] library_path Path to directory containing GENDF .asc files + //! \param[in] energy_bounds Energy group boundaries in [eV] + //! \param[in] energy_structure_name Optional name of energy structure + explicit GENDFLibrary( + const std::string& library_path, + const vector& energy_bounds, + const std::string& energy_structure_name = ""); + + // Methods + + //! Get cross-section for nuclide and MT + //! \param[in] nuclide Nuclide name (e.g., "U235") + //! \param[in] mt ENDF MT reaction number + //! \param[in] energy_bounds Energy group boundaries in [eV] + //! \return Vector of cross-section values (one per group) + vector get_xs( + const std::string& nuclide, + int mt, + const vector& energy_bounds); + + //! Check if nuclide is available in library + //! \param[in] nuclide Nuclide name + //! \return True if nuclide file exists + bool has_nuclide(const std::string& nuclide) const; + + //! Get list of available nuclides + //! \return Vector of nuclide names + vector available_nuclides() const; + + // Accessors + const vector& energy_bounds() const { return energy_bounds_; } + int n_groups() const { return n_groups_; } + const std::string& energy_structure() const { return energy_structure_; } + const std::string& library_path() const { return library_path_; } + +private: + // Data members + std::string library_path_; //!< Path to GENDF library directory + std::string energy_structure_; //!< Energy group structure name + vector energy_bounds_; //!< Energy group boundaries in [eV] + int n_groups_ {0}; //!< Number of energy groups + + //! Material cache: map nuclide name -> GENDFMaterial + std::unordered_map> material_cache_; + + //! Mutex for thread-safe cache access (mutable for const method compatibility) + //! Uses shared_mutex for read/write locking: multiple readers OR single writer + mutable std::shared_mutex cache_mutex_; + + //! File index: map nuclide name -> file path + std::unordered_map file_index_; + + // Private methods + + //! Load material (with caching) + //! \param[in] nuclide Nuclide name + //! \return Reference to cached material + GENDFMaterial& load_material(const std::string& nuclide); + + //! Build file index by scanning library directory + void build_file_index(); + + //! Get file path for nuclide + //! \param[in] nuclide Nuclide name + //! \return Full path to GENDF file + std::string get_file_path(const std::string& nuclide) const; +}; + +//============================================================================== +// Global variables +//============================================================================== + +namespace data { + +//! Map of library ID -> GENDFLibrary instance +extern std::unordered_map> gendf_libraries; + +//! Next available library ID +extern int n_gendf_libraries; + +} // namespace data + +//============================================================================== +// Non-member functions +//============================================================================== + +//! Strip leading zeros from mass number in nuclide name +//! \param[in] name Nuclide name with possible leading zeros (e.g., "Al027") +//! \return Name with leading zeros stripped (e.g., "Al27") +std::string strip_mass_leading_zeros(const std::string& name); + +//! Convert GENDF filename stem to OpenMC nuclide name +//! Handles metastable suffixes: mg->_m1, ng->_m2, og->_m3, pg->_m4, qg->_m5, g->ground +//! \param[in] stem GENDF filename stem (e.g., "U235g", "Am242mg") +//! \return OpenMC nuclide name (e.g., "U235", "Am242_m1") +std::string convert_gendf_to_openmc_name(const std::string& stem); + +//! Parse GENDF file using MF=3-only parser (7.6x faster) +//! \param[in] filename Path to GENDF .asc file +//! \param[out] xs_data Map of MT -> vector (cross-sections) +//! \param[out] energy_data Map of MT -> vector (energy boundaries) +//! \param[out] za Z*1000 + A +//! \param[out] zam Z*1000 + A + isomeric state +void parse_gendf_mf3_only( + const std::string& filename, + std::unordered_map>& xs_data, + std::unordered_map>& energy_data, + int& za, + int& zam); + +//! Parse GENDF file with validation and diagnostics (G9 fix) +//! \param[in] filename Path to GENDF .asc file +//! \param[in] options Parser validation options +//! \return Parse result with diagnostics +GENDFParseResult parse_gendf_validated( + const std::string& filename, + const GENDFParserOptions& options = GENDFParserOptions{}); + +//! Validate ZA value is physically reasonable +//! \param[in] za Z*1000 + A value +//! \param[out] error Error message if invalid +//! \return True if valid +bool validate_za(int za, std::string& error); + +} // namespace openmc + +#endif // OPENMC_GENDF_H diff --git a/include/openmc/version.h b/include/openmc/version.h new file mode 100644 index 00000000000..1d254f1c2e0 --- /dev/null +++ b/include/openmc/version.h @@ -0,0 +1,21 @@ +#ifndef OPENMC_VERSION_H +#define OPENMC_VERSION_H + +#include "openmc/array.h" + +namespace openmc { + +// OpenMC major, minor, and release numbers +// clang-format off +constexpr int VERSION_MAJOR {0}; +constexpr int VERSION_MINOR {0}; +constexpr int VERSION_RELEASE {0}; +constexpr bool VERSION_DEV {false}; +constexpr const char* VERSION_COMMIT_COUNT = ""; +constexpr const char* VERSION_COMMIT_HASH = "e8d61dd7e9ce918e98ec0b024b5a8bdf46f6bc36"; +constexpr std::array VERSION {VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE}; +// clang-format on + +} // namespace openmc + +#endif // OPENMC_VERSION_H diff --git a/lib64/cmake/OpenMC/OpenMCConfig.cmake b/lib64/cmake/OpenMC/OpenMCConfig.cmake new file mode 100644 index 00000000000..93caffb2a7b --- /dev/null +++ b/lib64/cmake/OpenMC/OpenMCConfig.cmake @@ -0,0 +1,37 @@ +get_filename_component(OpenMC_CMAKE_DIR "${CMAKE_CURRENT_LIST_FILE}" DIRECTORY) + +# Compute the install prefix from this file's location +get_filename_component(_OPENMC_PREFIX "${OpenMC_CMAKE_DIR}/../../.." ABSOLUTE) + +find_package(fmt CONFIG REQUIRED HINTS ${_OPENMC_PREFIX}) +find_package(pugixml CONFIG REQUIRED HINTS ${_OPENMC_PREFIX}) +find_package(xtl CONFIG REQUIRED HINTS ${_OPENMC_PREFIX}) +find_package(xtensor CONFIG REQUIRED HINTS ${_OPENMC_PREFIX}) +if(OFF) + find_package(DAGMC REQUIRED HINTS ) +endif() + +if(OFF) + include(FindPkgConfig) + list(APPEND CMAKE_PREFIX_PATH ) + set(PKG_CONFIG_USE_CMAKE_PREFIX_PATH True) + pkg_check_modules(LIBMESH REQUIRED >=1.7.0 IMPORTED_TARGET) +endif() + +find_package(PNG) + +if(NOT TARGET OpenMC::libopenmc) + include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake") +endif() + +if(on) + find_package(MPI REQUIRED) +endif() + +if(ON) + find_package(OpenMP REQUIRED) +endif() + +if(OFF AND NOT ${DAGMC_BUILD_UWUW}) + message(FATAL_ERROR "UWUW is enabled in OpenMC but the DAGMC installation discovered was not configured with UWUW.") +endif() diff --git a/lib64/cmake/OpenMC/OpenMCConfigVersion.cmake b/lib64/cmake/OpenMC/OpenMCConfigVersion.cmake new file mode 100644 index 00000000000..dc7293852fe --- /dev/null +++ b/lib64/cmake/OpenMC/OpenMCConfigVersion.cmake @@ -0,0 +1,11 @@ +set(PACKAGE_VERSION "0.0.0") + +# Check whether the requested PACKAGE_FIND_VERSION is compatible +if("${PACKAGE_VERSION}" VERSION_LESS "${PACKAGE_FIND_VERSION}") + set(PACKAGE_VERSION_COMPATIBLE FALSE) +else() + set(PACKAGE_VERSION_COMPATIBLE TRUE) + if ("${PACKAGE_VERSION}" VERSION_EQUAL "${PACKAGE_FIND_VERSION}") + set(PACKAGE_VERSION_EXACT TRUE) + endif() +endif() diff --git a/lib64/cmake/OpenMC/OpenMCTargets-relwithdebinfo.cmake b/lib64/cmake/OpenMC/OpenMCTargets-relwithdebinfo.cmake new file mode 100644 index 00000000000..d9b00607926 --- /dev/null +++ b/lib64/cmake/OpenMC/OpenMCTargets-relwithdebinfo.cmake @@ -0,0 +1,28 @@ +#---------------------------------------------------------------- +# Generated CMake target import file for configuration "RelWithDebInfo". +#---------------------------------------------------------------- + +# Commands may need to know the format version. +set(CMAKE_IMPORT_FILE_VERSION 1) + +# Import target "OpenMC::openmc" for configuration "RelWithDebInfo" +set_property(TARGET OpenMC::openmc APPEND PROPERTY IMPORTED_CONFIGURATIONS RELWITHDEBINFO) +set_target_properties(OpenMC::openmc PROPERTIES + IMPORTED_LOCATION_RELWITHDEBINFO "${_IMPORT_PREFIX}/bin/openmc" + ) + +list(APPEND _cmake_import_check_targets OpenMC::openmc ) +list(APPEND _cmake_import_check_files_for_OpenMC::openmc "${_IMPORT_PREFIX}/bin/openmc" ) + +# Import target "OpenMC::libopenmc" for configuration "RelWithDebInfo" +set_property(TARGET OpenMC::libopenmc APPEND PROPERTY IMPORTED_CONFIGURATIONS RELWITHDEBINFO) +set_target_properties(OpenMC::libopenmc PROPERTIES + IMPORTED_LOCATION_RELWITHDEBINFO "${_IMPORT_PREFIX}/lib64/libopenmc.so" + IMPORTED_SONAME_RELWITHDEBINFO "libopenmc.so" + ) + +list(APPEND _cmake_import_check_targets OpenMC::libopenmc ) +list(APPEND _cmake_import_check_files_for_OpenMC::libopenmc "${_IMPORT_PREFIX}/lib64/libopenmc.so" ) + +# Commands beyond this point should not need to know the version. +set(CMAKE_IMPORT_FILE_VERSION) diff --git a/lib64/cmake/OpenMC/OpenMCTargets.cmake b/lib64/cmake/OpenMC/OpenMCTargets.cmake new file mode 100644 index 00000000000..021e47148f4 --- /dev/null +++ b/lib64/cmake/OpenMC/OpenMCTargets.cmake @@ -0,0 +1,116 @@ +# Generated by CMake + +if("${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}" LESS 2.8) + message(FATAL_ERROR "CMake >= 2.8.12 required") +endif() +if(CMAKE_VERSION VERSION_LESS "2.8.12") + message(FATAL_ERROR "CMake >= 2.8.12 required") +endif() +cmake_policy(PUSH) +cmake_policy(VERSION 2.8.12...3.29) +#---------------------------------------------------------------- +# Generated CMake target import file. +#---------------------------------------------------------------- + +# Commands may need to know the format version. +set(CMAKE_IMPORT_FILE_VERSION 1) + +# Protect against multiple inclusion, which would fail when already imported targets are added once more. +set(_cmake_targets_defined "") +set(_cmake_targets_not_defined "") +set(_cmake_expected_targets "") +foreach(_cmake_expected_target IN ITEMS OpenMC::openmc OpenMC::libopenmc) + list(APPEND _cmake_expected_targets "${_cmake_expected_target}") + if(TARGET "${_cmake_expected_target}") + list(APPEND _cmake_targets_defined "${_cmake_expected_target}") + else() + list(APPEND _cmake_targets_not_defined "${_cmake_expected_target}") + endif() +endforeach() +unset(_cmake_expected_target) +if(_cmake_targets_defined STREQUAL _cmake_expected_targets) + unset(_cmake_targets_defined) + unset(_cmake_targets_not_defined) + unset(_cmake_expected_targets) + unset(CMAKE_IMPORT_FILE_VERSION) + cmake_policy(POP) + return() +endif() +if(NOT _cmake_targets_defined STREQUAL "") + string(REPLACE ";" ", " _cmake_targets_defined_text "${_cmake_targets_defined}") + string(REPLACE ";" ", " _cmake_targets_not_defined_text "${_cmake_targets_not_defined}") + message(FATAL_ERROR "Some (but not all) targets in this export set were already defined.\nTargets Defined: ${_cmake_targets_defined_text}\nTargets not yet defined: ${_cmake_targets_not_defined_text}\n") +endif() +unset(_cmake_targets_defined) +unset(_cmake_targets_not_defined) +unset(_cmake_expected_targets) + + +# Compute the installation prefix relative to this file. +get_filename_component(_IMPORT_PREFIX "${CMAKE_CURRENT_LIST_FILE}" PATH) +get_filename_component(_IMPORT_PREFIX "${_IMPORT_PREFIX}" PATH) +get_filename_component(_IMPORT_PREFIX "${_IMPORT_PREFIX}" PATH) +get_filename_component(_IMPORT_PREFIX "${_IMPORT_PREFIX}" PATH) +if(_IMPORT_PREFIX STREQUAL "/") + set(_IMPORT_PREFIX "") +endif() + +# Create imported target OpenMC::openmc +add_executable(OpenMC::openmc IMPORTED) + +set_target_properties(OpenMC::openmc PROPERTIES + INTERFACE_COMPILE_FEATURES "cxx_std_17" +) + +# Create imported target OpenMC::libopenmc +add_library(OpenMC::libopenmc SHARED IMPORTED) + +set_target_properties(OpenMC::libopenmc PROPERTIES + INTERFACE_COMPILE_DEFINITIONS "OPENMC_MPI" + INTERFACE_COMPILE_FEATURES "cxx_std_17" + INTERFACE_INCLUDE_DIRECTORIES "${_IMPORT_PREFIX}/include;/home/perry/Codes/HDF5/HDF5-1.10.9-Linux/HDF_Group/HDF5/1.10.9/include" + INTERFACE_LINK_LIBRARIES "hdf5::hdf5-shared;hdf5::hdf5_hl-shared;xtensor;fmt::fmt;dl;pugixml::pugixml;PNG::PNG;OpenMP::OpenMP_CXX;MPI::MPI_CXX" +) + +# Load information for each installed configuration. +file(GLOB _cmake_config_files "${CMAKE_CURRENT_LIST_DIR}/OpenMCTargets-*.cmake") +foreach(_cmake_config_file IN LISTS _cmake_config_files) + include("${_cmake_config_file}") +endforeach() +unset(_cmake_config_file) +unset(_cmake_config_files) + +# Cleanup temporary variables. +set(_IMPORT_PREFIX) + +# Loop over all imported files and verify that they actually exist +foreach(_cmake_target IN LISTS _cmake_import_check_targets) + if(CMAKE_VERSION VERSION_LESS "3.28" + OR NOT DEFINED _cmake_import_check_xcframework_for_${_cmake_target} + OR NOT IS_DIRECTORY "${_cmake_import_check_xcframework_for_${_cmake_target}}") + foreach(_cmake_file IN LISTS "_cmake_import_check_files_for_${_cmake_target}") + if(NOT EXISTS "${_cmake_file}") + message(FATAL_ERROR "The imported target \"${_cmake_target}\" references the file + \"${_cmake_file}\" +but this file does not exist. Possible reasons include: +* The file was deleted, renamed, or moved to another location. +* An install or uninstall procedure did not complete successfully. +* The installation package was faulty and contained + \"${CMAKE_CURRENT_LIST_FILE}\" +but not all the files it references. +") + endif() + endforeach() + endif() + unset(_cmake_file) + unset("_cmake_import_check_files_for_${_cmake_target}") +endforeach() +unset(_cmake_target) +unset(_cmake_import_check_targets) + +# This file does not depend on other imported targets which have +# been exported from the same project but in a separate export set. + +# Commands beyond this point should not need to know the version. +set(CMAKE_IMPORT_FILE_VERSION) +cmake_policy(POP) diff --git a/model.xml b/model.xml new file mode 100644 index 00000000000..e735a584cab --- /dev/null +++ b/model.xml @@ -0,0 +1,43 @@ + + + + + + + + + + + + + + + + + + + + + + + fixed source + 2000 + 10 + + + + 0.0 1.25 2.5 3.75 5.0 + 0.0 0.3141592653589793 0.6283185307179586 0.9424777960769379 1.2566370614359172 1.5707963267948966 1.8849555921538759 2.199114857512855 2.5132741228718345 2.827433388230814 3.141592653589793 3.4557519189487724 3.7699111843077517 4.084070449666731 4.39822971502571 4.71238898038469 5.026548245743669 5.340707511102648 5.654866776461628 5.969026041820607 6.283185307179586 + -5.0 -4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0 + 0.0 5.001 0.0 + + + 1 + + + 1 + flux + collision + + + diff --git a/openmc/data/data.py b/openmc/data/data.py index 5ecadd37be0..6e37378bfc4 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -199,6 +199,7 @@ '(n,3np)': (-3, -1), '(n,n2p)': (-2, -2), '(n,npa)': (-5, -3), + "(n,n')": (0, 0), '(n,gamma)': (1, 0), '(n,p)': (0, -1), '(n,d)': (-1, -1), diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 65b59582cf8..7ac5a9dda4d 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -66,6 +66,7 @@ REACTION_MT = {name: mt for mt, name in REACTION_NAME.items()} REACTION_MT['fission'] = 18 +REACTION_MT["(n,n')"] = 4 FISSION_MTS = (18, 19, 20, 21, 38) diff --git a/openmc/deplete/__init__.py b/openmc/deplete/__init__.py index 052e2245963..0e2f11adfc3 100644 --- a/openmc/deplete/__init__.py +++ b/openmc/deplete/__init__.py @@ -5,12 +5,43 @@ A depletion front-end tool. """ +# Energy group structure constants for isomeric branching +# These must be defined before submodule imports to avoid circular import +CCFE_709_N_GROUPS = 709 +UKAEA_1102_N_GROUPS = 1102 + +ISOMERIC_ENERGY_STRUCTURES = { + 'CCFE-709': CCFE_709_N_GROUPS, + 'UKAEA-1102': UKAEA_1102_N_GROUPS +} + + +def get_energy_structure_name(n_groups): + """Get energy structure name from number of groups. + + Parameters + ---------- + n_groups : int + Number of energy groups + + Returns + ------- + str or None + 'CCFE-709', 'UKAEA-1102', or None if unknown + """ + for name, count in ISOMERIC_ENERGY_STRUCTURES.items(): + if n_groups == count: + return name + return None + + from .nuclide import * from .chain import * from .openmc_operator import * from .coupled_operator import * from .independent_operator import * from .microxs import * +from .gendf import * from .reaction_rates import * from .atom_number import * from .stepresult import * diff --git a/openmc/deplete/_matrix_funcs.py b/openmc/deplete/_matrix_funcs.py index c7f7df76fbd..239e3bb0ad6 100644 --- a/openmc/deplete/_matrix_funcs.py +++ b/openmc/deplete/_matrix_funcs.py @@ -1,78 +1,78 @@ """Functions to form the special matrix for depletion""" -def celi_f1(chain, rates, fission_yields=None): - return (5 / 12 * chain.form_matrix(rates[0], fission_yields) - + 1 / 12 * chain.form_matrix(rates[1], fission_yields)) +def celi_f1(chain, rates, fission_yields=None, isomeric_branching=None): + return (5 / 12 * chain.form_matrix(rates[0], fission_yields, isomeric_branching) + + 1 / 12 * chain.form_matrix(rates[1], fission_yields, isomeric_branching)) -def celi_f2(chain, rates, fission_yields=None): - return (1 / 12 * chain.form_matrix(rates[0], fission_yields) - + 5 / 12 * chain.form_matrix(rates[1], fission_yields)) +def celi_f2(chain, rates, fission_yields=None, isomeric_branching=None): + return (1 / 12 * chain.form_matrix(rates[0], fission_yields, isomeric_branching) + + 5 / 12 * chain.form_matrix(rates[1], fission_yields, isomeric_branching)) -def cf4_f1(chain, rates, fission_yields=None): - return 1 / 2 * chain.form_matrix(rates, fission_yields) +def cf4_f1(chain, rates, fission_yields=None, isomeric_branching=None): + return 1 / 2 * chain.form_matrix(rates, fission_yields, isomeric_branching) -def cf4_f2(chain, rates, fission_yields=None): - return (-1 / 2 * chain.form_matrix(rates[0], fission_yields) - + chain.form_matrix(rates[1], fission_yields)) +def cf4_f2(chain, rates, fission_yields=None, isomeric_branching=None): + return (-1 / 2 * chain.form_matrix(rates[0], fission_yields, isomeric_branching) + + chain.form_matrix(rates[1], fission_yields, isomeric_branching)) -def cf4_f3(chain, rates, fission_yields=None): - return (1 / 4 * chain.form_matrix(rates[0], fission_yields) - + 1 / 6 * chain.form_matrix(rates[1], fission_yields) - + 1 / 6 * chain.form_matrix(rates[2], fission_yields) - - 1 / 12 * chain.form_matrix(rates[3], fission_yields)) +def cf4_f3(chain, rates, fission_yields=None, isomeric_branching=None): + return (1 / 4 * chain.form_matrix(rates[0], fission_yields, isomeric_branching) + + 1 / 6 * chain.form_matrix(rates[1], fission_yields, isomeric_branching) + + 1 / 6 * chain.form_matrix(rates[2], fission_yields, isomeric_branching) + - 1 / 12 * chain.form_matrix(rates[3], fission_yields, isomeric_branching)) -def cf4_f4(chain, rates, fission_yields=None): - return (-1 / 12 * chain.form_matrix(rates[0], fission_yields) - + 1 / 6 * chain.form_matrix(rates[1], fission_yields) - + 1 / 6 * chain.form_matrix(rates[2], fission_yields) - + 1 / 4 * chain.form_matrix(rates[3], fission_yields)) +def cf4_f4(chain, rates, fission_yields=None, isomeric_branching=None): + return (-1 / 12 * chain.form_matrix(rates[0], fission_yields, isomeric_branching) + + 1 / 6 * chain.form_matrix(rates[1], fission_yields, isomeric_branching) + + 1 / 6 * chain.form_matrix(rates[2], fission_yields, isomeric_branching) + + 1 / 4 * chain.form_matrix(rates[3], fission_yields, isomeric_branching)) -def rk4_f1(chain, rates, fission_yields=None): - return 1 / 2 * chain.form_matrix(rates, fission_yields) +def rk4_f1(chain, rates, fission_yields=None, isomeric_branching=None): + return 1 / 2 * chain.form_matrix(rates, fission_yields, isomeric_branching) -def rk4_f4(chain, rates, fission_yields=None): - return (1 / 6 * chain.form_matrix(rates[0], fission_yields) - + 1 / 3 * chain.form_matrix(rates[1], fission_yields) - + 1 / 3 * chain.form_matrix(rates[2], fission_yields) - + 1 / 6 * chain.form_matrix(rates[3], fission_yields)) +def rk4_f4(chain, rates, fission_yields=None, isomeric_branching=None): + return (1 / 6 * chain.form_matrix(rates[0], fission_yields, isomeric_branching) + + 1 / 3 * chain.form_matrix(rates[1], fission_yields, isomeric_branching) + + 1 / 3 * chain.form_matrix(rates[2], fission_yields, isomeric_branching) + + 1 / 6 * chain.form_matrix(rates[3], fission_yields, isomeric_branching)) -def leqi_f1(chain, inputs, fission_yields): - f1 = chain.form_matrix(inputs[0], fission_yields) - f2 = chain.form_matrix(inputs[1], fission_yields) +def leqi_f1(chain, inputs, fission_yields=None, isomeric_branching=None): + f1 = chain.form_matrix(inputs[0], fission_yields, isomeric_branching) + f2 = chain.form_matrix(inputs[1], fission_yields, isomeric_branching) dt_l, dt = inputs[2], inputs[3] return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 -def leqi_f2(chain, inputs, fission_yields=None): - f1 = chain.form_matrix(inputs[0], fission_yields) - f2 = chain.form_matrix(inputs[1], fission_yields) +def leqi_f2(chain, inputs, fission_yields=None, isomeric_branching=None): + f1 = chain.form_matrix(inputs[0], fission_yields, isomeric_branching) + f2 = chain.form_matrix(inputs[1], fission_yields, isomeric_branching) dt_l, dt = inputs[2], inputs[3] return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 -def leqi_f3(chain, inputs, fission_yields=None): - f1 = chain.form_matrix(inputs[0], fission_yields) - f2 = chain.form_matrix(inputs[1], fission_yields) - f3 = chain.form_matrix(inputs[2], fission_yields) +def leqi_f3(chain, inputs, fission_yields=None, isomeric_branching=None): + f1 = chain.form_matrix(inputs[0], fission_yields, isomeric_branching) + f2 = chain.form_matrix(inputs[1], fission_yields, isomeric_branching) + f3 = chain.form_matrix(inputs[2], fission_yields, isomeric_branching) dt_l, dt = inputs[3], inputs[4] return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1 + (dt ** 2 + 6 * dt * dt_l + 5 * dt_l ** 2) / (12 * dt_l * (dt + dt_l)) * f2 + dt_l / (12 * (dt + dt_l)) * f3) -def leqi_f4(chain, inputs, fission_yields=None): - f1 = chain.form_matrix(inputs[0], fission_yields) - f2 = chain.form_matrix(inputs[1], fission_yields) - f3 = chain.form_matrix(inputs[2], fission_yields) +def leqi_f4(chain, inputs, fission_yields=None, isomeric_branching=None): + f1 = chain.form_matrix(inputs[0], fission_yields, isomeric_branching) + f2 = chain.form_matrix(inputs[1], fission_yields, isomeric_branching) + f3 = chain.form_matrix(inputs[2], fission_yields, isomeric_branching) dt_l, dt = inputs[3], inputs[4] return (-dt ** 2 / (12 * dt_l * (dt + dt_l)) * f1 + (dt ** 2 + 2 * dt * dt_l + dt_l ** 2) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index ee742f22bc5..6727c02b1c9 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -733,7 +733,7 @@ def _timed_deplete(self, n, rates, dt, i=None, matrix_func=None): start = time.time() results = deplete( self._solver, self.chain, n, rates, dt, i, matrix_func, - self.transfer_rates, self.external_source_rates) + self.transfer_rates, self.external_source_rates, self.operator) return time.time() - start, results @abstractmethod @@ -871,6 +871,13 @@ def integrate( n = self.operator.initial_condition() t, self._i_res = self._get_start_data() + # Notify user if isomeric branching is enabled + if output and comm.rank == 0: + if (hasattr(self.operator, '_isomeric_branching') and + self.operator._isomeric_branching and + any(bool(d) for d in self.operator._isomeric_branching)): + print("[openmc.deplete] Depletion is using Isomeric Branching") + for i, (dt, source_rate) in enumerate(self): if output and comm.rank == 0: print(f"[openmc.deplete] t={t} s, dt={dt} s, source={source_rate}") diff --git a/openmc/deplete/atom_number.py b/openmc/deplete/atom_number.py index b24c048cc92..0c9e0470570 100644 --- a/openmc/deplete/atom_number.py +++ b/openmc/deplete/atom_number.py @@ -251,4 +251,4 @@ def set_density(self, total_density): """ for i, density_slice in enumerate(total_density): - self.set_mat_slice(i, density_slice) + self.set_mat_slice(i, density_slice) \ No newline at end of file diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index a835face72b..3261b740f1c 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -31,6 +31,7 @@ ReactionInfo = namedtuple('ReactionInfo', ('mts', 'secondaries')) REACTIONS = { + "(n,n')": ReactionInfo({4}, ()), '(n,2nd)': ReactionInfo({11}, ('H2',)), '(n,2n)': ReactionInfo(set(chain([16], range(875, 892))), ()), '(n,3n)': ReactionInfo({17}, ()), @@ -120,9 +121,192 @@ __all__ = ["Chain", "REACTIONS"] +def _parse_isomeric_state(nuclide: str) -> tuple: + """Parse nuclide into base name and isomeric level. + + Parameters + ---------- + nuclide : str + Nuclide name (e.g., 'Hf177', 'Hf177_m1', 'Hf177_m2') + + Returns + ------- + tuple + (base_name, isomeric_level) where level is 0 for ground, 1 for m1, etc. + + Examples + -------- + >>> _parse_isomeric_state('Hf177') + ('Hf177', 0) + >>> _parse_isomeric_state('Hf177_m1') + ('Hf177', 1) + >>> _parse_isomeric_state('Hf177_m2') + ('Hf177', 2) + """ + if '_m' in nuclide: + base, suffix = nuclide.rsplit('_m', 1) + try: + return base, int(suffix) + except ValueError: + # Malformed suffix, treat as ground state + return nuclide, 0 + return nuclide, 0 + + +def _load_isomeric_branching(root): + """Load isomeric branching data from XML root. + + Parameters + ---------- + root : xml.etree.ElementTree + Root of chain XML + + Returns + ------- + dict or None + Isomeric branching data or None if not present + """ + from openmc._xml import get_text + import numpy as np + + isomeric_data = {} + + for nuclide_elem in root.findall('nuclide'): + nuc_name = get_text(nuclide_elem, 'name') + if not nuc_name: + continue + + nuc_reactions = {} + + for reaction_elem in nuclide_elem.findall('reaction'): + rx_type = reaction_elem.get('type') + if not rx_type: + continue + + iso_elem = reaction_elem.find('isomeric_yields') + if iso_elem is None: + continue + + iso_type = iso_elem.get('type', 'energy_dependent') + + if iso_type == 'energy_dependent': + # Parse energies + energies_elem = iso_elem.find('energies') + if energies_elem is None or energies_elem.text is None: + continue + + try: + energies = np.array([float(e) for e in energies_elem.text.split()]) + except (ValueError, AttributeError): + continue + + # Parse targets + targets_elem = iso_elem.find('targets') + if targets_elem is None or targets_elem.text is None: + continue + targets = targets_elem.text.split() + + # Parse branching ratios + branching_elem = iso_elem.find('branching_ratios') + if branching_elem is None or branching_elem.text is None: + continue + + lines = [] + for line in branching_elem.text.strip().split('\n'): + line = line.strip() + if line and not line.startswith('