Build fine energy group structures

docs/source/pythonapi/mgxs.rst

@@ -20,6 +20,7 @@ Functions
     :template: myfunction.rst
 
     openmc.mgxs.convert_flux_groups
+    openmc.mgxs.build_fine_group_structure
 
 Classes
 +++++++

openmc/mgxs/__init__.py

@@ -1,6 +1,7 @@
 import numpy as np
 
-from openmc.mgxs.groups import EnergyGroups, convert_flux_groups
+from openmc.mgxs.groups import (
+    EnergyGroups, convert_flux_groups, build_fine_group_structure)
 from openmc.mgxs.library import Library
 from openmc.mgxs.mgxs import *
 from openmc.mgxs.mdgxs import *
@@ -24,6 +25,9 @@
   intermediate and thermal reactor applications ([SAR1990]_)
 - activation_ energy group structures "VITAMIN-J-42", "VITAMIN-J-175",
   "TRIPOLI-315", "LLNL-616", "CCFE-709_" and "UKAEA-1102_"
+- multi-group binning fine energy structures, Fusion Optimized Multi-Group "FOMG-16000_",
+  LWR U238 Optimized "VESTA-43000_" and general uniform lethargy "VESTA-100000_"
+  ([MORGAN2013]_, [HAECK2007]_)
 
 .. _CASMO: http://large.stanford.edu/courses/2013/ph241/dalvi1/docs/c5.physor2006.pdf
 .. _SCALE44: https://www-nds.iaea.org/publications/indc/indc-czr-0001.pdf
@@ -39,6 +43,9 @@
 .. _CCFE-709: https://fispact.ukaea.uk/wiki/CCFE-709_group_structure
 .. _UKAEA-1102: https://fispact.ukaea.uk/wiki/UKAEA-1102_group_structure
 .. _ECCO-1968: https://serpent.vtt.fi/mediawiki/index.php/ECCO_1968-group_structure
+.. _FOMG-16000: https://doi.org/10.1016/j.anucene.2012.11.028
+.. _VESTA-43000: https://roma.sckcen.be/ws/portalfiles/portal/4558859/An_Optimum_Approach_to_Monte_Carlo_Burn_Up.pdf
+.. _VESTA-100000: https://roma.sckcen.be/ws/portalfiles/portal/4558859/An_Optimum_Approach_to_Monte_Carlo_Burn_Up.pdf
 .. [SAR1990] Sartori, E., OECD/NEA Data Bank: Standard Energy Group Structures
    of Cross Section Libraries for Reactor Shielding, Reactor Cell and Fusion
    Neutronics Applications: VITAMIN-J, ECCO-33, ECCO-2000 and XMAS JEF/DOC-315
@@ -70,6 +77,13 @@
 .. [KIM2020] Kim, K.S., Ade, B., & Luciano, N.  (2020). Development
    of the MPACT 69-group Library for Magnox Reactor Analysis using VERA.
    Proceedings of International Conference on Physics of Reactors PHYSOR2020.
