From 671db024db057c17350666429890a805515effdc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 9 Nov 2017 14:01:26 -0600 Subject: [PATCH] Fix bugs with statepoint loading. Add properties to capi.Tally. --- docs/source/conf.py | 5 +- docs/source/pythonapi/capi.rst | 6 ++ openmc/capi/core.py | 22 ++++++- openmc/capi/tally.py | 117 +++++++++++++++++++++++++-------- src/input_xml.F90 | 2 +- src/state_point.F90 | 12 ++-- src/tallies/tally_header.F90 | 16 +++++ 7 files changed, 143 insertions(+), 37 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index 95930fedb..5f993b1a3 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -26,8 +26,9 @@ except ImportError: MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial', 'numpy.ctypeslib', 'scipy', 'scipy.sparse', 'scipy.interpolate', - 'scipy.integrate', 'scipy.optimize', 'scipy.special', 'h5py', - 'pandas', 'uncertainties', 'openmoc', 'openmc.data.reconstruct'] + 'scipy.integrate', 'scipy.optimize', 'scipy.special', + 'scipy.stats','h5py', 'pandas', 'uncertainties', 'openmoc', + 'openmc.data.reconstruct'] sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES) import numpy as np diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index 457f05320..44a755bbd 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -18,12 +18,18 @@ Functions openmc.capi.find_material openmc.capi.hard_reset openmc.capi.init + openmc.capi.iter_batches openmc.capi.keff openmc.capi.load_nuclide + openmc.capi.next_batch openmc.capi.plot_geometry openmc.capi.reset openmc.capi.run openmc.capi.run_in_memory + openmc.capi.simulation_init + openmc.capi.simulation_finalize + openmc.capi.source_bank + openmc.capi.statepoint_write Classes ------- diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 530ed5063..77c9ca7a4 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -122,7 +122,27 @@ def init(intracomm=None): def iter_batches(): - """Iterator over batches.""" + """Iterator over batches. + + This function returns a generator-iterator that allows Python code to be run + between batches in an OpenMC simulation. It should be used in conjunction + with :func:`openmc.capi.simulation_init` and + :func:`openmc.capi.simulation_finalize`. For example: + + .. code-block:: Python + + with openmc.capi.run_in_memory(): + openmc.capi.simulation_init() + for _ in openmc.capi.iter_batches(): + # Look at convergence of tallies, for example + ... + openmc.capi.simulation_finalize() + + See Also + -------- + openmc.capi.next_batch + + """ while True: # Run next batch retval = _dll.openmc_next_batch() diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index 14f9d1c1c..6725ee5c9 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -4,6 +4,7 @@ from weakref import WeakValueDictionary import numpy as np from numpy.ctypeslib import as_array +import scipy.stats from openmc.data.reaction import REACTION_NAME from . import _dll, Nuclide @@ -31,6 +32,9 @@ _dll.openmc_tally_get_filters.argtypes = [ c_int32, POINTER(POINTER(c_int32)), POINTER(c_int)] _dll.openmc_tally_get_filters.restype = c_int _dll.openmc_tally_get_filters.errcheck = _error_handler +_dll.openmc_tally_get_n_realizations.argtypes = [c_int32, POINTER(c_int32)] +_dll.openmc_tally_get_n_realizations.restype = c_int +_dll.openmc_tally_get_n_realizations.errcheck = _error_handler _dll.openmc_tally_get_nuclides.argtypes = [ c_int32, POINTER(POINTER(c_int)), POINTER(c_int)] _dll.openmc_tally_get_nuclides.restype = c_int @@ -70,6 +74,14 @@ _SCORES = { def global_tallies(): + """Mean and standard deviation of the mean for each global tally. + + Returns + ------- + list of tuple + For each global tally, a tuple of (mean, standard deviation) + + """ ptr = POINTER(c_double)() _dll.openmc_global_tallies(ptr) array = as_array(ptr, (4, 3)) @@ -113,10 +125,16 @@ class Tally(_FortranObjectWithID): ID of the tally filters : list List of tally filters + mean : numpy.ndarray + An array containing the sample mean for each bin nuclides : list of str List of nuclides to score results for + num_realizations : int + Number of realizations results : numpy.ndarray Array of tally results + std_dev : numpy.ndarray + An array containing the sample standard deviation for each bin """ __instances = WeakValueDictionary() @@ -169,21 +187,6 @@ class Tally(_FortranObjectWithID): _dll.openmc_tally_get_filters(self._index, filt_idx, n) return [_get_filter(filt_idx[i]) for i in range(n.value)] - @property - def nuclides(self): - nucs = POINTER(c_int)() - n = c_int() - _dll.openmc_tally_get_nuclides(self._index, nucs, n) - return [Nuclide(nucs[i]).name if nucs[i] > 0 else 'total' - for i in range(n.value)] - - @property - def results(self): - data = POINTER(c_double)() - shape = (c_int*3)() - _dll.openmc_tally_results(self._index, data, shape) - return as_array(data, tuple(shape[::-1])) - @filters.setter def filters(self, filters): # Get filter indices as int32_t[] @@ -192,12 +195,42 @@ class Tally(_FortranObjectWithID): _dll.openmc_tally_set_filters(self._index, n, indices) + @property + def mean(self): + n = self.num_realizations + sum_ = self.results[:, :, 1] + if n > 0: + return sum_ / n + else: + return sum_.copy() + + @property + def nuclides(self): + nucs = POINTER(c_int)() + n = c_int() + _dll.openmc_tally_get_nuclides(self._index, nucs, n) + return [Nuclide(nucs[i]).name if nucs[i] > 0 else 'total' + for i in range(n.value)] + @nuclides.setter def nuclides(self, nuclides): nucs = (c_char_p * len(nuclides))() nucs[:] = [x.encode() for x in nuclides] _dll.openmc_tally_set_nuclides(self._index, len(nuclides), nucs) + @property + def num_realizations(self): + n = c_int32() + _dll.openmc_tally_get_n_realizations(self._index, n) + return n.value + + @property + def results(self): + data = POINTER(c_double)() + shape = (c_int*3)() + _dll.openmc_tally_results(self._index, data, shape) + return as_array(data, tuple(shape[::-1])) + @property def scores(self): scores_as_int = POINTER(c_int)() @@ -223,21 +256,47 @@ class Tally(_FortranObjectWithID): scores_[:] = [x.encode() for x in scores] _dll.openmc_tally_set_scores(self._index, len(scores), scores_) - @classmethod - def new(cls, tally_id=None): - # Determine ID to assign - if tally_id is None: - try: - tally_id = max(tallies) + 1 - except ValueError: - tally_id = 1 + @property + def std_dev(self): + results = self.results + std_dev = np.empty(results.shape[:2]) + std_dev.fill(np.nan) - index = c_int32() - _dll.openmc_extend_tallies(1, index, None) - _dll.openmc_tally_set_type(index, b'generic') - tally = cls(index.value) - tally.id = tally_id - return tally + n = self.num_realizations + if n > 1: + # Get sum and sum-of-squares from results + sum_ = results[:, :, 1] + sum_sq = results[:, :, 2] + + # Determine non-zero entries + mean = sum_ / n + nonzero = np.abs(mean) > 0 + + # Calculate sample standard deviation of the mean + std_dev[nonzero] = np.sqrt( + (sum_sq[nonzero]/n - mean[nonzero]**2)/(n - 1)) + + return std_dev + + def ci_width(self, alpha=0.05): + """Confidence interval half-width based on a Student t distribution + + Parameters + ---------- + alpha : float + Significance level (one minus the confidence level!) + + Returns + ------- + float + Half-width of a two-sided (1 - :math:`alpha`) confidence interval + + """ + half_width = self.std_dev.copy() + n = self.num_realizations + if n > 1: + half_width *= scipy.stats.t.ppf(1 - alpha/2, n - 1) + return half_width class _TallyMapping(Mapping): diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 69b519ea5..766dbfa9f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -863,7 +863,7 @@ contains ! Preallocate space for keff and entropy by generation call k_generation % reserve(n_max_batches*gen_per_batch) - call entropy % initialize(n_max_batches*gen_per_batch) + call entropy % reserve(n_max_batches*gen_per_batch) ! Get the trigger information for keff if (check_for_node(node_base, "keff_trigger")) then diff --git a/src/state_point.F90 b/src/state_point.F90 index e93d6e472..0383150ac 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -632,6 +632,7 @@ contains subroutine load_state_point() integer :: i + integer :: n integer :: int_array(3) integer, allocatable :: array(:) integer(HID_T) :: file_id @@ -710,11 +711,14 @@ contains if (run_mode == MODE_EIGENVALUE) then call read_dataset(int_array(1), file_id, "n_inactive") call read_dataset(gen_per_batch, file_id, "generations_per_batch") - call read_dataset(k_generation % data(1:restart_batch*gen_per_batch), & - file_id, "k_generation") + + n = restart_batch*gen_per_batch + call k_generation % resize(n) + call read_dataset(k_generation % data(1:n), file_id, "k_generation") + if (entropy_on) then - call read_dataset(entropy % data(1:restart_batch*gen_per_batch), & - file_id, "entropy") + call entropy % resize(n) + call read_dataset(entropy % data(1:n), file_id, "entropy") end if call read_dataset(k_col_abs, file_id, "k_col_abs") call read_dataset(k_col_tra, file_id, "k_col_tra") diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 2812eef6e..1778474e0 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -526,6 +526,22 @@ contains end function openmc_tally_get_filters + function openmc_tally_get_n_realizations(index, n) result(err) bind(C) + ! Return the number of realizations for a tally + integer(C_INT32_T), value :: index + integer(C_INT32_T), intent(out) :: n + integer(C_INT) :: err + + if (index >= 1 .and. index <= size(tallies)) then + n = tallies(index) % obj % n_realizations + err = 0 + else + err = E_OUT_OF_BOUNDS + call set_errmsg('Index in tallies array is out of bounds.') + end if + end function openmc_tally_get_n_realizations + + function openmc_tally_get_nuclides(index, nuclides, n) result(err) bind(C) ! Return the list of nuclides assigned to a tally integer(C_INT32_T), value :: index