Flux Energy Group Conversion using lethargy-weighted redistribution (#3705)

Co-authored-by: GuySten <62616591+GuySten@users.noreply.github.com>

docs/source/pythonapi/mgxs.rst

@@ -11,6 +11,16 @@ Module Variables
 .. autodata:: openmc.mgxs.GROUP_STRUCTURES
     :annotation:
 
+Functions 
++++++++++
+
+.. autosummary::
+    :toctree: generated
+    :nosignatures:
+    :template: myfunction.rst
+
+    openmc.mgxs.convert_flux_groups
+
 Classes
 +++++++
 

openmc/mgxs/__init__.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from openmc.mgxs.groups import EnergyGroups
+from openmc.mgxs.groups import EnergyGroups, convert_flux_groups
 from openmc.mgxs.library import Library
 from openmc.mgxs.mgxs import *
 from openmc.mgxs.mdgxs import *

openmc/mgxs/groups.py

@@ -300,3 +300,139 @@ class EnergyGroups:
         # Assign merged edges to merged groups
         merged_groups.group_edges = list(merged_edges)
         return merged_groups
+
+
+def convert_flux_groups(flux, source_groups, target_groups):
+    """Convert flux spectrum between energy group structures.
+
+    Uses flux-per-unit-lethargy conservation, which assumes constant flux per
+    unit lethargy within each source group and distributes flux to target
+    groups proportionally to their lethargy width.
+
+    .. versionadded:: 0.15.4
+
+    Parameters
+    ----------
+    flux : Iterable of float
+        Flux values for source groups. Length must equal
+        source_groups.num_groups.
+    source_groups : EnergyGroups or str
+        Energy group structure of the input flux with boundaries in [eV].
+        Can be an EnergyGroups instance or the name of a group structure
+        (e.g., 'CCFE-709').
+    target_groups : EnergyGroups or str
+        Target energy group structure with boundaries in [eV]. Can be an
+        EnergyGroups instance or the name of a group structure
+        (e.g., 'UKAEA-1102').
+
+    Returns
+    -------
+    numpy.ndarray
+        Flux values for target groups. Total flux is conserved for
+        overlapping energy regions.
+
+    Raises
+    ------
+    TypeError
+        If source_groups or target_groups is not EnergyGroups or str
+    ValueError
+        If flux length doesn't match source_groups, or flux contains
+        negative, NaN, or infinite values
+
+    See Also
+    --------
+    EnergyGroups : Energy group structure class
+
+    Notes
+    -----
+    The assumption of constant flux per unit lethargy within each source
+    group is physically reasonable for most reactor spectra but is not
+    exact. For best accuracy, use source spectra with sufficiently fine
+    energy resolution.
+
+    Examples
+    --------
+    Convert FNS 709-group flux to UKAEA-1102 structure:
+
+    >>> import numpy as np
+    >>> flux_709 = np.load('tests/fns_flux_709.npy')
+    >>> flux_1102 = openmc.mgxs.convert_flux_groups(flux_709, 'CCFE-709', 'UKAEA-1102')
+
+    Convert using EnergyGroups instances:
+
+    >>> source = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0])
+    >>> target = openmc.mgxs.EnergyGroups([1.0, 5.0, 10.0, 50.0, 100.0])
+    >>> flux_target = openmc.mgxs.convert_flux_groups([1e8, 2e8], source, target)
+
+    References
+    ----------
+    .. [1] J. J. Duderstadt and L. J. Hamilton, "Nuclear Reactor Analysis,"
+       John Wiley & Sons, 1976.
+    .. [2] M. Fleming and J.-Ch. Sublet, "FISPACT-II User Manual,"
+       UKAEA-R(18)001, UK Atomic Energy Authority, 2018. See GRPCONVERT keyword.
+
+    """
+    # Handle string group structure names
+    if isinstance(source_groups, str):
+        source_groups = EnergyGroups(source_groups)
+    if isinstance(target_groups, str):
+        target_groups = EnergyGroups(target_groups)
+
+    # Type validation
+    cv.check_type('source_groups', source_groups, EnergyGroups)
+    cv.check_type('target_groups', target_groups, EnergyGroups)
+
+    # Convert flux to numpy array
+    flux = np.asarray(flux, dtype=np.float64)
+    if flux.ndim != 1:
+        raise ValueError(f'flux must be 1-dimensional, got shape {flux.shape}')
+
+    # Validate flux length matches source groups
+    if len(flux) != source_groups.num_groups:
+        raise ValueError(
+            f'Length of flux ({len(flux)}) must equal number of source '
+            f'groups ({source_groups.num_groups})'
+        )
+
+    # Check for invalid flux values
+    if np.any(np.isnan(flux)):
+        raise ValueError('flux contains NaN values')
+    if np.any(np.isinf(flux)):
+        raise ValueError('flux contains infinite values')
+    if np.any(flux < 0):
+        raise ValueError('flux values must be non-negative')
+
+    # Get energy edges
+    source_edges = source_groups.group_edges
+    target_edges = target_groups.group_edges
+    num_target = target_groups.num_groups
+
+    # Initialize output array
+    flux_target = np.zeros(num_target)
+
+    # Main conversion loop: distribute flux using lethargy weighting
+    for idx_src, flux_src in enumerate(flux):
+        if flux_src == 0:
+            continue
+
+        e_low_src = source_edges[idx_src]
+        e_high_src = source_edges[idx_src + 1]
+        lethargy_src = np.log(e_high_src / e_low_src)
+
+        for idx_tgt in range(num_target):
+            e_low_tgt = target_edges[idx_tgt]
+            e_high_tgt = target_edges[idx_tgt + 1]
+
+            # Skip non-overlapping groups
+            if e_high_tgt <= e_low_src or e_low_tgt >= e_high_src:
+                continue
+
+            # Calculate overlap region
+            e_low_overlap = max(e_low_src, e_low_tgt)
+            e_high_overlap = min(e_high_src, e_high_tgt)
+            lethargy_overlap = np.log(e_high_overlap / e_low_overlap)
+
+            # Distribute flux proportionally to lethargy fraction
+            flux_target[idx_tgt] += flux_src * (lethargy_overlap / lethargy_src)
+
+    return flux_target

tests/unit_tests/test_mgxs_convert_flux.py

@@ -0,0 +1,62 @@
+"""Tests for openmc.mgxs.convert_flux_groups function."""
+
+import numpy as np
+import pytest
+from pytest import approx
+
+import openmc.mgxs
+
+
+def test_coarse_to_fine():
+    """Test coarse to fine conversion with flux conservation."""
+    source = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0])
+    target = openmc.mgxs.EnergyGroups([1.0, 5.0, 10.0, 50.0, 100.0])
+    flux_source = np.array([1e8, 2e8])
+
+    flux_target = openmc.mgxs.convert_flux_groups(flux_source, source, target)
+
+    # Check conservation
+    assert np.sum(flux_target) == approx(np.sum(flux_source))
+    assert len(flux_target) == 4
+    assert np.all(flux_target >= 0)
+
+
+def test_fine_to_coarse():
+    """Test fine to coarse conversion (reverse direction)."""
+    source = openmc.mgxs.EnergyGroups([1.0, 5.0, 10.0, 50.0, 100.0])
+    target = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0])
+    flux_source = np.array([1e7, 2e7, 3e7, 4e7])
+
+    flux_target = openmc.mgxs.convert_flux_groups(flux_source, source, target)
+
+    assert np.sum(flux_target) == approx(np.sum(flux_source))
+    assert len(flux_target) == 2
+
+
+def test_lethargy_distribution():
+    """Test that flux is distributed by lethargy, not linear energy."""
+    # Single group from 1 to 100 eV
+    source = openmc.mgxs.EnergyGroups([1.0, 100.0])
+    # Split into two groups: 1-10 eV and 10-100 eV
+    target = openmc.mgxs.EnergyGroups([1.0, 10.0, 100.0])
+    flux_source = np.array([1e8])
+
+    flux_target = openmc.mgxs.convert_flux_groups(flux_source, source, target)
+
+    # Each target group spans one decade (ln(10) lethargy each)
+    # So flux should be split 50/50 by lethargy
+    assert flux_target[0] == approx(5e7)
+    assert flux_target[1] == approx(5e7)
+
+
+def test_fns_ccfe709_to_ukaea1102():
+    """Test CCFE-709 to UKAEA-1102 conversion with real FNS flux spectrum."""
+    from pathlib import Path
+    flux_file = Path(__file__).parent.parent / 'fns_flux_709.npy'
+    fns_flux_709 = np.load(flux_file)
+
+    flux_1102 = openmc.mgxs.convert_flux_groups(fns_flux_709, 'CCFE-709', 'UKAEA-1102')
+
+    assert len(flux_1102) == 1102
+    assert np.sum(flux_1102) == approx(np.sum(fns_flux_709), rel=1e-10)
+    assert np.all(flux_1102 >= 0)