+.. [HAECK2007] Haeck, W. (2007). An Optimum Approach to Monte Carlo Burnup.
+   PhD thesis, Ghent University (see also Haeck, W., & Verboomen, B.,
+   Nucl. Sci. Eng. 156, 180-196, 2007).
+.. [MORGAN2013] Morgan, L. W. G., Sublet, J. C., Haeck, W. (2013). Optimising
+   the nuclear data energy group structure used for fusion systems.
+   Joint International Conference on Supercomputing in Nuclear Applications +
+   Monte Carlo (SNA+MC 2013), Paris (France).
 """
 
 GROUP_STRUCTURES['CASMO-2'] = np.array([
@@ -1344,3 +1358,11 @@
     1.718869e7, 1.733253e7, 1.747757e7, 1.762383e7, 1.777131e7, 1.792002e7,
     1.806998e7, 1.822119e7, 1.837367e7, 1.852742e7, 1.868246e7, 1.883880e7,
     1.899644e7, 1.915541e7, 1.931570e7, 1.947734e7, 1.964033e7])
+# formulaic fine group energy structures
+GROUP_STRUCTURES['FOMG-16000'] = build_fine_group_structure(
+    [1e-5, 1e-3, 1e0, 2e6, 1.96e7], [1, 999, 14000, 1000],
+    spacing=['linear', 'linear', 'log', 'linear'])
+GROUP_STRUCTURES['VESTA-43000'] = build_fine_group_structure(
+    [1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1e0,  1e1, 1e2,  1e3,   1e4,  1e5,  1e6,  1e7, 2e7],
+    [1000, 1000, 1000, 1000, 1000, 4000, 4000, 10000, 10000, 4000, 4000, 1000, 1000])
+GROUP_STRUCTURES['VESTA-100000'] = build_fine_group_structure([1e-5, 2e7], 100000)

openmc/mgxs/groups.py

@@ -1,6 +1,7 @@
 from collections.abc import Iterable
 import copy
 from numbers import Real
+import operator
 
 import numpy as np
 
@@ -438,3 +439,82 @@ def convert_flux_groups(flux, source_groups, target_groups):
             flux_target[idx_tgt] += flux_src * (lethargy_overlap / lethargy_src)
 
     return flux_target
+
+
+def build_fine_group_structure(macro_bounds, bins_per_macro, spacing='log'):
+    """Generate energy group boundaries from macro-intervals and bins_per_macro.
+
+    The energy range is split into macro-intervals by ``macro_bounds``. Interval
+    ``i`` (from ``macro_bounds[i]`` to ``macro_bounds[i+1]``) is subdivided into
+    ``bins_per_macro[i]`` groups of equal lethargy width (``spacing='log'``) or
+    equal energy width (``spacing='linear'``). Ultra-fine structures for
+    multi-group binning of reaction rates are defined this way, e.g. the
+    VESTA-43000 group structure [HAECK2007]_ and the FOMG-16000 group structure
+    [MORGAN2013]_.
+
+    .. versionadded:: 0.16.1
+
+    Parameters
+    ----------
+    macro_bounds : iterable of float
+        Macro-group boundaries in [eV], non-negative and strictly increasing.
+    bins_per_macro : int or iterable of int
+        Number of fine energy groups in each macro-group.
+        A single integer applies to every macro-group interval.
+    spacing : {'log', 'linear'} or iterable of str, optional
+        Spacing within each macro-group: ``'log'`` for equal lethargy
+        width, ``'linear'`` for equal energy width. A single value applies to
+        every macro-group interval. Defaults to 'log'.
+
+    Returns
+    -------
+    numpy.ndarray
+        Group boundaries in [eV], ascending, with one edge more than the total
+        number of groups.
+
+    Example
+    --------
+    250 groups per decade up to 10 MeV plus 250 groups to 20 MeV (the
+    3250-group structure of [HAECK2007]_):
+
+    >>> VESTA3250 = openmc.mgxs.build_fine_group_structure([1e-5, 1e7, 2e7], [3000, 250])
+
+    """
+    macro_bounds = np.asarray(macro_bounds, dtype=float)
+    if macro_bounds.ndim != 1:
+        raise ValueError(
+            f'macro_bounds must be 1-dimensional, got shape {macro_bounds.shape}')
+    cv.check_length('macro_bounds', macro_bounds, 2)
+    if not np.all(np.isfinite(macro_bounds)):
+        raise ValueError(f'macro_bounds must be finite, got {macro_bounds}')
+    cv.check_greater_than('macro_bounds', macro_bounds[0], 0.0, equality=True)
+    cv.check_increasing('macro_bounds', macro_bounds)
+    num_macro_groups = macro_bounds.size - 1
+
+    # A scalar group count or spacing applies to every macro-group interval
+    if np.ndim(bins_per_macro) == 0:
+        bins_per_macro = [bins_per_macro] * num_macro_groups
+    if np.ndim(spacing) == 0:
+        spacing = [spacing] * num_macro_groups
+    cv.check_length('bins_per_macro', bins_per_macro, num_macro_groups, num_macro_groups)
+    cv.check_length('spacing', spacing, num_macro_groups, num_macro_groups)
+
+    edges = [macro_bounds[:1]]
+    for lo, hi, n, kind in zip(macro_bounds[:-1], macro_bounds[1:], bins_per_macro, spacing):
+        # Plain int, so that n + 1 cannot wrap for fixed-width numpy integers
+        n = operator.index(n)
+        if n < 1:
+            raise ValueError(f'Number of groups must be at least 1, got {n}')
+        cv.check_value('spacing', kind, ('log', 'linear'))
+        if kind == 'log':
+            if lo <= 0.0:
+                raise ValueError("A 'log' macro-interval needs a positive "
+                                 f"lower bound, got {lo} eV")
+            edges.append(np.geomspace(lo, hi, n + 1)[1:])
+        else:
+            edges.append(np.linspace(lo, hi, n + 1)[1:])
+
+    # An interval too narrow for its group count would repeat edges
+    edges = np.concatenate(edges)
+    cv.check_increasing('group boundaries', edges)
+    return edges

tests/unit_tests/test_mgxs_fine_group_structures.py

@@ -0,0 +1,174 @@
+"""Tests for openmc.mgxs.GROUP_STRUCTURES and build_fine_group_structure."""
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+import openmc
+from openmc.mgxs import GROUP_STRUCTURES, build_fine_group_structure
+
+REFERENCE_FILE = Path(__file__).with_name('fine_group_structures_reference.npz')
+
+# VESTA-43000 macro group bounds [eV] and number of fine group bins per macro group
+VESTA43000_MACRO_BOUNDS = [1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1., 10., 100., 1e3, 1e4, 1e5,
+                1e6, 1e7, 2e7]
+VESTA43000_BINS_PER_MACRO = [1000, 1000, 1000, 1000, 1000, 4000, 4000, 10000,
+                    10000, 4000, 4000, 1000, 1000]
+
+
+@pytest.fixture(scope='module')
+def reference():
+    with np.load(REFERENCE_FILE) as data:
+        return {name: data[name] for name in data.files}
+
+
+@pytest.mark.parametrize('name', sorted(GROUP_STRUCTURES))
+def test_group_structure_well_formed(name):
+    edges = GROUP_STRUCTURES[name]
+    assert isinstance(edges, np.ndarray)
+    assert edges.ndim == 1
+    assert edges.dtype == np.float64
+    assert np.all(np.isfinite(edges))
+    assert np.all(np.diff(edges) > 0.0)
+    assert edges[0] >= 0.0
+
+
+@pytest.mark.parametrize('spacing', ['log', 'linear'])
+def test_build_scalar_arguments(spacing):
+    macro_bounds = [1e-5, 1., 1e3, 2e7]
+    expected = build_fine_group_structure(macro_bounds, [10, 10, 10], [spacing] * 3)
+    edges = build_fine_group_structure(macro_bounds, 10, spacing)
+    np.testing.assert_array_equal(edges, expected)
+
+
+def test_build_mixed_spacing():
+    macro_bounds = [0., 1., 1e3, 2e7]
+    bins_per_macro = [4, 30, 7]
+    spacing = ['linear', 'log', 'linear']
+    edges = build_fine_group_structure(macro_bounds, bins_per_macro, spacing)
+    assert edges.size == sum(bins_per_macro) + 1
+
+    # Macro group bounds are reproduced exactly
+    idx = np.cumsum([0] + bins_per_macro)
+    np.testing.assert_array_equal(edges[idx], macro_bounds)
+
+    # Constant lethargy ('log') or energy ('linear') width within a macro group
+    for kind, start, stop in zip(spacing, idx[:-1], idx[1:]):
+        segment = edges[start:stop + 1]
+        if kind == 'log':
+            widths = np.diff(np.log(segment))
+        else:
+            widths = np.diff(segment)
+        np.testing.assert_allclose(widths, widths[0], rtol=1e-9)
+
+
+def test_build_linear_from_zero():
+    edges = build_fine_group_structure([0., 1.], 4, 'linear')
+    np.testing.assert_array_equal(edges, [0., 0.25, 0.5, 0.75, 1.])
+
+
+def test_build_uint8_bins_per_macro():
+    # np.uint8(255) + 1 would wrap to 0 if the count stayed a numpy scalar
+    edges = build_fine_group_structure([1., 10.], np.uint8(255))
+    assert edges.size == 256
+    assert edges[0] == 1.
+    assert edges[-1] == 10.
+
+
+def test_build_array_bins_per_macro_tuple_spacing():
+    edges = build_fine_group_structure(
+        [1., 10., 100.], np.array([3, 4], dtype=np.int64),
+        spacing=('log', 'linear'))
+    assert edges.size == 8
+    np.testing.assert_array_equal(edges[[0, 3, 7]], [1., 10., 100.])
+
+
+@pytest.mark.parametrize('macro_bounds, bins_per_macro, spacing', [
+    pytest.param([1., 10., 100.], [10], 'log', id='num-groups-length'),
+    pytest.param([1., 10., 100.], 10, ['log'], id='spacing-length'),
+    pytest.param([1., 100., 10.], 10, 'log', id='decreasing-bounds'),
+    pytest.param([1., 1., 10.], 10, 'log', id='repeated-bound'),
+    pytest.param([1.], 10, 'log', id='single-bound'),
+    pytest.param([[1., 10.]], 10, 'log', id='2d-bounds'),
+    pytest.param([-1., 1.], 2, 'linear', id='negative-bound'),
+    pytest.param([1., np.inf], 2, 'log', id='infinite-bound'),
+    pytest.param([1., 10.], 0, 'log', id='zero-groups'),
+    pytest.param([1., 10.], 10, 'cubic', id='unknown-spacing'),
+    pytest.param([0., 10.], 10, 'log', id='log-from-zero'),
+    pytest.param([1., 1. + 1e-12], 100000, 'log', id='too-narrow-macrobound'),
+])
+def test_build_invalid(macro_bounds, bins_per_macro, spacing):
+    with pytest.raises(ValueError):
+        build_fine_group_structure(macro_bounds, bins_per_macro, spacing)
+
+
+def test_build_non_integer_bins_per_macro():
+    with pytest.raises(TypeError):
+        build_fine_group_structure([1., 10.], 2.5)
+
+
+@pytest.mark.parametrize('name, rtol', [
+    ('VESTA-43000', 1e-12),
+    # FOMG edges were tabulated to 10 significant figures
+    ('FOMG-16000', 1e-9),
+])
+def test_reference_edges(reference, name, rtol):
+    edges = GROUP_STRUCTURES[name]
+    expected = reference[name]
+    assert edges.size == expected.size
+    assert edges[0] == expected[0]
+    assert edges[-1] == expected[-1]
+    np.testing.assert_allclose(edges, expected, rtol=rtol, atol=0)
+
+
+def test_vesta_43000():
+    edges = GROUP_STRUCTURES['VESTA-43000']
+    assert edges.size == 43001
+
+    # Macrogroup endpoints are exact; lethargy width ln(hi/lo)/n inside
+    idx = np.cumsum([0] + VESTA43000_BINS_PER_MACRO)
+    np.testing.assert_array_equal(edges[idx], VESTA43000_MACRO_BOUNDS)
+    for lo, hi, n, start, stop in zip(VESTA43000_MACRO_BOUNDS[:-1], VESTA43000_MACRO_BOUNDS[1:],
+                                      VESTA43000_BINS_PER_MACRO, idx[:-1], idx[1:]):
+        du = np.diff(np.log(edges[start:stop + 1]))
+        np.testing.assert_allclose(du, np.log(hi / lo) / n, rtol=1e-9)
+
+
+def test_fomg_16000():
+    edges = GROUP_STRUCTURES['FOMG-16000']
+    assert edges.size == 16001
+    assert edges[0] == 1e-5
+    assert edges[1000] == 1.0
+    assert edges[15000] == 2e6
+    assert edges[-1] == 1.96e7
+
+    # 1 meV bins from the second edge up to 1 eV
+    np.testing.assert_allclose(np.diff(edges[1:1001]), 1e-3, rtol=1e-9)
+    # Equal-lethargy bins from 1 eV to 2 MeV
+    np.testing.assert_allclose(np.diff(np.log(edges[1000:15001])),
+                               np.log(2e6) / 14000, rtol=1e-9)
+    # 17.6 keV bins from 2 MeV to 19.6 MeV
+    np.testing.assert_allclose(np.diff(edges[15000:]), 17.6e3, rtol=1e-9)
+
+
+def test_vesta_100000():
+    # Check the definition only
+    edges = GROUP_STRUCTURES['VESTA-100000']
+    assert edges.size == 100001
+    assert edges[0] == 1e-5
+    assert edges[-1] == 2e7
+    np.testing.assert_allclose(np.diff(np.log(edges)),
+                               np.log(2e7 / 1e-5) / 100000, rtol=1e-9)
+
+
+@pytest.mark.parametrize('name, n_groups', [
+    ('FOMG-16000', 16000),
+    ('VESTA-43000', 43000),
+    ('VESTA-100000', 100000),
+])
+def test_name_lookup(name, n_groups):
+    assert openmc.EnergyFilter.from_group_structure(name).num_bins == n_groups
+    assert openmc.mgxs.EnergyGroups(name).num_groups == n_groups
+    f = openmc.ParticleProductionFilter('photon', energies=name)
+    assert f.num_energy_bins == n_groups