From b39ee9091eb8b9dc52b62a793d14ce30c88e9eea Mon Sep 17 00:00:00 2001 From: Engineer Date: Sun, 23 Aug 2026 23:52:10 +0200 Subject: [PATCH 1/2] 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) From 56a8e5a971493b085668b1d3d0d4bf7cea2d4537 Mon Sep 17 00:00:00 2001 From: Engineer Date: Sun, 30 Aug 2026 16:55:50 +0200 Subject: [PATCH 2/2] fix: use floating-point accumulator in sum_functions for integer grids np.zeros_like(x) inherits the dtype of the union grid, so when the grid is integer-valued the accumulator cannot hold the floating-point results of combined functions such as Polynomial, raising UFuncTypeError. Initialize the accumulator with a float dtype, restored from the prior behavior of sum(f(x) for f in funcs) which promoted to float. Adds a regression test combining an integer-valued tabulated function with a polynomial (test_sum_functions_integer_grid). Co-authored-by: CAOShurong Signed-off-by: Engineer --- openmc/data/function.py | 7 +++++-- tests/unit_tests/test_function.py | 17 +++++++++++++++++ 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index baff0703af4..6597a0627a1 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -48,8 +48,11 @@ def sum_functions(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) + # tabulated range do not contribute to the sum. The accumulator uses + # a floating-point dtype since the combined functions (e.g., + # polynomials) may return values that cannot be represented + # losslessly in the dtype of an integer-valued grid. + y = np.zeros_like(x, dtype=float) for f in funcs: if isinstance(f, Tabulated1D): within = ((x >= f.x[0]) & (x <= f.x[-1])) | \ diff --git a/tests/unit_tests/test_function.py b/tests/unit_tests/test_function.py index 84105869814..6739c6f45e5 100644 --- a/tests/unit_tests/test_function.py +++ b/tests/unit_tests/test_function.py @@ -110,3 +110,20 @@ def test_sum_functions_polynomial(): 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) + + +def test_sum_functions_integer_grid(): + """Integer-valued union grids combine with floating-point functions. + + When the tabulated grid is integer-valued, the resulting union grid has an + integer dtype. Combining such a function with a polynomial must not raise + a casting error and should return a floating-point result. + """ + p = openmc.data.Polynomial((1.0, -0.5)) + f = openmc.data.Tabulated1D([2, 4], [10, 20]) + + s = openmc.data.sum_functions([f, p]) + assert np.array_equal(s.x, np.array([2, 4])) + assert s.y.dtype == np.float64 + assert s.y[0] == pytest.approx(10.0) + assert s.y[1] == pytest.approx(19.0)