From fba1392a9b54aff64050e9a6ac12f1af372c9eb2 Mon Sep 17 00:00:00 2001 From: Perry Date: Sat, 7 Feb 2026 20:18:13 -0800 Subject: [PATCH 1/7] Rename get_gendfxs_and_flux Make it consistent with get_microxs_and_flux --- openmc/deplete/microxs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 3ab95b730b7..7bbe262329a 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -250,7 +250,7 @@ def get_microxs_and_flux( return fluxes_with_energy, micros -def get_gendf_and_flux( +def get_gendfxs_and_flux( model: openmc.Model, domains: DomainTypes, gendf_library: PathLike | 'openmc.deplete.gendf.GENDFLibrary', From b71d527c03c66e425d914cba301062f5b52a48ac Mon Sep 17 00:00:00 2001 From: Perry Date: Sat, 7 Feb 2026 21:37:15 -0800 Subject: [PATCH 2/7] Limit flux to local materials for isomeric branching in MPI OMC-Deplete --- openmc/deplete/independent_operator.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index af35953e45c..3ad25b5d703 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -185,7 +185,6 @@ def __init__(self, else: self._flux_with_energy.append((flux_item, None)) - self.fluxes = fluxes super().__init__( materials=materials, @@ -197,6 +196,11 @@ def __init__(self, reduce_chain_level=reduce_chain_level, keep_isomeric_siblings=keep_isomeric_siblings) + # Filter _flux_with_energy to local materials for MPI + if len(self._flux_with_energy) != len(self.local_mats): + local_indices = [self._mat_index_map[m] for m in self.local_mats] + self._flux_with_energy = [self._flux_with_energy[i] for i in local_indices] + # Store parameters for isomeric branching setup self._require_isomeric_branching = require_isomeric_branching # require_isomeric_branching, maybe can be wholly removed self._gendf_library = gendf_library From 600aebd0d889add13438f3d71c886af41f5a4e61 Mon Sep 17 00:00:00 2001 From: Perry Date: Sat, 7 Feb 2026 22:20:29 -0800 Subject: [PATCH 3/7] Simplify GENDFLibrary API Remove energy_structure and backend params from GENDFLibrary. As elsewhere, energy is now auto-detected from the GENDF files. GENDF backend selection is automated. Changes to be committed: modified: openmc/deplete/coupled_operator.py modified: openmc/deplete/gendf.py modified: openmc/deplete/helpers.py modified: openmc/deplete/independent_operator.py modified: openmc/deplete/microxs.py modified: openmc/lib/gendf.py modified: tests/unit_tests/test_gendf_elis_mapping.py modified: tests/unit_tests/test_isomeric_branching_weighting.py modified: tools/add_gendf_isomeric_branching_to_chain.py --- openmc/deplete/coupled_operator.py | 10 -- openmc/deplete/gendf.py | 146 ++++-------------- openmc/deplete/helpers.py | 26 ++-- openmc/deplete/independent_operator.py | 31 ++-- openmc/deplete/microxs.py | 8 +- openmc/lib/gendf.py | 11 ++ tests/unit_tests/test_gendf_elis_mapping.py | 13 +- .../test_isomeric_branching_weighting.py | 58 ++++--- .../add_gendf_isomeric_branching_to_chain.py | 16 +- 9 files changed, 98 insertions(+), 221 deletions(-) diff --git a/openmc/deplete/coupled_operator.py b/openmc/deplete/coupled_operator.py index 593597be230..f5a9d9fae49 100644 --- a/openmc/deplete/coupled_operator.py +++ b/openmc/deplete/coupled_operator.py @@ -381,19 +381,9 @@ def _setup_isomeric_branching(self): self._isomeric_branching = None return - # Check if energy structure was detected - energy_structure = getattr(self, '_isomeric_energy_structure', None) - if energy_structure is None: - # Could not determine energy structure from chain - self._isomeric_helper = None - self._isomeric_branching = None - return - - # All conditions met - create the isomeric branching helper self._isomeric_helper = IsomericBranchingHelper( self.chain, self._gendf_library, - energy_structure=energy_structure ) def _calculate_isomeric_branching(self): diff --git a/openmc/deplete/gendf.py b/openmc/deplete/gendf.py index ab7b5369991..ce75e25e2d1 100644 --- a/openmc/deplete/gendf.py +++ b/openmc/deplete/gendf.py @@ -360,12 +360,9 @@ class _PythonGENDFLibrary: ---------- library_path : path-like Path to directory containing GENDF .asc files - energy_structure : str, optional - Name of the energy group structure. Must be one of 'CCFE-709' or - 'UKAEA-1102'. Default is 'UKAEA-1102'. validate_energy_grid : bool, optional If True, validate that each GENDF file's energy grid matches the - specified energy structure. Default is True. + auto-detected energy structure. Default is True. decay_file : path-like, optional Path to ENDF decay library for ELIS-based isomeric state mapping. Can be a directory of decay files or a single concatenated file. @@ -438,17 +435,15 @@ class _PythonGENDFLibrary: def __init__( self, library_path: PathLike, - energy_structure: str = 'UKAEA-1102', validate_energy_grid: bool = True, decay_file: Optional[PathLike] = None, elis_rtol: float = ELIS_RTOL, elis_atol: float = ELIS_ATOL, skip_zero_elis_metastables: bool = True, - mapping_mode: str = 'elis' + mapping_mode: str = 'elis', + _energy_structure: Optional[str] = None ): - # Validate inputs check_type('library_path', library_path, (str, Path)) - check_value('energy_structure', energy_structure, SUPPORTED_GROUP_STRUCTURES) check_type('validate_energy_grid', validate_energy_grid, bool) check_type('elis_rtol', elis_rtol, float) check_type('elis_atol', elis_atol, float) @@ -463,6 +458,11 @@ def __init__( raise ValueError( f"GENDF library path must be a directory: {self.library_path}") + if _energy_structure is not None: + energy_structure = _energy_structure + else: + energy_structure = detect_energy_structure(self.library_path) + self.energy_structure = energy_structure self.energy_bounds = GROUP_STRUCTURES[energy_structure].copy() self.n_groups = len(self.energy_bounds) - 1 @@ -2320,9 +2320,7 @@ def __repr__(self) -> str: def GENDFLibrary( library_path: PathLike, - energy_structure: str = 'UKAEA-1102', validate_energy_grid: bool = True, - backend: str = 'auto', decay_file: Optional[PathLike] = None, elis_rtol: float = ELIS_RTOL, elis_atol: float = ELIS_ATOL, @@ -2331,54 +2329,38 @@ def GENDFLibrary( ): """Create a GENDF cross-section library with automatic backend selection. - This factory function automatically selects the best available backend - for GENDF cross-section loading: - - **C++ backend**: Fast implementation (2-5x speedup) if OpenMC was built - with C++ GENDF support - - **Python backend**: Pure Python fallback implementation + Energy structure is auto-detected from the GENDF files. The backend + (C++ or Python) is selected automatically: C++ if available and no + decay_file is provided, otherwise Python. Parameters ---------- library_path : path-like Path to directory containing GENDF .asc files - energy_structure : str, optional - Name of the energy group structure. Must be one of 'CCFE-709' or - 'UKAEA-1102'. Default is 'UKAEA-1102'. validate_energy_grid : bool, optional If True, validate that each GENDF file's energy grid matches the - specified energy structure. Default is True. Only used with Python + detected energy structure. Default is True. Only used with Python backend. - backend : {'auto', 'cpp', 'python'}, optional - Backend selection: - - 'auto': Automatically use C++ if available, otherwise Python (default) - - 'cpp': Force C++ backend (raises error if not available) - - 'python': Force Python backend decay_file : path-like, optional Path to ENDF decay library for ELIS-based isomeric state mapping. Can be a directory of decay files or a single concatenated file. When provided, enables accurate mapping of GENDF MF=10 metastable products to OpenMC ``_m{n}`` naming based on excitation energy matching. **Highly recommended** for isomeric branching workflows. - Only used with Python backend. elis_rtol : float, optional Relative tolerance for ELIS matching (default: 0.01 = 1%). - Only used with Python backend when decay_file is provided. elis_atol : float, optional Absolute tolerance in eV for ELIS matching (default: 100.0 eV). - Only used with Python backend when decay_file is provided. skip_zero_elis_metastables : bool, optional If True, skip metastable states with ELIS=0 in decay library (likely - data errors). Default is True. Only used with Python backend when - decay_file is provided. + data errors). Default is True. mapping_mode : {'elis', 'lfs_order'}, optional Isomeric state mapping mode: - 'elis' (default): Use excitation energy (ELIS) matching between GENDF MF=10 products and decay library. Most accurate method. - 'lfs_order': Use FISPACT-like positional mapping where the 1st metastable LFS maps to _m1, 2nd to _m2, etc. Useful for validation - testing against FISPACT-II. **Warning**: LFS order may not match - LISO order for some nuclides (e.g., Ag116). - Only used with Python backend. + testing against FISPACT-II. Returns ------- @@ -2387,17 +2369,14 @@ def GENDFLibrary( Raises ------ - ValueError - If requested backend is not available FileNotFoundError If library_path does not exist ValueError - If energy_structure is not supported + If energy structure cannot be detected from GENDF files Examples -------- - >>> # Auto-select fastest backend (without ELIS mapping) - >>> lib = GENDFLibrary('/path/to/JEFF40-GENDF/', 'UKAEA-1102') + >>> lib = GENDFLibrary('/path/to/JEFF40-GENDF/') >>> xs = lib.get_xs('U235', 102, lib.energy_bounds) >>> >>> # With ELIS-based isomeric mapping (recommended for branching) @@ -2405,49 +2384,6 @@ def GENDFLibrary( ... '/path/to/JEFF40-GENDF/', ... decay_file='/path/to/JEFF40-decay/' ... ) - >>> branching = lib.get_branching_ratios('Ir191', 102) - >>> print(branching.products) # ['Ir192', 'Ir192_m1', 'Ir192_m2'] - >>> - >>> # Cross-library usage with relaxed tolerance - >>> lib = GENDFLibrary( - ... '/path/to/TENDL2017-GENDF/', - ... decay_file='/path/to/UKDD12_decay.dat', - ... elis_rtol=0.10 # Allow 10% tolerance for cross-library mismatches - ... ) - >>> - >>> # Force Python backend - >>> lib_py = GENDFLibrary('/path/to/JEFF40-GENDF/', backend='python') - >>> - >>> # Force C++ backend (error if not available) - >>> lib_cpp = GENDFLibrary('/path/to/JEFF40-GENDF/', backend='cpp') - - Notes - ----- - The C++ backend provides significant performance improvements: - - 2-5x faster GENDF file parsing - - Lower memory overhead - - Better caching performance - - Both backends provide identical interfaces and results, so code using - GENDFLibrary does not need to change when switching backends. - - **ELIS-Based Mapping** (Python backend only, when decay_file is provided): - - Instead of assuming LFS order equals LISO order, the library uses excitation - energy (ELIS) matching: - - 1. For each GENDF MF=10 metastable product, calculate ELIS = QM - QI - 2. Look up matching metastable state in decay library by ELIS - 3. Use decay library's LISO value for ``_m{n}`` naming - - This handles cases like Ir191(n,gamma)->Ir192 where GENDF uses LFS=3,15 - but decay library correctly identifies these as LISO=1,2 (m1, m2). - - .. versionadded:: 0.15.3 - C++ backend and automatic selection - - .. versionadded:: 0.15.3 - ELIS-based isomeric state mapping (decay_file, elis_rtol, elis_atol) See Also -------- @@ -2455,48 +2391,19 @@ def GENDFLibrary( lookup_liso : Find LISO for given excitation energy DecayState : Data class for nuclear state information """ - # Validate backend choice - if backend not in ('auto', 'cpp', 'python'): - raise ValueError( - f"Invalid backend '{backend}'. Must be 'auto', 'cpp', or 'python'") - - # Validate mapping_mode - if mapping_mode not in ('elis', 'lfs_order'): + library_path = Path(library_path) + if not library_path.exists(): + raise FileNotFoundError( + f"GENDF library path does not exist: {library_path}") + if not library_path.is_dir(): raise ValueError( - f"Invalid mapping_mode '{mapping_mode}'. Must be 'elis' or 'lfs_order'") - - # decay_file is optional - only needed for isomeric branching (MF=10) - # Library can still be used for available_nuclides() and XS retrieval without it + f"GENDF library path must be a directory: {library_path}") - # Determine which backend to use - use_cpp = False - if backend == 'cpp': - if _CppGENDFLibrary is None: - raise ValueError( - "C++ backend requested but not available. " - "Ensure OpenMC was built with C++ GENDF support.") - use_cpp = True - elif backend == 'auto': - # If decay_file is provided, force Python backend for ELIS mapping - # C++ backend doesn't support ELIS mapping yet - if decay_file is not None: - use_cpp = False - else: - use_cpp = _CppGENDFLibrary is not None - # else backend == 'python', use_cpp remains False + energy_structure = detect_energy_structure(library_path) - # Warn if C++ backend requested but isomeric mapping parameters provided - if use_cpp and decay_file is not None: - warnings.warn( - "C++ backend does not support isomeric state mapping (neither 'elis' " - "nor 'lfs_order' mode). decay_file and mapping_mode parameters will " - "be ignored. Use backend='python' to enable isomeric state mapping.", - UserWarning - ) + use_cpp = _CppGENDFLibrary is not None and decay_file is None - # Create appropriate backend if use_cpp: - # Use C++ backend energy_bounds = GROUP_STRUCTURES[energy_structure] return _CppGENDFLibrary( str(library_path), @@ -2504,16 +2411,15 @@ def GENDFLibrary( energy_structure ) else: - # Use Python backend return _PythonGENDFLibrary( library_path, - energy_structure, validate_energy_grid, decay_file, elis_rtol, elis_atol, skip_zero_elis_metastables, - mapping_mode + mapping_mode, + _energy_structure=energy_structure ) diff --git a/openmc/deplete/helpers.py b/openmc/deplete/helpers.py index 3fd49054a75..c1aa8b28685 100644 --- a/openmc/deplete/helpers.py +++ b/openmc/deplete/helpers.py @@ -14,10 +14,9 @@ from numpy import dot, zeros, newaxis, asarray import numpy as np -from openmc.mgxs import GROUP_STRUCTURES from openmc.mpi import comm -from openmc.checkvalue import check_type, check_greater_than, check_value +from openmc.checkvalue import check_type, check_greater_than from openmc.data import JOULE_PER_EV, REACTION_MT from openmc.exceptions import OpenMCError from openmc.deplete.gendf import REACTION_TO_MT @@ -1177,9 +1176,7 @@ class IsomericBranchingHelper: Chain containing isomeric branching data gendf_library : GENDFLibrary GENDF library for on-the-fly multigroup cross-section lookup. - Required for σ×φ-weighted isomeric branching. - energy_structure : {'CCFE-709', 'UKAEA-1102'}, optional - Required energy group structure. Default is 'UKAEA-1102'. + Energy structure is derived from the library. Attributes ---------- @@ -1188,7 +1185,7 @@ class IsomericBranchingHelper: isomeric_data : dict or None Isomeric branching data from the chain energy_structure : str - Name of the required energy group structure + Name of the energy group structure (from gendf_library) expected_energies : numpy.ndarray Expected energy bin boundaries in eV n_groups : int @@ -1201,9 +1198,8 @@ def __init__( self, chain: 'Chain', gendf_library, - energy_structure: str = 'UKAEA-1102' ) -> None: - """Initialize with a depletion chain and energy structure. + """Initialize with a depletion chain and GENDF library. Parameters ---------- @@ -1211,26 +1207,22 @@ def __init__( Chain containing isomeric branching data gendf_library : GENDFLibrary GENDF library for on-the-fly multigroup cross-section lookup. - energy_structure : {'CCFE-709', 'UKAEA-1102'}, optional - Energy group structure. Default is 'UKAEA-1102'. + Energy structure is read from the library. Raises ------ ValueError - If energy_structure is not supported + If gendf_library is None """ - check_value('energy_structure', energy_structure, - ['CCFE-709', 'UKAEA-1102']) - if gendf_library is None: raise ValueError("gendf_library is required for IsomericBranchingHelper") self.chain: 'Chain' = chain self.isomeric_data: Optional[Dict] = chain.isomeric_branching - self.energy_structure: str = energy_structure - self.expected_energies: np.ndarray = GROUP_STRUCTURES[energy_structure].copy() - self.n_groups: int = len(self.expected_energies) - 1 self.gendf_library = gendf_library + self.energy_structure: str = gendf_library.energy_structure + self.expected_energies: np.ndarray = gendf_library.energy_bounds.copy() + self.n_groups: int = len(self.expected_energies) - 1 def weighted_branching_ratios( self, diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index 3ad25b5d703..27c732c575f 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -23,7 +23,6 @@ from .results import Results from .helpers import (ChainFissionHelper, ConstantFissionYieldHelper, SourceRateHelper, IsomericBranchingHelper) -from . import get_energy_structure_name class IndependentOperator(OpenMCOperator): @@ -378,6 +377,13 @@ def _handle_issue(message): ) return True # Signal to return early + if self._gendf_library is None: + if _handle_issue( + "Isomeric branching data is present in chain but no GENDF " + "library was provided." + ): + return + # Check flux spectra and energy bins if len(self.fluxes) == 0 or self._energy_bins is None or len(self._energy_bins) == 0: if _handle_issue( @@ -387,27 +393,9 @@ def _handle_issue(message): ): return - # Detect energy structure from first material with flux data - energy_structure = None - for flux_spectrum, energy in self._flux_with_energy: - if energy is not None and isinstance(flux_spectrum, np.ndarray): - n_groups = len(flux_spectrum) - energy_structure = get_energy_structure_name(n_groups) - if energy_structure is not None: - break - - if energy_structure is None: - if _handle_issue( - "Could not determine energy structure from flux spectra. " - "Isomeric branching requires CCFE-709 (709 groups) or " - "UKAEA-1102 (1102 groups) energy structure." - ): - return - helper = IsomericBranchingHelper( self.chain, self._gendf_library, - energy_structure=energy_structure, ) self._isomeric_branching = [] @@ -417,8 +405,9 @@ def _handle_issue(message): if energy is None or not isinstance(flux_spectrum, np.ndarray): raise RuntimeError( f"Material {i} is missing flux spectrum or energy information. " - f"All materials must have flux spectra with {energy_structure} " - f"energy structure when using energy-dependent isomeric branching." + f"All materials must have flux spectra with " + f"{helper.energy_structure} energy structure when using " + f"energy-dependent isomeric branching." ) # Calculate σ×φ-weighted branching diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index 7bbe262329a..d3ee64d701c 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -25,7 +25,7 @@ from ..utility_funcs import h5py_file_or_group import openmc.lib from openmc.mpi import comm -from .gendf import GENDFLibrary, _PythonGENDFLibrary, _CppGENDFLibrary, detect_energy_structure +from .gendf import GENDFLibrary, _PythonGENDFLibrary, _CppGENDFLibrary # Build tuple of valid GENDF backend types for isinstance checks _GENDF_TYPES = (_PythonGENDFLibrary,) @@ -319,8 +319,7 @@ def get_gendfxs_and_flux( """ # Handle GENDF library input if isinstance(gendf_library, (str, Path)): - detected = detect_energy_structure(gendf_library) - gendf_library = GENDFLibrary(gendf_library, energy_structure=detected) + gendf_library = GENDFLibrary(gendf_library) elif not isinstance(gendf_library, _GENDF_TYPES): raise TypeError( f"gendf_library must be a path or GENDFLibrary instance, " @@ -625,8 +624,7 @@ def from_multigroup_flux_with_gendf( """ # Handle GENDF library input if isinstance(gendf_library, (str, Path)): - detected = detect_energy_structure(gendf_library) - gendf_library = GENDFLibrary(gendf_library, energy_structure=detected) + gendf_library = GENDFLibrary(gendf_library) elif not isinstance(gendf_library, _GENDF_TYPES): raise TypeError( f"gendf_library must be a path or GENDFLibrary instance, " diff --git a/openmc/lib/gendf.py b/openmc/lib/gendf.py index d135c88c4b6..5248045e204 100644 --- a/openmc/lib/gendf.py +++ b/openmc/lib/gendf.py @@ -154,6 +154,7 @@ def __init__(self, library_path: str, energy_bounds: np.ndarray, """ # Convert inputs self._energy_bounds = np.asarray(energy_bounds, dtype=np.float64) + self._energy_structure = energy_structure_name n_bounds = len(self._energy_bounds) # Prepare C arrays @@ -227,6 +228,16 @@ def energy_bins(self) -> np.ndarray: """ return self.energy_bounds + @property + def energy_structure(self) -> Optional[str]: + """Name of the energy group structure (e.g., 'CCFE-709', 'UKAEA-1102').""" + return self._energy_structure + + @property + def decay_lookup(self): + """C++ backend does not support decay file lookup.""" + return None + def has_nuclide(self, nuclide: str) -> bool: """ Check if nuclide is available in library. diff --git a/tests/unit_tests/test_gendf_elis_mapping.py b/tests/unit_tests/test_gendf_elis_mapping.py index a586f41ca21..4e35afc0444 100644 --- a/tests/unit_tests/test_gendf_elis_mapping.py +++ b/tests/unit_tests/test_gendf_elis_mapping.py @@ -317,9 +317,7 @@ def test_gendf_library_with_elis_mapping(jeff33_gendf_path, jeff33_decay_path): # Create library with ELIS mapping lib = GENDFLibrary( jeff33_gendf_path, - energy_structure='CCFE-709', decay_file=jeff33_decay_path, - backend='python' ) # Verify decay_lookup is populated @@ -331,13 +329,7 @@ def test_library_without_decay_file(jeff33_gendf_path): """Test that GENDFLibrary works without decay_file but has no decay_lookup.""" from openmc.deplete.gendf import GENDFLibrary - # decay_file is optional - library works without it but can't do isomeric mapping - lib = GENDFLibrary( - jeff33_gendf_path, - energy_structure='CCFE-709', - backend='python' - # decay_file not provided -> decay_lookup will be None - ) + lib = GENDFLibrary(jeff33_gendf_path) # Without decay_file, decay_lookup is None assert lib.decay_lookup is None @@ -347,12 +339,9 @@ def test_library_with_decay_file_works(jeff33_gendf_path, jeff33_decay_path): """Test that GENDFLibrary works with decay_file.""" from openmc.deplete.gendf import GENDFLibrary - # Create library with decay_file (required) lib = GENDFLibrary( jeff33_gendf_path, - energy_structure='CCFE-709', decay_file=jeff33_decay_path, - backend='python' ) # decay_lookup should be populated diff --git a/tests/unit_tests/test_isomeric_branching_weighting.py b/tests/unit_tests/test_isomeric_branching_weighting.py index a144bae230f..829842d9332 100644 --- a/tests/unit_tests/test_isomeric_branching_weighting.py +++ b/tests/unit_tests/test_isomeric_branching_weighting.py @@ -9,6 +9,7 @@ import pytest from unittest.mock import Mock, MagicMock +from openmc.mgxs import GROUP_STRUCTURES from openmc.deplete.helpers import IsomericBranchingHelper from openmc.deplete.chain import Chain @@ -34,7 +35,8 @@ def create_mock_chain_with_isomeric_data(): return chain -def create_mock_gendf_library(n_groups, energy_bins, in_range_mask=None): +def create_mock_gendf_library(n_groups, energy_bins, in_range_mask=None, + energy_structure='CCFE-709'): """Create a mock GENDF library for testing. Parameters @@ -46,8 +48,13 @@ def create_mock_gendf_library(n_groups, energy_bins, in_range_mask=None): in_range_mask : numpy.ndarray, optional Boolean mask for groups with non-zero XS. If None, uses isomeric data energy range (1e5 to 1e7 eV). + energy_structure : str, optional + Energy structure name. Default is 'CCFE-709'. """ mock_gendf = Mock() + mock_gendf.energy_structure = energy_structure + mock_gendf.energy_bounds = energy_bins.copy() + mock_gendf.n_groups = n_groups if in_range_mask is None: # Isomeric data energy range is 1e5 to 1e7 eV @@ -71,7 +78,9 @@ def test_no_extrapolation_below_threshold(): """Verify BR array is zero below isomeric data range.""" chain = create_mock_chain_with_isomeric_data() mock_gendf = Mock() - helper = IsomericBranchingHelper(chain, mock_gendf, energy_structure='CCFE-709') + mock_gendf.energy_structure = 'CCFE-709' + mock_gendf.energy_bounds = GROUP_STRUCTURES['CCFE-709'].copy() + helper = IsomericBranchingHelper(chain, mock_gendf) # Create arrays simulating energy ranges n_groups = 10 @@ -95,7 +104,9 @@ def test_no_extrapolation_above_range(): """Verify BR array is zero above isomeric data range.""" chain = create_mock_chain_with_isomeric_data() mock_gendf = Mock() - helper = IsomericBranchingHelper(chain, mock_gendf, energy_structure='CCFE-709') + mock_gendf.energy_structure = 'CCFE-709' + mock_gendf.energy_bounds = GROUP_STRUCTURES['CCFE-709'].copy() + helper = IsomericBranchingHelper(chain, mock_gendf) n_groups = 10 ratios = np.array([0.9, 0.8, 0.7]) @@ -116,7 +127,9 @@ def test_in_range_values_preserved(): """Verify in-range BR values are correctly mapped.""" chain = create_mock_chain_with_isomeric_data() mock_gendf = Mock() - helper = IsomericBranchingHelper(chain, mock_gendf, energy_structure='CCFE-709') + mock_gendf.energy_structure = 'CCFE-709' + mock_gendf.energy_bounds = GROUP_STRUCTURES['CCFE-709'].copy() + helper = IsomericBranchingHelper(chain, mock_gendf) n_groups = 5 ratios = np.array([0.9, 0.8, 0.7]) @@ -143,14 +156,13 @@ def test_sigma_phi_weighting_basic(): chain = create_mock_chain_with_isomeric_data() # Use CCFE-709 structure (709 groups) - from openmc.mgxs import GROUP_STRUCTURES energy_bins = GROUP_STRUCTURES['CCFE-709'] n_groups = 709 # Create mock GENDF library mock_gendf = create_mock_gendf_library(n_groups, energy_bins) - helper = IsomericBranchingHelper(chain, mock_gendf, energy_structure='CCFE-709') + helper = IsomericBranchingHelper(chain, mock_gendf) # Create flux spectrum - uniform flux_spectrum = np.ones(n_groups) @@ -172,7 +184,6 @@ def test_skip_when_nuclide_not_in_gendf(): """Verify empty dict returned when nuclide not in GENDF library.""" chain = create_mock_chain_with_isomeric_data() - from openmc.mgxs import GROUP_STRUCTURES energy_bins = GROUP_STRUCTURES['CCFE-709'] n_groups = 709 @@ -180,9 +191,11 @@ def test_skip_when_nuclide_not_in_gendf(): # GENDF library that raises KeyError for Ag109 mock_gendf = Mock() + mock_gendf.energy_structure = 'CCFE-709' + mock_gendf.energy_bounds = energy_bins.copy() mock_gendf.get_xs = Mock(side_effect=KeyError("Ag109 not found")) - helper = IsomericBranchingHelper(chain, mock_gendf, energy_structure='CCFE-709') + helper = IsomericBranchingHelper(chain, mock_gendf) result = helper.weighted_branching_ratios(flux_spectrum, energy_bins) @@ -196,7 +209,7 @@ def test_gendf_library_required(): # gendf_library=None should raise ValueError at construction with pytest.raises(ValueError, match="gendf_library is required"): - IsomericBranchingHelper(chain, None, energy_structure='CCFE-709') + IsomericBranchingHelper(chain, None) def test_with_gendf_succeeds(): @@ -207,14 +220,13 @@ def test_with_gendf_succeeds(): """ chain = create_mock_chain_with_isomeric_data() - from openmc.mgxs import GROUP_STRUCTURES energy_bins = GROUP_STRUCTURES['CCFE-709'] n_groups = 709 # Create mock GENDF library mock_gendf = create_mock_gendf_library(n_groups, energy_bins) - helper = IsomericBranchingHelper(chain, mock_gendf, energy_structure='CCFE-709') + helper = IsomericBranchingHelper(chain, mock_gendf) flux_spectrum = np.ones(n_groups) @@ -241,7 +253,6 @@ def test_error_on_group_mismatch(): """Verify ValueError raised when GENDF groups don't match flux.""" chain = create_mock_chain_with_isomeric_data() - from openmc.mgxs import GROUP_STRUCTURES energy_bins = GROUP_STRUCTURES['CCFE-709'] n_groups = 709 @@ -249,11 +260,11 @@ def test_error_on_group_mismatch(): # GENDF library that returns wrong number of groups mock_gendf = Mock() + mock_gendf.energy_structure = 'CCFE-709' + mock_gendf.energy_bounds = energy_bins.copy() mock_gendf.get_xs = Mock(return_value=np.ones(500)) # Wrong number! - helper = IsomericBranchingHelper( - chain, energy_structure='CCFE-709', gendf_library=mock_gendf - ) + helper = IsomericBranchingHelper(chain, mock_gendf) with pytest.raises(ValueError, match="group structure mismatch"): helper.weighted_branching_ratios(flux_spectrum, energy_bins) @@ -263,7 +274,6 @@ def test_zero_weight_sum_returns_empty(): """Verify empty dict when all reaction rate is below threshold.""" chain = create_mock_chain_with_isomeric_data() - from openmc.mgxs import GROUP_STRUCTURES energy_bins = GROUP_STRUCTURES['CCFE-709'] n_groups = 709 @@ -276,11 +286,11 @@ def test_zero_weight_sum_returns_empty(): gendf_xs[650:] = 1.0 # Only fast XS mock_gendf = Mock() + mock_gendf.energy_structure = 'CCFE-709' + mock_gendf.energy_bounds = energy_bins.copy() mock_gendf.get_xs = Mock(return_value=gendf_xs) - helper = IsomericBranchingHelper( - chain, energy_structure='CCFE-709', gendf_library=mock_gendf - ) + helper = IsomericBranchingHelper(chain, mock_gendf) result = helper.weighted_branching_ratios(flux_spectrum, energy_bins) @@ -338,9 +348,11 @@ def test_coupled_operator_isomeric_enabled_with_gendf(): # Create minimal mock to test the methods mock_coupled = Mock(spec=CoupledOperator) - mock_coupled._gendf_library = Mock() # Has GENDF library + mock_gendf = Mock() + mock_gendf.energy_structure = 'CCFE-709' + mock_gendf.energy_bounds = GROUP_STRUCTURES['CCFE-709'].copy() + mock_coupled._gendf_library = mock_gendf mock_coupled._rate_helper = Mock(spec=DirectWithFluxHelper) # Correct helper type - mock_coupled._isomeric_energy_structure = 'CCFE-709' mock_coupled.chain = create_mock_chain_with_isomeric_data() # Should enable isomeric branching @@ -369,9 +381,7 @@ def test_coupled_operator_calculate_isomeric_branching(): # Create actual helper (not mock) chain = create_mock_chain_with_isomeric_data() - mock_coupled._isomeric_helper = IsomericBranchingHelper( - chain, energy_structure='CCFE-709', gendf_library=mock_gendf - ) + mock_coupled._isomeric_helper = IsomericBranchingHelper(chain, mock_gendf) # Mock rate helper with flux spectrum mock_rate_helper = Mock(spec=DirectWithFluxHelper) diff --git a/tools/add_gendf_isomeric_branching_to_chain.py b/tools/add_gendf_isomeric_branching_to_chain.py index bad0c6275de..b2138ac32b0 100644 --- a/tools/add_gendf_isomeric_branching_to_chain.py +++ b/tools/add_gendf_isomeric_branching_to_chain.py @@ -31,7 +31,7 @@ from openmc.deplete import Chain from openmc.deplete.gendf import ( - GENDFLibrary, detect_energy_structure, REACTION_TO_MT + GENDFLibrary, REACTION_TO_MT ) @@ -1543,7 +1543,6 @@ def main(endf_gxs_dir, base_chain_file, output_chain_file, mt_list=None, verbose=True, isomer_mapping_log_file=None, renormalization_log_file=None, - energy_structure='auto', prune_nn_prime_self_loops=False, suppress_single_target_yields=False): """ @@ -1597,17 +1596,10 @@ def main(endf_gxs_dir, base_chain_file, output_chain_file, print(f" Loaded {len(chain.nuclides)} nuclides") # Step 2: Detect energy structure - print("\nStep 2: Detecting energy structure...") - if energy_structure == 'auto': - energy_structure = detect_energy_structure(endf_gxs_dir) - print(f" Energy structure: {energy_structure}") - - # Step 3: Load GENDF library and extract branching - print("\nStep 3: Loading GENDF library...") + # Step 2-3: Load GENDF library (energy structure auto-detected) + print("\nStep 2: Loading GENDF library...") lib_kwargs = { - 'energy_structure': energy_structure, 'validate_energy_grid': False, - 'use_fast_parser': True, 'skip_zero_elis_metastables': skip_zero_elis_metastables, 'decay_file': decay_file, 'elis_rtol': elis_rtol, @@ -1618,6 +1610,7 @@ def main(endf_gxs_dir, base_chain_file, output_chain_file, print(f" ELIS tolerance: rtol={elis_rtol} ({elis_rtol*100:.0f}%), atol={elis_atol} eV") lib = GENDFLibrary(endf_gxs_dir, **lib_kwargs) + print(f" Energy structure: {lib.energy_structure}") print(f" Available nuclides: {len(lib.available_nuclides())}") print("\nStep 4: Extracting MF=10 branching data...") @@ -1790,7 +1783,6 @@ def main(endf_gxs_dir, base_chain_file, output_chain_file, elis_atol=args.atol, verbose=verbose, isomer_mapping_log_file=log_file, - energy_structure='auto', prune_nn_prime_self_loops=args.prune_nn_prime_self_loops, suppress_single_target_yields=args.suppress_single_target_yields ) From e981c02f02b72bf74743db096314b2a791a93faf Mon Sep 17 00:00:00 2001 From: Perry Date: Sun, 8 Feb 2026 17:56:39 -0800 Subject: [PATCH 4/7] Reduce MPI per-rank memory in IndependentOperator, isomeric version MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit In __init__, filter fluxes, cross_sections, and _flux_with_energy to local materials only and remap _mat_index_map to local sequential indices. This reduces per-rank memory from O(all_materials) to O(local_materials). In _update_materials, replace the broadcast loop with local-only negative-density zeroing. Unlike CoupledOperator, IndependentOperator has no C library to synchronize material compositions with — the only required action is zeroing negative densities produced by CRAM. Works in and around the isomeric branching implementation. --- openmc/deplete/independent_operator.py | 55 +++++++++----------------- 1 file changed, 19 insertions(+), 36 deletions(-) diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py index 27c732c575f..45c7e14b386 100644 --- a/openmc/deplete/independent_operator.py +++ b/openmc/deplete/independent_operator.py @@ -195,13 +195,18 @@ def __init__(self, reduce_chain_level=reduce_chain_level, keep_isomeric_siblings=keep_isomeric_siblings) - # Filter _flux_with_energy to local materials for MPI - if len(self._flux_with_energy) != len(self.local_mats): + # Filter fluxes, cross sections, and flux-with-energy to local + # materials only (MPI RAM savings: O(all_mats) -> O(local_mats)) + if comm.size > 1: local_indices = [self._mat_index_map[m] for m in self.local_mats] + self.fluxes = [self.fluxes[i] for i in local_indices] + self.cross_sections = [self.cross_sections[i] for i in local_indices] self._flux_with_energy = [self._flux_with_energy[i] for i in local_indices] + self._mat_index_map = { + lm: i for i, lm in enumerate(self.local_mats)} # Store parameters for isomeric branching setup - self._require_isomeric_branching = require_isomeric_branching # require_isomeric_branching, maybe can be wholly removed + self._require_isomeric_branching = require_isomeric_branching self._gendf_library = gendf_library # Setup isomeric branching after initialization @@ -584,36 +589,14 @@ def __call__(self, vec, source_rate): return copy.deepcopy(op_result) def _update_materials(self): - """Updates material compositions in OpenMC on all processes.""" - - for rank in range(comm.size): - number_i = comm.bcast(self.number, root=rank) - - for mat in number_i.materials: - nuclides = [] - densities = [] - for nuc in number_i.nuclides: - if nuc in self.nuclides_with_data: - val = 1.0e-24 * number_i.get_atom_density(mat, nuc) - - # If nuclide is zero, do not add to the problem. - if val > 0.0: - if self.round_number: - val_magnitude = np.floor(np.log10(val)) - val_scaled = val / 10**val_magnitude - val_round = round(val_scaled, 8) - - val = val_round * 10**val_magnitude - - nuclides.append(nuc) - densities.append(val) - else: - # Only output warnings if values are significantly - # negative. CRAM does not guarantee positive - # values. - if val < -1.0e-21: - print(f'WARNING: nuclide {nuc} in material' - f'{mat} is negative (density = {val}' - - ' atom/b-cm)') - number_i[mat, nuc] = 0.0 + """Zero out negative nuclide densities on local materials.""" + for mat in self.number.materials: + for nuc in self.number.nuclides: + if nuc in self.nuclides_with_data: + val = 1.0e-24 * self.number.get_atom_density(mat, nuc) + if val < 0.0: + if val < -1.0e-21: + print(f'WARNING: nuclide {nuc} in material ' + f'{mat} is negative (density = {val}' + ' atom/b-cm)') + self.number[mat, nuc] = 0.0 From 629e94823032755d8b6184fb05fe66c44e48f5dd Mon Sep 17 00:00:00 2001 From: Perry Date: Sun, 8 Feb 2026 18:04:03 -0800 Subject: [PATCH 5/7] energy_bounds parameter in _PythonGENDFLibrary.get_xs() now defaults to the library's own energy structure (matching C++ backend API). Add _energy_validated flag to skip redundant np.allclose() on 1103-element arrays after the first successful validation. --- openmc/deplete/gendf.py | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/openmc/deplete/gendf.py b/openmc/deplete/gendf.py index ce75e25e2d1..36902661a28 100644 --- a/openmc/deplete/gendf.py +++ b/openmc/deplete/gendf.py @@ -484,6 +484,9 @@ def __init__( n_nuclides = len(self.decay_lookup) n_states = sum(len(states) for states in self.decay_lookup.values()) + # Skip redundant energy validation after first successful check + self._energy_validated = False + # Cache for loaded ENDF materials # Key: nuclide name (OpenMC format), Value: endf.Material object OR dict self._material_cache = {} @@ -984,7 +987,7 @@ def get_xs( self, nuclide_name: str, mt: int, - energy_bounds: np.ndarray, + energy_bounds: Optional[np.ndarray] = None, strict_alignment: bool = True ) -> np.ndarray: """Get group-averaged cross-section for a nuclide and reaction. @@ -995,8 +998,9 @@ def get_xs( Nuclide name in OpenMC format (e.g., 'U235', 'Ac225') mt : int ENDF MT number for the reaction - energy_bounds : numpy.ndarray - Energy group boundaries in eV. Must match the library's energy structure. + energy_bounds : numpy.ndarray, optional + Energy group boundaries in eV. If None, uses the library's + energy structure. If provided, must match the library's structure. strict_alignment : bool, optional If True (default), raise ValueError when GENDF energy grid for threshold reactions cannot be exactly aligned with library energy @@ -1024,13 +1028,15 @@ def get_xs( only contain data above the threshold. The start and end energies must align exactly with library group boundaries for accurate placement. """ - # Validate energy bounds - # May need to relax GENDF_RTOL_MATCH if too strict for some GENDF files - if not np.allclose(energy_bounds, self.energy_bounds, - rtol=GENDF_RTOL_MATCH, atol=GENDF_ATOL): - raise ValueError( - f"Provided energy bounds do not match library energy structure " - f"'{self.energy_structure}'") + if energy_bounds is None: + energy_bounds = self.energy_bounds + elif not self._energy_validated: + if not np.allclose(energy_bounds, self.energy_bounds, + rtol=GENDF_RTOL_MATCH, atol=GENDF_ATOL): + raise ValueError( + f"Provided energy bounds do not match library energy " + f"structure '{self.energy_structure}'") + self._energy_validated = True # Load material material = self._load_material(nuclide_name) From 523499736f0e9bd96e5daaf4534b62e9d5c92f0d Mon Sep 17 00:00:00 2001 From: Perry Date: Sun, 8 Feb 2026 18:33:29 -0800 Subject: [PATCH 6/7] Add get_all_xs() and _SparseXSTable for batch XS extraction MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The current per-reaction get_xs() loop generates ~35,000 KeyError exceptions per depletion step (most nuclides lack most reactions) and re-loads the GENDF material on each call. get_all_xs() loads the material once and returns only the reactions that exist, eliminating the exception overhead. _SparseXSTable stores only the non-zero (nuclide, reaction) XS as a contiguous (nnz, n_groups) matrix. This enables collapsing all XS in a single BLAS dgemv call instead of N_nuc × N_rxn individual dot products, and reduces memory from ~370 MB dense to ~44 MB sparse. Extract threshold alignment logic into _extract_xs() shared by get_xs() and get_all_xs(). C++ backend wrapper loops get_xs() per MT from a caller-provided list. _build_sparse_xs_table() builder assembles the table from get_all_xs() calls. Not yet wired into the depletion workflow. Changes to be committed: modified: openmc/deplete/gendf.py modified: openmc/deplete/microxs.py modified: openmc/lib/gendf.py --- openmc/deplete/gendf.py | 171 ++++++++++++++++++++++---------------- openmc/deplete/microxs.py | 89 ++++++++++++++++++++ openmc/lib/gendf.py | 31 +++++++ 3 files changed, 220 insertions(+), 71 deletions(-) diff --git a/openmc/deplete/gendf.py b/openmc/deplete/gendf.py index 36902661a28..ffcc2836786 100644 --- a/openmc/deplete/gendf.py +++ b/openmc/deplete/gendf.py @@ -1047,97 +1047,126 @@ def get_xs( f"Reaction MT={mt} not found for {nuclide_name} in GENDF library. " f"Available reactions: {[mt for mf,mt in material.section_data.keys() if mf==3]}") - # Extract cross-section data xs_data = material.section_data[3, mt] + return self._extract_xs(xs_data, nuclide_name, mt, strict_alignment) + def get_all_xs( + self, + nuclide_name: str, + strict_alignment: bool = True + ) -> dict[int, np.ndarray]: + """Get all available cross-sections for a nuclide. + + Loads the material once and extracts all MF=3 reactions in a single + pass. More efficient than calling :meth:`get_xs` per reaction. + + Parameters + ---------- + nuclide_name : str + Nuclide name in OpenMC format (e.g., 'U235', 'Ac225') + strict_alignment : bool, optional + If True (default), raise ValueError for threshold alignment + failures. See :meth:`get_xs` for details. + + Returns + ------- + dict of int to numpy.ndarray + Mapping of MT number to group-averaged cross-section array. + Each array has length n_groups. + """ + material = self._load_material(nuclide_name) + result = {} + for (mf, mt), xs_data in material.section_data.items(): + if mf != 3: + continue + if 'sigma' not in xs_data: + continue + result[mt] = self._extract_xs( + xs_data, nuclide_name, mt, strict_alignment) + return result + + def _extract_xs(self, xs_data, nuclide_name, mt, strict_alignment): + """Extract and align XS from a single MF=3 section. Returns a copy.""" if 'sigma' not in xs_data: raise ValueError( f"No 'sigma' data in MF=3, MT={mt} for {nuclide_name}") sigma = xs_data['sigma'] - - # Handle threshold reactions that may have partial energy coverage gendf_energies = sigma.x gendf_xs = sigma.y - # Check if this is a full or partial energy range + # Full energy range if len(gendf_energies) == len(self.energy_bounds): - # Full energy range - use first n_groups values directly - # Return view instead of copy for performance (values are read-only in usage) - xs_values = gendf_xs[:self.n_groups] + return gendf_xs[:self.n_groups].copy() + + # Partial energy range (threshold reaction) + xs_values = np.zeros(self.n_groups) + + start_energy = gendf_energies[0] + end_energy = gendf_energies[-1] + + # Try exact match for start boundary + start_matches = np.where(np.isclose( + self.energy_bounds, start_energy, + rtol=GENDF_RTOL_MATCH, atol=GENDF_ATOL))[0] + + if len(start_matches) == 1: + start_idx = start_matches[0] + elif len(start_matches) > 1: + warnings.warn( + f"Multiple energy boundary matches for {nuclide_name} MT={mt} " + f"start energy {start_energy:.6e} eV. Using first match.", + UserWarning) + start_idx = start_matches[0] else: - # Partial energy range (threshold reaction) - # Find where the GENDF energies fit in the full energy structure - xs_values = np.zeros(self.n_groups) - - start_energy = gendf_energies[0] - end_energy = gendf_energies[-1] - - # Try exact match for start boundary - start_matches = np.where(np.isclose( - self.energy_bounds, start_energy, - rtol=GENDF_RTOL_MATCH, atol=GENDF_ATOL))[0] - - if len(start_matches) == 1: - start_idx = start_matches[0] - elif len(start_matches) > 1: - # Ambiguous match - should not happen with unique energy bounds + # No exact match found + nearest_idx = np.argmin(np.abs(self.energy_bounds - start_energy)) + nearest_energy = self.energy_bounds[nearest_idx] + relative_diff = abs(start_energy - nearest_energy) / max(start_energy, 1e-10) + + if strict_alignment: + raise ValueError( + f"Cannot align GENDF energy grid for {nuclide_name} MT={mt}. " + f"GENDF starts at {start_energy:.6e} eV, " + f"nearest library boundary is {nearest_energy:.6e} eV " + f"(relative difference: {relative_diff:.2e}). " + f"This may indicate incompatible energy structures. " + f"Set strict_alignment=False to use nearest-group alignment.") + else: warnings.warn( - f"Multiple energy boundary matches for {nuclide_name} MT={mt} " - f"start energy {start_energy:.6e} eV. Using first match.", + f"Energy alignment uncertainty for {nuclide_name} MT={mt}: " + f"GENDF starts at {start_energy:.6e} eV, " + f"using nearest boundary {nearest_energy:.6e} eV " + f"(relative difference: {relative_diff:.2e}). " + f"Cross-section placement may be off by one energy group.", UserWarning) - start_idx = start_matches[0] - else: - # No exact match found - nearest_idx = np.argmin(np.abs(self.energy_bounds - start_energy)) - nearest_energy = self.energy_bounds[nearest_idx] - relative_diff = abs(start_energy - nearest_energy) / max(start_energy, 1e-10) - + start_idx = nearest_idx + + # Calculate number of GENDF groups and end index + n_gendf_groups = len(gendf_energies) - 1 + end_idx = min(start_idx + n_gendf_groups, self.n_groups) + + # Validate end boundary (sanity check for contiguous data) + if end_idx < self.n_groups: + expected_end = self.energy_bounds[end_idx] + end_rtol = GENDF_RTOL_MATCH * 10 # 1e-5 for end boundary + if not np.isclose(expected_end, end_energy, rtol=end_rtol, atol=GENDF_ATOL): + relative_diff_end = abs(end_energy - expected_end) / max(end_energy, 1e-10) if strict_alignment: raise ValueError( - f"Cannot align GENDF energy grid for {nuclide_name} MT={mt}. " - f"GENDF starts at {start_energy:.6e} eV, " - f"nearest library boundary is {nearest_energy:.6e} eV " - f"(relative difference: {relative_diff:.2e}). " - f"This may indicate incompatible energy structures. " - f"Set strict_alignment=False to use nearest-group alignment.") + f"GENDF energy range for {nuclide_name} MT={mt} does not " + f"align with library structure. End energy mismatch: " + f"GENDF {end_energy:.6e} eV vs expected {expected_end:.6e} eV " + f"(relative difference: {relative_diff_end:.2e})") else: warnings.warn( - f"Energy alignment uncertainty for {nuclide_name} MT={mt}: " - f"GENDF starts at {start_energy:.6e} eV, " - f"using nearest boundary {nearest_energy:.6e} eV " - f"(relative difference: {relative_diff:.2e}). " - f"Cross-section placement may be off by one energy group.", + f"End energy mismatch for {nuclide_name} MT={mt}: " + f"GENDF {end_energy:.6e} eV vs expected {expected_end:.6e} eV " + f"(relative difference: {relative_diff_end:.2e})", UserWarning) - start_idx = nearest_idx - - # Calculate number of GENDF groups and end index - n_gendf_groups = len(gendf_energies) - 1 - end_idx = min(start_idx + n_gendf_groups, self.n_groups) - - # Validate end boundary also matches (sanity check for contiguous data) - if end_idx < self.n_groups: - expected_end = self.energy_bounds[end_idx] - end_rtol = GENDF_RTOL_MATCH * 10 # 1e-5 for end boundary sanity check - if not np.isclose(expected_end, end_energy, rtol=end_rtol, atol=GENDF_ATOL): - relative_diff_end = abs(end_energy - expected_end) / max(end_energy, 1e-10) - if strict_alignment: - raise ValueError( - f"GENDF energy range for {nuclide_name} MT={mt} does not " - f"align with library structure. End energy mismatch: " - f"GENDF {end_energy:.6e} eV vs expected {expected_end:.6e} eV " - f"(relative difference: {relative_diff_end:.2e})") - else: - warnings.warn( - f"End energy mismatch for {nuclide_name} MT={mt}: " - f"GENDF {end_energy:.6e} eV vs expected {expected_end:.6e} eV " - f"(relative difference: {relative_diff_end:.2e})", - UserWarning) - - # Fill in the cross-section values for the available energy range - n_values = min(end_idx - start_idx, len(gendf_xs)) - xs_values[start_idx:start_idx + n_values] = gendf_xs[:n_values] + n_values = min(end_idx - start_idx, len(gendf_xs)) + xs_values[start_idx:start_idx + n_values] = gendf_xs[:n_values] return xs_values def has_nuclide(self, nuclide_name: str) -> bool: diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index d3ee64d701c..c5c568db793 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -6,6 +6,7 @@ from __future__ import annotations from collections.abc import Sequence +from dataclasses import dataclass from pathlib import Path import shutil from tempfile import TemporaryDirectory @@ -425,6 +426,94 @@ def get_gendfxs_and_flux( return fluxes_with_energy, micros +@dataclass +class _SparseXSTable: + """Sparse storage of GENDF cross-sections for vectorized flux collapse. + + Stores only non-zero (nuclide, reaction) pairs. xs_matrix has shape + (nnz, n_groups); nuc_indices and rxn_indices map rows to positions + in the dense (n_nuclides, n_reactions) result. + """ + nuclides: list[str] + reactions: list[str] + n_groups: int + xs_matrix: np.ndarray + nuc_indices: np.ndarray + rxn_indices: np.ndarray + + def collapse(self, phi_norm: np.ndarray) -> np.ndarray: + """Collapse to one-group XS. phi_norm must sum to 1.""" + collapsed_sparse = self.xs_matrix @ phi_norm + result = np.zeros((len(self.nuclides), len(self.reactions))) + result[self.nuc_indices, self.rxn_indices] = collapsed_sparse + return result + + +def _build_sparse_xs_table( + gendf_library, + nuclides: list[str], + reactions: list[str], + mts: list[int] +) -> _SparseXSTable: + """Build a sparse XS table from a GENDF library. + + Parameters + ---------- + gendf_library : _PythonGENDFLibrary or _CppGENDFLibrary + GENDF library instance + nuclides : list of str + Nuclide names to include + reactions : list of str + Reaction names (parallel to mts) + mts : list of int + MT numbers corresponding to reactions + + Returns + ------- + _SparseXSTable + Sparse table ready for vectorized collapse + """ + if len(reactions) != len(mts): + raise ValueError( + f"reactions ({len(reactions)}) and mts ({len(mts)}) " + f"must have same length") + + n_groups = gendf_library.n_groups + rows = [] + nuc_idx_list = [] + rxn_idx_list = [] + + mt_to_rxn_idx = {mt: i for i, mt in enumerate(mts)} + is_cpp = (_CppGENDFLibrary is not None + and isinstance(gendf_library, _CppGENDFLibrary)) + + for nuc_idx, nuc in enumerate(nuclides): + if is_cpp: + all_xs = gendf_library.get_all_xs(nuc, mts=mts) + else: + all_xs = gendf_library.get_all_xs(nuc) + for mt, xs_arr in all_xs.items(): + if mt not in mt_to_rxn_idx: + continue + rows.append(xs_arr) + nuc_idx_list.append(nuc_idx) + rxn_idx_list.append(mt_to_rxn_idx[mt]) + + if rows: + xs_matrix = np.vstack(rows) + else: + xs_matrix = np.empty((0, n_groups)) + + return _SparseXSTable( + nuclides=nuclides, + reactions=reactions, + n_groups=n_groups, + xs_matrix=xs_matrix, + nuc_indices=np.array(nuc_idx_list, dtype=np.int32), + rxn_indices=np.array(rxn_idx_list, dtype=np.int32), + ) + + class MicroXS: """Microscopic cross section data for use in transport-independent depletion. diff --git a/openmc/lib/gendf.py b/openmc/lib/gendf.py index 5248045e204..9aa65ff75f5 100644 --- a/openmc/lib/gendf.py +++ b/openmc/lib/gendf.py @@ -11,6 +11,7 @@ from . import _dll from .error import _error_handler +from openmc.exceptions import OpenMCError # ============================================================================= @@ -391,6 +392,36 @@ def get_xs(self, nuclide: str, mt: int, return xs_array + def get_all_xs(self, nuclide: str, + mts: Optional[list[int]] = None) -> dict[int, np.ndarray]: + """Get cross-sections for all requested reactions for a nuclide. + + Parameters + ---------- + nuclide : str + Nuclide name (e.g., 'U235', 'Pu239') + mts : list of int, optional + MT numbers to retrieve. If None, raises ValueError (C++ backend + cannot discover available reactions). + + Returns + ------- + dict of int to numpy.ndarray + Mapping of MT number to cross-section array. MTs not found + in the library are silently skipped. + """ + if mts is None: + raise ValueError( + "C++ backend requires explicit mts list; cannot discover " + "available reactions. Use Python backend or pass mts.") + result = {} + for mt in mts: + try: + result[mt] = self.get_xs(nuclide, mt) + except OpenMCError: + pass + return result + def __repr__(self): return (f"GENDFLibrary(lib_id={self._lib_id}, " f"n_groups={self.n_groups})") From affd8aa54397945131c0b7732d6e4667fd191382 Mon Sep 17 00:00:00 2001 From: Perry Date: Sun, 8 Feb 2026 20:26:20 -0800 Subject: [PATCH 7/7] Use _SparseXSTable for vectorized GENDF flux collapse MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace the per-domain, per-nuclide, per-reaction get_xs() loop with a single _SparseXSTable build followed by vectorized matrix-vector collapse. Exploit domain-independence of GENDF cross-sections: extract XS once, collapse per domain via single BLAS dgemv instead of per-nuclide Python dot products. Eliminates D× redundant XS extraction. Reduces peak memory. In get_gendfxs_and_flux(), the previous code called from_multigroup_flux_with_gendf() D times (once per domain), each performing N_nuc × N_rxn individual get_xs() calls with try/except. GENDF cross-sections are domain-independent, so all D calls extracted identical data. The sparse table is now built once and reused for all domains, reducing D × N_nuc × N_rxn get_xs() calls to N_nuc get_all_xs() calls plus D BLAS dgemv collapses. In from_multigroup_flux_with_gendf(), the nested for-nuc/for-mt/try-except loop is replaced with _build_sparse_xs_table() + table.collapse(), giving single-domain callers the same vectorized path. This new workflow Update MockGENDFLibrary in test_gendf_nuclide_filtering.py to implement n_groups and get_all_xs() required by _build_sparse_xs_table(). Changes to be committed: modified: openmc/deplete/microxs.py modified: tests/unit_tests/test_gendf_nuclide_filtering.py Untracked files: tests/unit_tests/test_sparse_xs_table_validation.py --- openmc/deplete/microxs.py | 60 +++++++++---------- .../test_gendf_nuclide_filtering.py | 9 ++- 2 files changed, 35 insertions(+), 34 deletions(-) diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py index c5c568db793..a0292c44e1e 100644 --- a/openmc/deplete/microxs.py +++ b/openmc/deplete/microxs.py @@ -338,10 +338,12 @@ def get_gendfxs_and_flux( reactions = chain.reactions # Get nuclides from chain, filtered to those with GENDF data + available_nuclides = gendf_library.available_nuclides_set() if not nuclides: - available_nuclides = gendf_library.available_nuclides_set() nuclides = [nuc.name for nuc in chain.nuclides if nuc.name in available_nuclides] + else: + nuclides = [nuc for nuc in nuclides if nuc in available_nuclides] # Set up the flux tallies energy_filter = openmc.EnergyFilter(energies) @@ -409,14 +411,23 @@ def get_gendfxs_and_flux( # Create list where each item corresponds to one domain fluxes = list(flux.squeeze((1, 2))) - # Generate microscopic cross sections using GENDF library - micros = [MicroXS.from_multigroup_flux_with_gendf( - multigroup_flux=flux_i, - gendf_library=gendf_library, - chain_file=chain_file, - nuclides=nuclides, - reactions=reactions - ) for flux_i in fluxes] + # Build sparse XS table once (GENDF XS are domain-independent) + mts = [REACTION_MT[name] for name in reactions] + table = _build_sparse_xs_table(gendf_library, nuclides, reactions, mts) + + # Collapse per domain + micros = [] + for flux_i in fluxes: + flux_arr = np.asarray(flux_i, dtype=float) + flux_sum = flux_arr.sum() + if flux_sum == 0.0: + micros.append(MicroXS( + np.zeros((len(nuclides), len(reactions), 1)), + nuclides, reactions)) + continue + collapsed = table.collapse(flux_arr / flux_sum) + micros.append(MicroXS(collapsed[:, :, np.newaxis], + nuclides, reactions)) # Reset tallies model.tallies = original_tallies @@ -746,34 +757,17 @@ def from_multigroup_flux_with_gendf( reactions = chain.reactions mts = [REACTION_MT[name] for name in reactions] - # Create 3D array for microscopic cross sections - microxs_arr = np.zeros((len(nuclides), len(mts), 1)) - # If flux is zero, safely return zero cross sections - multigroup_flux = np.array(multigroup_flux) + multigroup_flux = np.asarray(multigroup_flux, dtype=float) if (flux_sum := multigroup_flux.sum()) == 0.0: - return cls(microxs_arr, nuclides, reactions) - - # Normalize multigroup flux - multigroup_flux /= flux_sum + return cls(np.zeros((len(nuclides), len(mts), 1)), + nuclides, reactions) - # For each nuclide and reaction, get GENDF XS and compute flux-weighted value - for nuc_index, nuc in enumerate(nuclides): - for mt_index, mt in enumerate(mts): - try: - # Get group-averaged cross-section from GENDF - xs_group = gendf_library.get_xs(nuc, mt, energies) + # Build sparse table and collapse with normalized flux + table = _build_sparse_xs_table(gendf_library, nuclides, reactions, mts) + collapsed = table.collapse(multigroup_flux / flux_sum) - # Collapse cross-section - # Since flux is normalized: <σ> = Σ σ_g * φ_g - microxs_arr[nuc_index, mt_index, 0] = np.sum(xs_group * multigroup_flux) - - except (KeyError, openmc.exceptions.OpenMCError): - # Reaction not available in GENDF for this nuclide - # Leave as zero (e.g., H1 doesn't have MT=16 n,2n) - pass - - return cls(microxs_arr, nuclides, reactions) + return cls(collapsed[:, :, np.newaxis], nuclides, reactions) @classmethod def from_csv(cls, csv_file, **kwargs): diff --git a/tests/unit_tests/test_gendf_nuclide_filtering.py b/tests/unit_tests/test_gendf_nuclide_filtering.py index a5f8e836f1f..1c518a17ec0 100644 --- a/tests/unit_tests/test_gendf_nuclide_filtering.py +++ b/tests/unit_tests/test_gendf_nuclide_filtering.py @@ -28,6 +28,7 @@ def __init__(self, available_nuclides): self.energy_structure = 'CCFE-709' # CCFE-709 has 710 energy boundaries (709 groups) self.energy_bounds = np.logspace(-5, np.log10(2e7), 710) + self.n_groups = 709 def available_nuclides_set(self): return self._available @@ -36,9 +37,15 @@ def get_xs(self, nuclide, mt, energies): """Return mock cross-section data.""" if nuclide not in self._available: raise KeyError(f"Nuclide {nuclide} not in GENDF library") - # Return 709 group values (for CCFE-709) return np.ones(709) * 1.0 # 1 barn for all groups + def get_all_xs(self, nuclide): + """Return all reactions for a nuclide.""" + if nuclide not in self._available: + raise KeyError(f"Nuclide {nuclide} not in GENDF library") + # MT 102 = (n,gamma), MT 18 = fission + return {102: np.ones(709) * 1.0, 18: np.ones(709) * 1.0} + class MockChain: """Mock chain with configurable nuclides."""