diff --git a/docs/source/pythonapi/mgxs.rst b/docs/source/pythonapi/mgxs.rst index 4141aa0a09a..10411be87dc 100644 --- a/docs/source/pythonapi/mgxs.rst +++ b/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 +++++++ diff --git a/openmc/mgxs/__init__.py b/openmc/mgxs/__init__.py index 5adbac9c458..519820ad0dc 100644 --- a/openmc/mgxs/__init__.py +++ b/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) diff --git a/openmc/mgxs/groups.py b/openmc/mgxs/groups.py index dbf9a25ed5c..acc7caf4b9c 100644 --- a/openmc/mgxs/groups.py +++ b/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 diff --git a/tests/unit_tests/fine_group_structures_reference.npz b/tests/unit_tests/fine_group_structures_reference.npz new file mode 100644 index 00000000000..91a1f3584c7 Binary files /dev/null and b/tests/unit_tests/fine_group_structures_reference.npz differ diff --git a/tests/unit_tests/test_mgxs_fine_group_structures.py b/tests/unit_tests/test_mgxs_fine_group_structures.py new file mode 100644 index 00000000000..7f22f06c618 --- /dev/null +++ b/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