From eed9e193200841297c1e526f0d905a1aa6a12102 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 4 Dec 2020 08:23:52 -0800 Subject: [PATCH 1/7] Provide ResultsList.get_times: BOS and EOL points in some units Simple method that just returns the points that define the depletion schedule. Simple unit conversion is performed, defaulting to returning time in days. But supports s (seconds), d (days), h (hours) and min (minutes). Would be useful to add a burnup unit too, but want to handle that separately --- openmc/deplete/results_list.py | 38 +++++++++++++++++++++++++++++++++- 1 file changed, 37 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index f9f50b822..2bb72d6ff 100644 --- a/openmc/deplete/results_list.py +++ b/openmc/deplete/results_list.py @@ -2,7 +2,7 @@ import h5py import numpy as np from .results import Results, VERSION_RESULTS -from openmc.checkvalue import check_filetype_version, check_value +from openmc.checkvalue import check_filetype_version, check_value, check_type __all__ = ["ResultsList"] @@ -192,3 +192,39 @@ class ResultsList(list): for ix, res in enumerate(items): times[ix] = res.proc_time return times + + def get_times(self, time_units="d") -> np.ndarray: + """Return the points in time that define the depletion schedule + + Parameters + ---------- + time_units : {"s", "d", "h", "min"}, optional + Return the vector in these units. Default is to + convert to days + + Returns + ------- + numpy.ndarray + 1-D vector of time points + + """ + check_type("time_units", time_units, str) + + times = np.fromiter( + (r.time[0] for r in self), + dtype=self[0].time.dtype, + count=len(self), + ) + + if time_units == "d": + times /= (60 * 60 * 24) + elif time_units == "h": + times /= (60 * 60) + elif time_units == "min": + times /= 60 + elif time_units != "s": + raise ValueError( + f'Unable to set "time_units" to {time_units} since it is not ' + 'in ("s", "d", "min", "h")' + ) + return times From 08923ce94d4946f01bf425a4216d2884c004b069 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 4 Dec 2020 08:28:25 -0800 Subject: [PATCH 2/7] Provide ResultsList.get_step_where: closest index to a given time Closes #1713 by finding the index closet to a user-requested point in time. Support for various units is provided through ResultsList.get_times (s, d, min, h) and there is support for absolute and relative tolerances. --- openmc/deplete/results_list.py | 65 ++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index 2bb72d6ff..390724fc8 100644 --- a/openmc/deplete/results_list.py +++ b/openmc/deplete/results_list.py @@ -1,3 +1,7 @@ +import numbers +import bisect +import math + import h5py import numpy as np @@ -228,3 +232,64 @@ class ResultsList(list): 'in ("s", "d", "min", "h")' ) return times + + def get_step_where( + self, time, time_units="d", atol=1e-6, rtol=1e-3 + ) -> int: + """Return the index closest to a given point in time + + In the event ``time`` lies exactly between two points, the + lower index will be returned. It is possible that the index + will be at most one past the point in time requested, but only + according to tolerances requested. + + Parameters + ---------- + time : float + Desired point in time + time_units : {"s", "d", "min", "h"}, optional + Units on ``time``. Default: days + atol : float, optional + Absolute tolerance (in ``time_units``) if ``time`` is not + found. A value of ``None`` is converted to infinity. + rtol : float, optional + Relative tolerance if ``time`` is not found. A value of + ``None`` is converted to infinity. + + Returns + ------- + int + + """ + check_type("time", time, numbers.Real) + if atol is not None: + check_type("atol", atol, numbers.Real) + else: + atol = math.inf + if rtol is not None: + check_type("rtol", rtol, numbers.Real) + else: + rtol = math.inf + + times = self.get_times(time_units) + + if times[0] < time < times[-1]: + ix = bisect.bisect_left(times, time) + if ix == times.size: + ix -= 1 + # Bisection will place us either directly on the point + # or one-past the first value less than time + elif time - times[ix - 1] <= times[ix] - time: + ix -= 1 + elif times[0] >= time: + ix = 0 + elif time >= times[-1]: + ix = times.size - 1 + + if math.isclose(time, times[ix], rel_tol=rtol, abs_tol=atol): + return ix + + raise ValueError( + f"A value of {time} {time_units} was not found given absolute and " + f"relative tolerances {atol} and {rtol}." + ) From 833529a880187981e4a3f8b548275eda485e4a71 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 4 Dec 2020 08:28:52 -0800 Subject: [PATCH 3/7] Test ResultsList.get_step_where --- tests/unit_tests/test_deplete_resultslist.py | 63 ++++++++++++++++++++ 1 file changed, 63 insertions(+) diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index 9b40a1ac9..8bfabd395 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -65,3 +65,66 @@ def test_get_eigenvalue(res): np.testing.assert_allclose(t, t_ref) np.testing.assert_allclose(k[:, 0], k_ref) np.testing.assert_allclose(k[:, 1], u_ref) + + +@pytest.mark.parametrize("unit", ("s", "d", "min", "h")) +def test_get_steps(unit): + # Make a ResultsList full of near-empty Result instances + # Just fill out a time schedule + results = openmc.deplete.ResultsList() + # Time in units of unit + times = np.linspace(0, 100, num=5) + if unit == "d": + conversion_to_seconds = 60 * 60 * 24 + elif unit == "h": + conversion_to_seconds = 60 * 60 + elif unit == "min": + conversion_to_seconds = 60 + else: + conversion_to_seconds = 1 + + for ix in range(times.size): + res = openmc.deplete.Results() + res.time = times[ix:ix + 1] * conversion_to_seconds + results.append(res) + + for expected, value in enumerate(times): + actual = results.get_step_where( + value, time_units=unit, atol=0, rtol=0) + assert actual == expected, (value, results[actual].time[0]) + + with pytest.raises(ValueError): + # Emulate a result file with a non-zero initial point in time + # as in starting from a restart + results.get_step_where(times[0] - 1, time_units=unit, atol=0, rtol=0) + + with pytest.raises(ValueError): + results.get_step_where(times[-1] + 1, time_units=unit, atol=0, rtol=0) + + # Grab intermediate points with a small offset + delta = (times[1] - times[0]) + offset = delta * 0.1 + for expected, value in enumerate(times[1:-1], start=1): + # Shoot a little low and a little high + for mult in {1, -1}: + target = value + mult * offset + # Compare using absolute and relative tolerances + actual = results.get_step_where( + target, time_units=unit, atol=offset * 2, rtol=None) + assert actual == expected, ( + target, times[actual], times[expected], offset) + + actual = results.get_step_where( + target, time_units=unit, atol=None, rtol=offset / value) + assert actual == expected, ( + target, times[actual], times[expected], offset) + # Check that the lower index is returned for the exact mid-point + target = value + delta * 0.5 + actual = results.get_step_where( + target, time_units=unit, atol=delta, rtol=delta / value) + assert actual == expected + + # Shoot way over with no tolerance -> just give closest value + actual = results.get_step_where( + times[-1] * 100, time_units=unit, atol=None, rtol=None) + assert actual == times.size - 1 From 3905e7a8a05bc72e211b8f72d345361df38220d8 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 9 Dec 2020 08:36:49 -0800 Subject: [PATCH 4/7] Update tests/unit_tests/test_deplete_resultslist.py Co-authored-by: Paul Romano --- tests/unit_tests/test_deplete_resultslist.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index 8bfabd395..c96108cfd 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -106,7 +106,7 @@ def test_get_steps(unit): offset = delta * 0.1 for expected, value in enumerate(times[1:-1], start=1): # Shoot a little low and a little high - for mult in {1, -1}: + for mult in (1, -1): target = value + mult * offset # Compare using absolute and relative tolerances actual = results.get_step_where( From eb5121b14a8d8fb709c06e3b06324908de1d09b4 Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Wed, 9 Dec 2020 08:41:30 -0800 Subject: [PATCH 5/7] Require atol / rtol to be real in ResultsList.get_step_where Previous implementation converted a value of None to infinity. --- openmc/deplete/results_list.py | 18 +++++++----------- tests/unit_tests/test_deplete_resultslist.py | 7 ++++--- 2 files changed, 11 insertions(+), 14 deletions(-) diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index 390724fc8..02e55f8f8 100644 --- a/openmc/deplete/results_list.py +++ b/openmc/deplete/results_list.py @@ -243,6 +243,9 @@ class ResultsList(list): will be at most one past the point in time requested, but only according to tolerances requested. + Passing ``atol=math.inf`` and ``rtol=math.inf`` will return + the closest index to the requested point. + Parameters ---------- time : float @@ -251,10 +254,9 @@ class ResultsList(list): Units on ``time``. Default: days atol : float, optional Absolute tolerance (in ``time_units``) if ``time`` is not - found. A value of ``None`` is converted to infinity. + found. rtol : float, optional - Relative tolerance if ``time`` is not found. A value of - ``None`` is converted to infinity. + Relative tolerance if ``time`` is not found. Returns ------- @@ -262,14 +264,8 @@ class ResultsList(list): """ check_type("time", time, numbers.Real) - if atol is not None: - check_type("atol", atol, numbers.Real) - else: - atol = math.inf - if rtol is not None: - check_type("rtol", rtol, numbers.Real) - else: - rtol = math.inf + check_type("atol", atol, numbers.Real) + check_type("rtol", rtol, numbers.Real) times = self.get_times(time_units) diff --git a/tests/unit_tests/test_deplete_resultslist.py b/tests/unit_tests/test_deplete_resultslist.py index c96108cfd..392d5d0d5 100644 --- a/tests/unit_tests/test_deplete_resultslist.py +++ b/tests/unit_tests/test_deplete_resultslist.py @@ -1,6 +1,7 @@ """Tests the ResultsList class""" from pathlib import Path +from math import inf import numpy as np import pytest @@ -110,12 +111,12 @@ def test_get_steps(unit): target = value + mult * offset # Compare using absolute and relative tolerances actual = results.get_step_where( - target, time_units=unit, atol=offset * 2, rtol=None) + target, time_units=unit, atol=offset * 2, rtol=inf) assert actual == expected, ( target, times[actual], times[expected], offset) actual = results.get_step_where( - target, time_units=unit, atol=None, rtol=offset / value) + target, time_units=unit, atol=inf, rtol=offset / value) assert actual == expected, ( target, times[actual], times[expected], offset) # Check that the lower index is returned for the exact mid-point @@ -126,5 +127,5 @@ def test_get_steps(unit): # Shoot way over with no tolerance -> just give closest value actual = results.get_step_where( - times[-1] * 100, time_units=unit, atol=None, rtol=None) + times[-1] * 100, time_units=unit, atol=inf, rtol=inf) assert actual == times.size - 1 From dfd4624804cee4e8db140a4ae832f2c4856dabde Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 11 Dec 2020 08:54:50 -0800 Subject: [PATCH 6/7] Remove f-strings from results_list.py --- openmc/deplete/results_list.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index 02e55f8f8..708fb83cd 100644 --- a/openmc/deplete/results_list.py +++ b/openmc/deplete/results_list.py @@ -228,8 +228,8 @@ class ResultsList(list): times /= 60 elif time_units != "s": raise ValueError( - f'Unable to set "time_units" to {time_units} since it is not ' - 'in ("s", "d", "min", "h")' + 'Unable to set "time_units" to {} since it is not ' + 'in ("s", "d", "min", "h")'.format(time_units) ) return times @@ -286,6 +286,7 @@ class ResultsList(list): return ix raise ValueError( - f"A value of {time} {time_units} was not found given absolute and " - f"relative tolerances {atol} and {rtol}." + "A value of {} {} was not found given absolute and " + "relative tolerances {} and {}.".format( + time, time_units, atol, rtol) ) From d72978ddf79c12a7ad6f06230f4746d812c95c5d Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Thu, 24 Dec 2020 07:12:12 -0800 Subject: [PATCH 7/7] Apply suggestions from code review Co-authored-by: Paul Romano --- openmc/deplete/results_list.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/openmc/deplete/results_list.py b/openmc/deplete/results_list.py index 708fb83cd..bd02d576e 100644 --- a/openmc/deplete/results_list.py +++ b/openmc/deplete/results_list.py @@ -200,6 +200,9 @@ class ResultsList(list): def get_times(self, time_units="d") -> np.ndarray: """Return the points in time that define the depletion schedule + + .. versionadded:: 0.12.1 + Parameters ---------- time_units : {"s", "d", "h", "min"}, optional @@ -246,6 +249,9 @@ class ResultsList(list): Passing ``atol=math.inf`` and ``rtol=math.inf`` will return the closest index to the requested point. + + .. versionadded:: 0.12.1 + Parameters ---------- time : float