From b39ee9091eb8b9dc52b62a793d14ce30c88e9eea Mon Sep 17 00:00:00 2001 From: Engineer Date: Sun, 23 Aug 2026 23:52:10 +0200 Subject: [PATCH] fix: make Tabulated1D array evaluation consistent with scalar path Evaluating openmc.data.Tabulated1D on an array containing values outside the tabulated range returned zeros for those points, while scalar evaluation returns the value at the nearest tabulated endpoint. Assign boundary values to out-of-range points in the array evaluation path so that both paths agree. sum_functions is also updated to evaluate each tabulated component only where it is defined, which preserves the behavior of combined functions (e.g., fission energy release components) whose tabulated components cover different incident energy ranges. Fixes: #4041 Signed-off-by: Engineer --- openmc/data/function.py | 24 ++++++- tests/unit_tests/test_function.py | 112 ++++++++++++++++++++++++++++++ 2 files changed, 134 insertions(+), 2 deletions(-) create mode 100644 tests/unit_tests/test_function.py diff --git a/openmc/data/function.py b/openmc/data/function.py index c5914f513d5..baff0703af4 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -46,8 +46,18 @@ def sum_functions(funcs): # Take the union of all energies (sorted) x = reduce(np.union1d, xs) - # Evaluate each function and add together - y = sum(f(x) for f in funcs) + # Evaluate each function and add together. Tabulated functions are + # only evaluated where they are defined; values beyond a function's + # tabulated range do not contribute to the sum. + y = np.zeros_like(x) + for f in funcs: + if isinstance(f, Tabulated1D): + within = ((x >= f.x[0]) & (x <= f.x[-1])) | \ + np.isclose(x, f.x[0], atol=1e-14) | \ + np.isclose(x, f.x[-1], atol=1e-14) + y[within] += f(x[within]) + else: + y += f(x) return Tabulated1D(x, y) else: # If no tabulated functions are present, we need to combine the @@ -111,6 +121,10 @@ class Tabulated1D(Function1D): >>> [f(xi) for xi in numpy.linspace(0, 10, 5)] [4.0, 4.25, 4.5, 4.75, 5.0] + When evaluated outside the tabulated range, the value at the nearest + tabulated endpoint is returned. This holds whether the argument is a + scalar or an array of values. + Parameters ---------- x : Iterable of float @@ -202,6 +216,12 @@ def __call__(self, x): y[contained] = (yi*np.exp(np.log(xk/xi)/np.log(xi1/xi) *np.log(yi1/yi))) + # Assign boundary values to points that lie outside the tabulated + # range so that array evaluation is consistent with the scalar + # interpolation path, which returns the value at the nearest endpoint. + y[idx < 0] = self.y[0] + y[idx > len(self.x) - 2] = self.y[-1] + # In some cases, x values might be outside the tabulated region due only # to precision, so we check if they're close and set them equal if so. y[np.isclose(x, self.x[0], atol=1e-14)] = self.y[0] diff --git a/tests/unit_tests/test_function.py b/tests/unit_tests/test_function.py new file mode 100644 index 00000000000..84105869814 --- /dev/null +++ b/tests/unit_tests/test_function.py @@ -0,0 +1,112 @@ +#!/usr/bin/env python + +import numpy as np +import pytest + +import openmc.data + + +def test_tabulated1d_scalar_array_consistency(): + """Scalar and array evaluation should agree everywhere, including + outside the tabulated domain (issue #4041).""" + f = openmc.data.Tabulated1D([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]) + + # Out-of-range scalars return nearest endpoint value + assert f(-100.0) == 4.0 + assert f(0.999) == 4.0 + assert f(3.001) == 6.0 + assert f(1e300) == 6.0 + + xs = np.array([-100.0, 0.0, 0.999, 1.0, 1.5, 2.0, 2.5, + 3.0, 3.000000001, 100.0]) + assert np.all(f(xs) == np.array([f(x) for x in xs])) + + +@pytest.mark.parametrize('interp', [1, 2, 3, 4, 5]) +def test_tabulated1d_interpolation_schemes(interp): + """All ENDF interpolation schemes give identical scalar/array results.""" + f = openmc.data.Tabulated1D([1.0, 2.0, 3.0], [1.0, 4.0, 9.0], + breakpoints=[3], interpolation=[interp]) + + xs = np.array([-10.0, 0.5, 1.0, 1.1, 1.5, 1.9, 2.0, 2.5, 2.99, 3.0, 50.0]) + expected = np.array([f(x) for x in xs]) + assert np.allclose(f(xs), expected, rtol=1e-14) + + # Values within the tabulated range should not be zero-filled + assert np.all(f(xs[2:-1]) != 0.0) + + if interp == 1: + # Histogram + assert f(1.5) == 1.0 + elif interp == 2: + # Linear-linear + assert f(1.5) == pytest.approx(2.5) + elif interp == 4: + # Log-linear + assert f(1.5) == pytest.approx(2.0) + + +def test_tabulated1d_multiple_regions(): + """Interpolation regions are respected by both evaluation paths.""" + f = openmc.data.Tabulated1D( + [1.0, 2.0, 3.0, 4.0, 5.0], [10.0, 20.0, 30.0, 40.0, 80.0], + breakpoints=[3, 5], interpolation=[2, 1]) + + assert f(1.5) == pytest.approx(15.0) + assert f(3.5) == 30.0 # histogram region + + xs = np.array([0.0, 1.5, 2.5, 3.0, 3.5, 4.9, 5.0, 6.0]) + assert np.all(f(xs) == np.array([f(x) for x in xs])) + + +def test_tabulated1d_endpoints(): + """Exact endpoints and precision-level perturbations return endpoint + values.""" + f = openmc.data.Tabulated1D([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]) + + assert np.array_equal(f(np.array([1.0, 3.0])), np.array([4.0, 6.0])) + assert f(1.0) == 4.0 + assert f(3.0) == 6.0 + + # Slightly beyond the upper endpoint due to floating point precision + assert f(np.array([3.0*(1.0 + 1e-15)]))[0] == 6.0 + assert f(np.array([3.0 + 1e-6]))[0] == 6.0 + + +def test_tabulated1d_multidimensional_input(): + """Array shape is preserved through evaluation.""" + f = openmc.data.Tabulated1D([1.0, 2.0, 3.0], [4.0, 5.0, 6.0]) + + x = np.array([[0.0, 1.5], [2.5, 9.0]]) + y = f(x) + assert y.shape == (2, 2) + assert np.array_equal(y, np.array([[4.0, 4.5], [5.5, 6.0]])) + + +def test_sum_functions_partial_domains(): + """Combining tabulated functions with differing domains leaves zeros + where a component is undefined.""" + f1 = openmc.data.Tabulated1D([1.0, 2.0, 3.0], [10.0, 20.0, 30.0]) + f2 = openmc.data.Tabulated1D([2.0, 3.0, 4.0], [100.0, 200.0, 400.0]) + + s = openmc.data.sum_functions([f1, f2]) + assert isinstance(s, openmc.data.Tabulated1D) + assert np.array_equal(s.x, np.array([1.0, 2.0, 3.0, 4.0])) + + # Where both functions are defined, they add; outside a function's own + # domain its contribution is zero + assert s.y[0] == pytest.approx(10.0) # f2 undefined at 1.0 + assert s.y[1] == pytest.approx(120.0) + assert s.y[2] == pytest.approx(230.0) + assert s.y[3] == pytest.approx(400.0) # f1 undefined at 4.0 + + +def test_sum_functions_polynomial(): + """Polynomial contributions apply over the full union grid.""" + p = openmc.data.Polynomial((1.0, -0.5)) + f = openmc.data.Tabulated1D([2.0, 4.0], [10.0, 20.0]) + + s = openmc.data.sum_functions([f, p]) + assert np.array_equal(s.x, np.array([2.0, 4.0])) + assert s.y[0] == pytest.approx(10.0) + assert s.y[1] == pytest.approx(19.0)