From 756a43f1c563276af124b520a974389f8b286114 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 28 Jan 2017 11:34:29 -0500 Subject: [PATCH 01/29] Added ability to convert the isotropic or angle representation from one to the other for XSdata objects and MGXSLibrary objects --- openmc/mgxs/mgxs.py | 2 +- openmc/mgxs_library.py | 196 ++++++++++++++++++++++++++++++++++++++++- 2 files changed, 194 insertions(+), 4 deletions(-) diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index de90d4397..9f1ac59d1 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -205,7 +205,7 @@ class MGXS(object): clone._name = self.name clone._rxn_type = self.rxn_type clone._by_nuclide = self.by_nuclide - clone._nuclides = copy.deepcopy(self._nuclides) + clone._nuclides = copy.deepcopy(self._nuclides, memo) clone._domain = self.domain clone._domain_type = self.domain_type clone._energy_groups = copy.deepcopy(self.energy_groups, memo) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 4fb347b73..09f7e421b 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1,4 +1,4 @@ -from collections import Iterable +import copy from numbers import Real, Integral import os @@ -185,6 +185,52 @@ class XSdata(object): self._inverse_velocity = len(temperatures) * [None] self._xs_shapes = None + def __deepcopy__(self, memo): + existing = memo.get(id(self)) + + # If this is the first time we have tried to copy this object, copy it + if existing is None: + clone = type(self).__new__(type(self)) + clone._name = self.name + clone._energy_groups = copy.deepcopy(self.energy_groups, memo) + clone._num_delayed_groups = self.num_delayed_groups + clone._temperatures = copy.deepcopy(self.temperatures, memo) + clone._drepresentation = self.representation + clone._atomic_weight_ratio = self._atomic_weight_ratio + clone._fissionable = self._fissionable + clone._scatter_format = self._scatter_format + clone._order = self._order + clone._num_polar = self._num_polar + clone._num_azimuthal = self._num_azimuthal + clone._total = copy.deepcopy(self._total, memo) + clone._absorption = copy.deepcopy(self._absorption, memo) + clone._scatter_matrix = copy.deepcopy(self._scatter_matrix, memo) + clone._multiplicity_matrix = \ + copy.deepcopy(self._multiplicity_matrix, memo) + clone._fission = copy.deepcopy(self._fission, memo) + clone._nu_fission = copy.deepcopy(self._nu_fission, memo) + clone._prompt_nu_fission = \ + copy.deepcopy(self._prompt_nu_fission, memo) + clone._delayed_nu_fission = \ + copy.deepcopy(self._delayed_nu_fission, memo) + clone._kappa_fission = copy.deepcopy(self._kappa_fission, memo) + clone._chi = copy.deepcopy(self._chi, memo) + clone._chi_prompt = copy.deepcopy(self._chi_prompt, memo) + clone._chi_delayed = copy.deepcopy(self._chi_delayed, memo) + clone._beta = copy.deepcopy(self._beta, memo) + clone._decay_rate = copy.deepcopy(self._decay_rate, memo) + clone._inverse_velocity = \ + copy.deepcopy(self._inverse_velocity, memo) + clone._xs_shapes = copy.deepcopy(self._xs_shapes, memo) + + memo[id(self)] = clone + + return clone + + # If this object has been copied before, return the first copy made + else: + return existing + @property def name(self): return self._name @@ -348,7 +394,7 @@ class XSdata(object): @representation.setter def representation(self, representation): - # Check it is of valid type. + # Check it is of valid value. check_value('representation', representation, _REPRESENTATIONS) self._representation = representation @@ -1622,7 +1668,8 @@ class XSdata(object): """ - check_type('inverse_velocity', inverse_velocity, openmc.mgxs.InverseVelocity) + check_type('inverse_velocity', inverse_velocity, + openmc.mgxs.InverseVelocity) check_value('energy_groups', inverse_velocity.energy_groups, [self.energy_groups]) check_value('domain_type', inverse_velocity.domain_type, @@ -1634,6 +1681,89 @@ class XSdata(object): self._inverse_velocity[i] = inverse_velocity.get_xs( nuclides=nuclide, xs_type=xs_type, subdomains=subdomain) + def convert_representation(self, target_representation, num_polar=None, + num_azimuthal=None): + """Produce a new XSdata object with the same data, but converted to the + new representation + + Parameters + ---------- + target_representation : {'isotropic', 'angle'} + Representation of the MGXS (isotropic or angle-dependent flux + weighting). + num_polar : int, optional + Number of equal width angular bins that the polar angular + domain is subdivided into. This is required when + :param:`target_representation` is "angle". + num_azimuthal : int, optional + Number of equal width angular bins that the azimuthal angular + domain is subdivided into. This is required when + :param:`target_representation` is "angle". + + Returns + ------- + openmc.XSdata + Multi-group cross section data with the same data as self, but + represented as specified in :param:`target_representation`. + + """ + + check_value('target_representation', target_representation, + _REPRESENTATIONS) + if target_representation == 'angle': + check_type('num_polar', num_polar, Integral) + check_type('num_azimuthal', num_azimuthal, Integral) + check_greater_than('num_polar', num_polar, 0) + check_greater_than('num_azimuthal', num_azimuthal, 0) + + xsdata = copy.deepcopy(self) + if target_representation == self.representation: + # Check to make sure the num_polar and num_azimuthal values match + if target_representation == 'angle': + if num_polar != self.num_polar or num_azimuthal != self.num_azimuthal: + raise NotImplementedError("XSdata.convert_representation " + "cannot translate between " + "`angle` representations with " + "different angle bin structures") + # Nothing to do + return xsdata + + types = ['total', 'absorption', 'fission', 'nu_fission', + 'scatter_matrix', 'multiplicity_matrix', 'prompt_nu_fission', + 'delayed_nu_fission', 'kappa_fission', 'chi', 'chi_prompt', + 'chi_delayed', 'beta', 'decay_rate', 'inverse_velocity'] + + xsdata.representation = target_representation + # We have different actions depending on the representation conversion + if target_representation == 'isotropic': + xsdata.num_polar = None + xsdata.num_azimuthal = None + for i, temp in enumerate(xsdata.temperatures): + for xs in types: + # Get the original data + orig_data = getattr(self, '_' + xs)[i] + # Since we are going from angle to isotropic, the current + # data should be just averaged over the angle bins + new_data = orig_data.mean(axis=(0, 1)) + setter = getattr(xsdata, 'set_' + xs) + setter(new_data, temp) + + elif target_representation == 'angle': + xsdata.num_polar = num_polar + xsdata.num_azimuthal = num_azimuthal + for i, temp in enumerate(xsdata.temperatures): + for xs in types: + # Get the original data + orig_data = getattr(self, '_' + xs)[i] + # Since we are going from isotropic to angle, the current + # data should be just copied for every polar/azimuthal bin + new_shape = (num_polar, num_azimuthal) + orig_data.shape + new_data = np.resize(orig_data, new_shape) + setter = getattr(xsdata, 'set_' + xs) + setter(new_data, temp) + + return xsdata + def to_hdf5(self, file): """Write XSdata to an HDF5 file @@ -1995,6 +2125,24 @@ class MGXSLibrary(object): self.num_delayed_groups = num_delayed_groups self._xsdatas = [] + def __deepcopy__(self, memo): + existing = memo.get(id(self)) + + # If this is the first time we have tried to copy this object, copy it + if existing is None: + clone = type(self).__new__(type(self)) + clone._energy_groups = copy.deepcopy(self.energy_groups, memo) + clone._num_delayed_groups = self.num_delayed_groups + clone._xsdatas = copy.deepcopy(self.xsdatas, memo) + + memo[id(self)] = clone + + return clone + + # If this object has been copied before, return the first copy made + else: + return existing + @property def energy_groups(self): return self._energy_groups @@ -2099,6 +2247,48 @@ class MGXSLibrary(object): result = xsdata return result + def convert_representation(self, target_representation, num_polar=None, + num_azimuthal=None): + """Produce a new MGXSLibrary object with the same data, but converted + to the new representation + + Parameters + ---------- + target_representation : {'isotropic', 'angle'} + Representation of the MGXS (isotropic or angle-dependent flux + weighting). + num_polar : int, optional + Number of equal width angular bins that the polar angular + domain is subdivided into. This is required when + :param:`target_representation` is "angle". + num_azimuthal : int, optional + Number of equal width angular bins that the azimuthal angular + domain is subdivided into. This is required when + :param:`target_representation` is "angle". + + Returns + ------- + openmc.MGXSLibrary + Multi-group Library with the same data as self, but represented as + specified in :param:`target_representation`. + + """ + + check_value('target_representation', target_representation, + _REPRESENTATIONS) + if target_representation == 'angle': + check_type('num_polar', num_polar, Integral) + check_type('num_azimuthal', num_azimuthal, Integral) + check_greater_than('num_polar', num_polar, 0) + check_greater_than('num_azimuthal', num_azimuthal, 0) + + library = copy.deepcopy(self) + for i, xsdata in enumerate(self.xsdatas): + library.xsdatas[i] = \ + xsdata.convert_representation(target_representation, + num_polar, num_azimuthal) + return library + def export_to_hdf5(self, filename='mgxs.h5'): """Create an hdf5 file that can be used for a simulation. From 1cd15696c17eb1c4bb63926aad5b0a4efc2a7ed3 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 28 Jan 2017 13:44:38 -0500 Subject: [PATCH 02/29] fixes to convert_representation and MGXSLibrary implementation of convert_scatter_format --- openmc/mgxs_library.py | 83 ++++++++++++++++++++++++++++++++---------- 1 file changed, 63 insertions(+), 20 deletions(-) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 09f7e421b..5a23e6dbe 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1736,29 +1736,33 @@ class XSdata(object): xsdata.representation = target_representation # We have different actions depending on the representation conversion if target_representation == 'isotropic': - xsdata.num_polar = None - xsdata.num_azimuthal = None - for i, temp in enumerate(xsdata.temperatures): - for xs in types: - # Get the original data - orig_data = getattr(self, '_' + xs)[i] - # Since we are going from angle to isotropic, the current - # data should be just averaged over the angle bins - new_data = orig_data.mean(axis=(0, 1)) - setter = getattr(xsdata, 'set_' + xs) - setter(new_data, temp) - + # This is not needed for the correct functionality, but these + # values are changed back to None for clarity + xsdata._num_polar = None + xsdata._num_azimuthal = None elif target_representation == 'angle': xsdata.num_polar = num_polar xsdata.num_azimuthal = num_azimuthal - for i, temp in enumerate(xsdata.temperatures): - for xs in types: - # Get the original data - orig_data = getattr(self, '_' + xs)[i] - # Since we are going from isotropic to angle, the current - # data should be just copied for every polar/azimuthal bin - new_shape = (num_polar, num_azimuthal) + orig_data.shape - new_data = np.resize(orig_data, new_shape) + + # Reset XSdata.xs_shapes to accomodate the new shape + # import pdb; pdb.set_trace() + xsdata._xs_shapes = None + xsdata.xs_shapes + + for i, temp in enumerate(xsdata.temperatures): + for xs in types: + # Get the original data + orig_data = getattr(self, '_' + xs)[i] + if orig_data is not None: + if target_representation == 'isotropic': + # Since we are going from angle to isotropic, the + # current data is just averaged over the angle bins + new_data = orig_data.mean(axis=(0, 1)) + elif target_representation == 'angle': + # Since we are going from isotropic to angle, the + # current data is just copied for every angle bin + new_shape = (num_polar, num_azimuthal) + orig_data.shape + new_data = np.resize(orig_data, new_shape) setter = getattr(xsdata, 'set_' + xs) setter(new_data, temp) @@ -2289,6 +2293,45 @@ class MGXSLibrary(object): num_polar, num_azimuthal) return library + def convert_scatter_format(self, target_format, target_order): + """Produce a new MGXSLibrary object with the same data, but converted + to the new scatter format and order + + Parameters + ---------- + target_format : {'isotropic', 'angle'} + Representation of the MGXS (isotropic or angle-dependent flux + weighting). + target_order : int + Either the Legendre target_order, number of bins, or number of + points used to describe the angular distribution associated with + each group-to-group transfer probability. + + Returns + ------- + openmc.MGXSLibrary + Multi-group Library with the same data as self, but with the + scatter format represented as specified in :param:`target_format` + and :param:`target_order`. + + """ + + check_value('target_format', target_format, _SCATTER_TYPES) + check_type('target_order', target_order, Integral) + if target_format == 'legendre': + check_greater_than('target_order', target_order, Integral, 0, + equality=True) + else: + check_greater_than('target_order', target_order, Integral, 0, + equality=True) + + library = copy.deepcopy(self) + for i, xsdata in enumerate(self.xsdatas): + library.xsdatas[i] = \ + xsdata.convert_scatter_format(target_format, target_order) + + return library + def export_to_hdf5(self, filename='mgxs.h5'): """Create an hdf5 file that can be used for a simulation. From 1c367b049c77cf945565187fb384a5a1b49be0dc Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 29 Jan 2017 15:25:00 -0500 Subject: [PATCH 03/29] Finished convert_scatter_format; tested, seems to be great! --- openmc/mgxs_library.py | 206 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 193 insertions(+), 13 deletions(-) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 5a23e6dbe..b740b4dae 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -14,9 +14,13 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \ # Supported incoming particle MGXS angular treatment representations _REPRESENTATIONS = ['isotropic', 'angle'] +# Supported scattering angular distribution representations _SCATTER_TYPES = ['tabular', 'legendre', 'histogram'] +# List of MGXS dimension types _XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[DG][G]", "[DG][G']", "[DG][G][G']"] +# Number of mu points for conversion between scattering formats +_NMU = 257 class XSdata(object): @@ -195,7 +199,7 @@ class XSdata(object): clone._energy_groups = copy.deepcopy(self.energy_groups, memo) clone._num_delayed_groups = self.num_delayed_groups clone._temperatures = copy.deepcopy(self.temperatures, memo) - clone._drepresentation = self.representation + clone._representation = self.representation clone._atomic_weight_ratio = self._atomic_weight_ratio clone._fissionable = self._fissionable clone._scatter_format = self._scatter_format @@ -1745,9 +1749,7 @@ class XSdata(object): xsdata.num_azimuthal = num_azimuthal # Reset XSdata.xs_shapes to accomodate the new shape - # import pdb; pdb.set_trace() xsdata._xs_shapes = None - xsdata.xs_shapes for i, temp in enumerate(xsdata.temperatures): for xs in types: @@ -1761,13 +1763,194 @@ class XSdata(object): elif target_representation == 'angle': # Since we are going from isotropic to angle, the # current data is just copied for every angle bin - new_shape = (num_polar, num_azimuthal) + orig_data.shape + new_shape = (num_polar, num_azimuthal) + \ + orig_data.shape new_data = np.resize(orig_data, new_shape) setter = getattr(xsdata, 'set_' + xs) setter(new_data, temp) return xsdata + def convert_scatter_format(self, target_format, target_order=None): + """Produce a new MGXSLibrary object with the same data, but converted + to the new scatter format and order + + Parameters + ---------- + target_format : {'tabular', 'legendre', 'histogram'} + Representation of the scattering angle distribution + target_order : int + Either the Legendre target_order, number of bins, or number of + points used to describe the angular distribution associated with + each group-to-group transfer probability + + Returns + ------- + openmc.XSdata + Multi-group cross section data with the same data as self, but + represented as specified in :param:`target_format`. + + """ + + from scipy.interpolate import interp1d + from scipy.integrate import quad, simps + from scipy.special import eval_legendre + + check_value('target_format', target_format, _SCATTER_TYPES) + check_type('target_order', target_order, Integral) + if target_format == 'legendre': + check_greater_than('target_order', target_order, 0, equality=True) + else: + check_greater_than('target_order', target_order, 0) + + xsdata = copy.deepcopy(self) + xsdata.scatter_format = target_format + xsdata.order = target_order + # Reset XSdata.xs_shapes to accomodate the new shape + xsdata._xs_shapes = None + xsdata.xs_shapes + + # We have to accomodate the following possibilities: + # histogram -> tabular w/ same or diff order + # histogram -> legendre + # histogram -> histogram w/ same or diff order + + for i, temp in enumerate(xsdata.temperatures): + orig_data = self._scatter_matrix[i] + new_shape = orig_data.shape[:-1] + (xsdata.num_orders,) + new_data = np.zeros(new_shape) + if self.scatter_format == 'legendre': + if target_format == 'legendre': + # Then we are changing orders and only need to change + # dimensionality of the mu data and pad/truncate as needed + order = min(xsdata.num_orders, self.num_orders) + new_data[..., :order] = orig_data[..., :order] + elif target_format == 'tabular': + mu = np.linspace(-1, 1, xsdata.num_orders) + # Create the Legendre series and either do a point-wise + # evaluation (tabular) or bin-wise integral (histogram) + new_data[..., :] = \ + np.sum((l + 0.5) * eval_legendre(l, mu) * orig_data + for l in range(self.num_orders)) + elif target_format == 'histogram': + mu = np.linspace(-1, 1, xsdata.num_orders + 1) + # This code will be written to utilize the vectorized + # integration capabilities instead of having an isotropic + # and angle representation path. + for h_bin in range(xsdata.num_orders): + mu_fine = np.linspace(mu[i], mu[i + 1], _NMU) + new_data[..., h_bin] = \ + simps(np.sum((l + 0.5) * + eval_legendre(l, mu_fine) * + orig_data + for l in range(self.num_orders)), + mu_fine) + + # Remove the very small results from numerical precision + # issues (allowing conversions to be reproduced exactly) + new_data[..., np.abs(new_data) < 1.E-10] = 0. + + elif self.scatter_format == 'tabular': + mu_self = np.linspace(-1, 1, self.num_orders) + if target_format == 'legendre': + # find the legendre coefficients via integration. To best + # use the vectorized integration capabilities of scipy, + # this will be done with fixed sample integration routines. + mu_fine = np.linspace(-1, 1, _NMU) + y = [interp1d(mu_self, orig_data)(mu_fine) * + eval_legendre(l, mu_fine) + for l in range(xsdata.num_orders)] + for l in range(xsdata.num_orders): + new_data[..., l] = (l + 0.5) * simps(y[l], mu_fine) + + # Remove the very small results from numerical precision + # issues (allowing conversions to be reproduced exactly) + new_data[..., np.abs(new_data) < 1.E-10] = 0. + elif target_format == 'tabular': + # Simply use an interpolating function to get the new data + mu = np.linspace(-1, 1, xsdata.num_orders) + new_data[..., :] = interp1d(mu_self, orig_data)(mu) + elif target_format == 'histogram': + # Use an interpolating function to do the bin-wise + # integrals + mu = np.linspace(-1, 1, xsdata.num_orders + 1) + + # Like the tabular -> legendre path above, this code will + # be written to utilize the vectorized integration + # capabilities instead of having an isotropic and + # angle representation path. + interp = interp1d(mu_self, orig_data) + for h_bin in range(xsdata.num_orders): + mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU) + new_data[..., h_bin] = simps(interp(mu_fine), mu_fine) + + # Remove the very small results from numerical precision + # issues (allowing conversions to be reproduced exactly) + new_data[..., np.abs(new_data) < 1.E-10] = 0. + elif self.scatter_format == 'histogram': + # The histogram format does not have enough information to + # convert to the other forms without inducing an error. + # We will make the assumption that the left-edge of the bin + # has the value of the bin. The mu=1 value will be extrapolated + # from the previous two points. + mu_self = np.linspace(-1, 1, self.num_orders + 1) + tab_data = np.zeros(orig_data.shape[:-1] + + (orig_data.shape[-1] + 1,)) + tab_data[..., :-1] = orig_data + tab_data[..., -1] = (orig_data[..., -2] - + orig_data[..., -3]) + orig_data[..., -2] + # Now get the distribution normalization factor.such + # This factor will be such that its average value is the + # group -> group cross section (which is the sum of the + # original data) + norm = np.sum(orig_data, axis=-1) / np.mean(tab_data, axis=-1) + norm = np.nan_to_num(norm) + for i in range(tab_data.shape[-1]): + tab_data[..., i] *= norm[...] + + # We now have a tabular distribution in tab_data on mu_self + # with a normalization in norm. We now proceed just like the + # tabular branch above. + if target_format == 'legendre': + # find the legendre coefficients via integration. To best + # use the vectorized integration capabilities of scipy, + # this will be done with fixed sample integration routines. + mu_fine = np.linspace(-1, 1, _NMU) + y = [interp1d(mu_self, tab_data)(mu_fine) * + eval_legendre(l, mu_fine) + for l in range(xsdata.num_orders)] + for l in range(xsdata.num_orders): + new_data[..., l] = (l + 0.5) * simps(y[l], mu_fine) + + # Remove the very small results from numerical precision + # issues (allowing conversions to be reproduced exactly) + new_data[..., np.abs(new_data) < 1.E-10] = 0. + elif target_format == 'tabular': + # Simply use an interpolating function to get the new data + mu = np.linspace(-1, 1, xsdata.num_orders) + new_data[..., :] = interp1d(mu_self, tab_data)(mu) + elif target_format == 'histogram': + # Use an interpolating function to do the bin-wise + # integrals + mu = np.linspace(-1, 1, xsdata.num_orders + 1) + + # Like the tabular -> legendre path above, this code will + # be written to utilize the vectorized integration + # capabilities instead of having an isotropic and + # angle representation path. + interp = interp1d(mu_self, tab_data) + for h_bin in range(xsdata.num_orders): + mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU) + new_data[..., h_bin] = simps(interp(mu_fine), mu_fine) + + # Remove the very small results from numerical precision + # issues (allowing conversions to be reproduced exactly) + new_data[..., np.abs(new_data) < 1.E-10] = 0. + + xsdata.set_scatter_matrix(new_data, temp) + + return xsdata + def to_hdf5(self, file): """Write XSdata to an HDF5 file @@ -1891,7 +2074,7 @@ class XSdata(object): elif self.representation == 'angle': matrix = \ self._scatter_matrix[i][p, a, g_in, :, 0] - elif self.scatter_format == 'histogram': + else: if self.representation == 'isotropic': matrix = \ np.sum(self._scatter_matrix[i][g_in, :, :], @@ -2299,13 +2482,12 @@ class MGXSLibrary(object): Parameters ---------- - target_format : {'isotropic', 'angle'} - Representation of the MGXS (isotropic or angle-dependent flux - weighting). + target_format : {'tabular', 'legendre', 'histogram'} + Representation of the scattering angle distribution target_order : int Either the Legendre target_order, number of bins, or number of points used to describe the angular distribution associated with - each group-to-group transfer probability. + each group-to-group transfer probability Returns ------- @@ -2319,11 +2501,9 @@ class MGXSLibrary(object): check_value('target_format', target_format, _SCATTER_TYPES) check_type('target_order', target_order, Integral) if target_format == 'legendre': - check_greater_than('target_order', target_order, Integral, 0, - equality=True) + check_greater_than('target_order', target_order, 0, equality=True) else: - check_greater_than('target_order', target_order, Integral, 0, - equality=True) + check_greater_than('target_order', target_order, 0) library = copy.deepcopy(self) for i, xsdata in enumerate(self.xsdatas): From 7e1700227c6bf7dc43e122813f2f12423e71f234 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 5 Feb 2017 15:56:23 -0500 Subject: [PATCH 04/29] Fixes to scatter format conversion as well as initial implementation of conversion of MG library to CE --- openmc/mgxs_library.py | 253 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 239 insertions(+), 14 deletions(-) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index b740b4dae..6f46f4887 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1793,7 +1793,7 @@ class XSdata(object): """ from scipy.interpolate import interp1d - from scipy.integrate import quad, simps + from scipy.integrate import simps from scipy.special import eval_legendre check_value('target_format', target_format, _SCATTER_TYPES) @@ -1827,24 +1827,34 @@ class XSdata(object): new_data[..., :order] = orig_data[..., :order] elif target_format == 'tabular': mu = np.linspace(-1, 1, xsdata.num_orders) - # Create the Legendre series and either do a point-wise - # evaluation (tabular) or bin-wise integral (histogram) - new_data[..., :] = \ - np.sum((l + 0.5) * eval_legendre(l, mu) * orig_data - for l in range(self.num_orders)) + # Evaluate the legendre on the tabular grid within mu + for imu in range(len(mu)): + new_data[..., imu] = \ + np.sum((l + 0.5) * eval_legendre(l, mu[imu]) * + orig_data[..., l] + for l in range(self.num_orders)) elif target_format == 'histogram': mu = np.linspace(-1, 1, xsdata.num_orders + 1) # This code will be written to utilize the vectorized # integration capabilities instead of having an isotropic # and angle representation path. for h_bin in range(xsdata.num_orders): - mu_fine = np.linspace(mu[i], mu[i + 1], _NMU) - new_data[..., h_bin] = \ - simps(np.sum((l + 0.5) * - eval_legendre(l, mu_fine) * - orig_data - for l in range(self.num_orders)), - mu_fine) + mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU) + table_shape = new_data.shape[:-1] + (_NMU,) + table_fine = np.zeros(table_shape) + for imu in range(len(mu_fine)): + table_fine[..., imu] = \ + np.sum((l + 0.5) * + eval_legendre(l, mu_fine[imu]) * + orig_data[..., l] + for l in range(self.num_orders)) + new_data[..., h_bin] = simps(table_fine, mu_fine) + # new_data[..., h_bin] = \ + # simps(np.sum((l + 0.5) * + # eval_legendre(l, mu_fine) * + # orig_data[..., l] + # for l in range(self.num_orders)), + # mu_fine) # Remove the very small results from numerical precision # issues (allowing conversions to be reproduced exactly) @@ -1899,7 +1909,7 @@ class XSdata(object): tab_data[..., :-1] = orig_data tab_data[..., -1] = (orig_data[..., -2] - orig_data[..., -3]) + orig_data[..., -2] - # Now get the distribution normalization factor.such + # Now get the distribution normalization factor. # This factor will be such that its average value is the # group -> group cross section (which is the sum of the # original data) @@ -1951,6 +1961,198 @@ class XSdata(object): return xsdata + def convert_to_continuous_energy(self, temperature=294.): + """Converts the XSdata object to an equivalent + openmc.data.IncidentNeutron object. + + Parameters + ---------- + temperature : float + Temperature of dataset to print; defaults to 294K + + Returns + ------- + openmc.data.IncidentNeutron + The continuous-energy IncidentNeutron data library equivalent + to the MGXS data in self + """ + + # Check if this can be performed successfully + if self.representation == 'angle': + raise ValueError("Cannot convert angle-dependent MGXS; convert to " + "an isotropic representation first") + required_types = ['absorption', 'scatter_matrix'] + for type_check in required_types: + if getattr(self, '_' + type_check)[0] is None: + raise ValueError(type_check + ' data is required') + if temperature not in self.temperatures: + raise ValueError("Invalid temperature") + + def convert_xs(group_edges, values): + cexs = openmc.data.Tabulated1D(group_edges, + np.append(values[::-1], + [values[0]]), + breakpoints=[len(group_edges)], + interpolation=[1]) + return cexs + + # Build required metadata + kTs = self.temperatures[:] * openmc.data.K_BOLTZMANN + if self.atomic_weight_ratio: + awr = self.atomic_weight_ratio + else: + awr = 1. + + # Get the temperature index + iT = np.where(self.temperatures == temperature)[0][0] + data = openmc.data.IncidentNeutron(self.name, 1, 1, 0, awr, kTs) + strT = "{}K".format(int(round(self.temperatures[iT]))) + data.energy = {strT: self.energy_groups.group_edges} + energy_midpoints = (self.energy_groups.group_edges[1:] + + self.energy_groups.group_edges[:-1]) / 2. + + # Elastic MGXS: must explicitly be 0 barns to avoid incorrect + # elastic_scatter calculation with data available to us in MG library. + el_rxn = openmc.data.Reaction(2) + el_rxn.xs[strT] = \ + convert_xs(self.energy_groups.group_edges, + np.zeros(self.energy_groups.num_groups)) + data.reactions[2] = el_rxn + + # Absorption MGXS + abs_rxn = openmc.data.Reaction(102) + abs_rxn.xs[strT] = convert_xs(self.energy_groups.group_edges, + np.subtract(self._absorption[iT], + self._fission[iT])) + data.reactions[102] = abs_rxn + + if self._nu_fission[iT] is not None and self._fission[iT] is not None: + fiss_rxn = openmc.data.Reaction(18) + fiss_rxn.xs[strT] = convert_xs(self.energy_groups.group_edges, + self._fission[iT]) + + # Get nu_fission and chi from the presence of both, or the + # nu_fission matrix + if self._nu_fission[iT].shape == self.xs_shapes["[G][G']"]: + nu_fiss = np.sum(self._nu_fission[iT], axis=1)[::-1] + chi = self._nu_fission[iT][::-1] / nu_fiss + else: + nu_fiss = self._nu_fission[iT][::-1] + chi = np.reshape(np.tile(self._chi[iT][::-1], + self.energy_groups.num_groups), + (self.energy_groups.num_groups, + self.energy_groups.num_groups)) + nu = np.divide(nu_fiss, self._fission[iT][::-1]) + + # Build a histogram distribution to represent the yield for the + # incoming groups + prod = openmc.data.Product() + prod.yield_ = \ + openmc.data.Tabulated1D(self.energy_groups.group_edges[:-1], + nu, breakpoints=[len(nu)], + interpolation=[1]) + + # Now build the outgoing energy distribution using discrete energy + # lines within each of the outgoing groups + chi_eouts = [] + for g in range(self.energy_groups.num_groups): + chi_eouts.append(openmc.stats.Discrete(energy_midpoints, + chi[g, :])) + # Ensure the distribution CDF starts with 0 + chi_eouts[-1].c = np.cumsum(chi[g, :]) - chi[g, 0] + + # Create the continuous distribution for the fission + # The angular distribution will remain None (thus isotropic) + chi_distrib = openmc.data.ContinuousTabular( + [len(chi)], [1], self.energy_groups.group_edges[:-1], + chi_eouts) + prod.distribution = \ + [openmc.data.UncorrelatedAngleEnergy(energy=chi_distrib)] + fiss_rxn.products = [prod] + fiss_rxn.center_of_mass = False + data.reactions[18] = fiss_rxn + + # Scattering Data + # First convert the data to a tabular representation + if self.scatter_format == 'tabular': + tabular = self + else: + tabular = self.convert_scatter_format('tabular', 33) + + # Calculate the isotropic scattering matrix and use that to find the + # total scattering x/s and outgoing energy distributions + isotropic_matrix = np.mean(tabular._scatter_matrix[iT], axis=-1)[::-1, + ::-1] + scatt_xs = np.sum(isotropic_matrix, axis=1) + energy = np.zeros((self.energy_groups.num_groups, + self.energy_groups.num_groups)) + for gin in range(self.energy_groups.num_groups): + energy[gin, :] = isotropic_matrix[gin, :] / scatt_xs[gin] + + # Get the anisotropic but normalized angular distribution + distrib = np.zeros((self.energy_groups.num_groups, + self.energy_groups.num_groups, + tabular.num_orders)) + for gin in range(self.energy_groups.num_groups): + for gout in range(self.energy_groups.num_groups): + distrib[gin, gout, :] = \ + np.divide(tabular._scatter_matrix[iT][gin, gout, :], + isotropic_matrix[gin, gout]) + distrib = np.nan_to_num(distrib) + + # Incorporate the scattering multiplication, if required + scatt_prod = openmc.data.Product() + if self._multiplicity_matrix[iT]: + yield_ = self._multiplicity_matrix[iT][::-1, ::-1] + scatt_prod.yield_ = \ + openmc.data.Tabulated1D(self.energy_groups.group_edges[:-1], + yield_, breakpoints=[len(yield_)], + interpolation=[1]) + + # Now build the outgoing energy distribution using discrete energy + # lines within each of the outgoing groups + scatt_eouts = [] + for g in range(self.energy_groups.num_groups): + scatt_eouts.append( + openmc.stats.Discrete(energy_midpoints, energy[g, :])) + # Ensure the distribution CDF starts with 0 + scatt_eouts[-1].c = np.cumsum(energy[g, :]) - energy[g, 0] + scatt_eouts.append(scatt_eouts[-1]) + + # Build the angular distributions associated with each outgoing energy + # group + mu = np.linspace(-1., 1., tabular.num_orders) + scatt_angles = [] + for gin in range(self.energy_groups.num_groups): + scatt_angles.append([]) + for gout in range(self.energy_groups.num_groups): + scatt_angles[gin].append( + openmc.stats.Tabular(mu, distrib[gin, gout, :])) + # Ensure the distribution CDF starts with 0 + scatt_angles[gin][-1].c = \ + np.cumsum(distrib[gin, gout, :]) - distrib[gin, gout, 0] + scatt_angles.append(scatt_angles[-1]) + + # Combine the energy and angle distributions in to a correlated + # angle/energy object + scatt_prod.distribution = \ + [openmc.data.CorrelatedAngleEnergy( + breakpoints=[len(scatt_eouts)], interpolation=[1], + energy=self.energy_groups.group_edges[:-1], + energy_out=scatt_eouts, mu=scatt_angles)] + + # Finally build the reaction with the just-calculated information + # This will be set to the (n,2n) reaction, though any scattering + # reaction would suffice + scatt_rxn = openmc.data.Reaction(16) + scatt_rxn.xs[strT] = \ + convert_xs(self.energy_groups.group_edges, scatt_xs[::-1]) + scatt_rxn.products = [scatt_prod] + scatt_rxn.center_of_mass = False + data.reactions[16] = scatt_rxn + + return data + def to_hdf5(self, file): """Write XSdata to an HDF5 file @@ -2512,6 +2714,29 @@ class MGXSLibrary(object): return library + def convert_to_continuous_energy(self, h5_filename='ce_mgxs.h5', + library_filename='ce_mgxs.xml'): + """Converts the MGXSLibrary object to an equivalent + library of openmc.data.IncidentNeutron objects + + Parameters + ---------- + h5_filename : str + HDF5 file to write with all the files + library_filename : str + cross_sections.xml file describing the HDF5 file + + """ + + library = openmc.data.DataLibrary() + data = [] + for i, xsdata in enumerate(self.xsdatas): + data.append(xsdata.convert_to_continuous_energy()) + data[-1].export_to_hdf5(h5_filename) + + library.register_file(h5_filename) + library.export_to_xml(library_filename) + def export_to_hdf5(self, filename='mgxs.h5'): """Create an hdf5 file that can be used for a simulation. From 2987ed88221663400dba872878b5178118021021 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 12 Feb 2017 13:18:52 -0500 Subject: [PATCH 05/29] Fix for floating point errors (rarely) causing crashes as no outgoing group is found --- openmc/mgxs_library.py | 6 ------ src/scattdata_header.F90 | 6 +++--- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 6f46f4887..b2dee3948 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1849,12 +1849,6 @@ class XSdata(object): orig_data[..., l] for l in range(self.num_orders)) new_data[..., h_bin] = simps(table_fine, mu_fine) - # new_data[..., h_bin] = \ - # simps(np.sum((l + 0.5) * - # eval_legendre(l, mu_fine) * - # orig_data[..., l] - # for l in range(self.num_orders)), - # mu_fine) # Remove the very small results from numerical precision # issues (allowing conversions to be reproduced exactly) diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index 3d9df0bbf..511e2a237 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -522,7 +522,7 @@ contains gout = this % gmin(gin) prob = this % energy(gin) % data(gout) - do while (prob < xi) + do while ((prob < xi) .and. (gout < this % gmax(gin))) gout = gout + 1 prob = prob + this % energy(gin) % data(gout) end do @@ -568,7 +568,7 @@ contains gout = this % gmin(gin) prob = this % energy(gin) % data(gout) - do while (prob < xi) + do while ((prob < xi) .and. (gout < this % gmax(gin))) gout = gout + 1 prob = prob + this % energy(gin) % data(gout) end do @@ -605,7 +605,7 @@ contains gout = this % gmin(gin) prob = this % energy(gin) % data(gout) - do while (prob < xi) + do while ((prob < xi) .and. (gout < this % gmax(gin))) gout = gout + 1 prob = prob + this % energy(gin) % data(gout) end do From 8093f3654de5decaf4f90ec18a429a49f2bfcf64 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 12 Feb 2017 20:36:39 -0500 Subject: [PATCH 06/29] Tested scattering rigorously which resulted in some changes --- openmc/mgxs_library.py | 73 +++++++++++++++++++++--------------------- 1 file changed, 37 insertions(+), 36 deletions(-) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index b2dee3948..14fb3dddc 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1865,7 +1865,7 @@ class XSdata(object): eval_legendre(l, mu_fine) for l in range(xsdata.num_orders)] for l in range(xsdata.num_orders): - new_data[..., l] = (l + 0.5) * simps(y[l], mu_fine) + new_data[..., l] = simps(y[l], mu_fine) # Remove the very small results from numerical precision # issues (allowing conversions to be reproduced exactly) @@ -1893,38 +1893,30 @@ class XSdata(object): new_data[..., np.abs(new_data) < 1.E-10] = 0. elif self.scatter_format == 'histogram': # The histogram format does not have enough information to - # convert to the other forms without inducing an error. - # We will make the assumption that the left-edge of the bin - # has the value of the bin. The mu=1 value will be extrapolated - # from the previous two points. - mu_self = np.linspace(-1, 1, self.num_orders + 1) - tab_data = np.zeros(orig_data.shape[:-1] + - (orig_data.shape[-1] + 1,)) - tab_data[..., :-1] = orig_data - tab_data[..., -1] = (orig_data[..., -2] - - orig_data[..., -3]) + orig_data[..., -2] - # Now get the distribution normalization factor. - # This factor will be such that its average value is the - # group -> group cross section (which is the sum of the - # original data) - norm = np.sum(orig_data, axis=-1) / np.mean(tab_data, axis=-1) - norm = np.nan_to_num(norm) - for i in range(tab_data.shape[-1]): - tab_data[..., i] *= norm[...] + # convert to the other forms without inducing some amount of + # error. We will make the assumption that the center of the bin + # has the value of the bin. The mu=-1 and 1 points will be + # extrapolated from the shape. + mu_midpoint = np.linspace(-1, 1, self.num_orders, + endpoint=False) + mu_midpoint += (mu_midpoint[1] - mu_midpoint[0]) * 0.5 + interp = interp1d(mu_midpoint, orig_data, + fill_value='extrapolate') + # Now get the distribution normalization factor to take from + # an integral quantity to a point-wise quantity + norm = float(self.num_orders) / 2.0 - # We now have a tabular distribution in tab_data on mu_self - # with a normalization in norm. We now proceed just like the - # tabular branch above. + # We now have a tabular distribution in tab_data on mu_self. + # We now proceed just like the tabular branch above. if target_format == 'legendre': # find the legendre coefficients via integration. To best # use the vectorized integration capabilities of scipy, # this will be done with fixed sample integration routines. mu_fine = np.linspace(-1, 1, _NMU) - y = [interp1d(mu_self, tab_data)(mu_fine) * - eval_legendre(l, mu_fine) + y = [interp(mu_fine) * norm * eval_legendre(l, mu_fine) for l in range(xsdata.num_orders)] for l in range(xsdata.num_orders): - new_data[..., l] = (l + 0.5) * simps(y[l], mu_fine) + new_data[..., l] = simps(y[l], mu_fine) # Remove the very small results from numerical precision # issues (allowing conversions to be reproduced exactly) @@ -1932,7 +1924,7 @@ class XSdata(object): elif target_format == 'tabular': # Simply use an interpolating function to get the new data mu = np.linspace(-1, 1, xsdata.num_orders) - new_data[..., :] = interp1d(mu_self, tab_data)(mu) + new_data[..., :] = interp(mu) * norm elif target_format == 'histogram': # Use an interpolating function to do the bin-wise # integrals @@ -1942,15 +1934,14 @@ class XSdata(object): # be written to utilize the vectorized integration # capabilities instead of having an isotropic and # angle representation path. - interp = interp1d(mu_self, tab_data) for h_bin in range(xsdata.num_orders): mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU) - new_data[..., h_bin] = simps(interp(mu_fine), mu_fine) + new_data[..., h_bin] = \ + norm * simps(interp(mu_fine), mu_fine) # Remove the very small results from numerical precision # issues (allowing conversions to be reproduced exactly) new_data[..., np.abs(new_data) < 1.E-10] = 0. - xsdata.set_scatter_matrix(new_data, temp) return xsdata @@ -2107,10 +2098,18 @@ class XSdata(object): # lines within each of the outgoing groups scatt_eouts = [] for g in range(self.energy_groups.num_groups): + nz = np.nonzero(energy[g, :]) + lo = nz[0][0] + hi = nz[0][-1] + 1 + # scatt_eouts.append( + # openmc.stats.Discrete(energy_midpoints[lo: hi], + # energy[g, lo:hi])) scatt_eouts.append( - openmc.stats.Discrete(energy_midpoints, energy[g, :])) + openmc.stats.Tabular(self.energy_groups.group_edges[:-1], + energy[g, :], + interpolation='histogram')) # Ensure the distribution CDF starts with 0 - scatt_eouts[-1].c = np.cumsum(energy[g, :]) - energy[g, 0] + scatt_eouts[-1].c = np.cumsum(energy[g, lo:hi]) - energy[g, lo] scatt_eouts.append(scatt_eouts[-1]) # Build the angular distributions associated with each outgoing energy @@ -2120,11 +2119,13 @@ class XSdata(object): for gin in range(self.energy_groups.num_groups): scatt_angles.append([]) for gout in range(self.energy_groups.num_groups): - scatt_angles[gin].append( - openmc.stats.Tabular(mu, distrib[gin, gout, :])) - # Ensure the distribution CDF starts with 0 - scatt_angles[gin][-1].c = \ - np.cumsum(distrib[gin, gout, :]) - distrib[gin, gout, 0] + if energy[gin, gout] > 0.: + scatt_angles[gin].append( + openmc.stats.Tabular(mu, distrib[gin, gout, :])) + # Ensure the distribution CDF starts with 0 + scatt_angles[gin][-1].c = \ + np.cumsum(distrib[gin, gout, :]) - \ + distrib[gin, gout, 0] scatt_angles.append(scatt_angles[-1]) # Combine the energy and angle distributions in to a correlated From 13bbabd8cfa388593ea51ee2ab256480b66b8d22 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 13 Feb 2017 19:51:25 -0500 Subject: [PATCH 07/29] old code cleanup --- openmc/mgxs_library.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 14fb3dddc..b4482d60b 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -2098,18 +2098,12 @@ class XSdata(object): # lines within each of the outgoing groups scatt_eouts = [] for g in range(self.energy_groups.num_groups): - nz = np.nonzero(energy[g, :]) - lo = nz[0][0] - hi = nz[0][-1] + 1 - # scatt_eouts.append( - # openmc.stats.Discrete(energy_midpoints[lo: hi], - # energy[g, lo:hi])) scatt_eouts.append( openmc.stats.Tabular(self.energy_groups.group_edges[:-1], energy[g, :], interpolation='histogram')) # Ensure the distribution CDF starts with 0 - scatt_eouts[-1].c = np.cumsum(energy[g, lo:hi]) - energy[g, lo] + scatt_eouts[-1].c = np.cumsum(energy[g, :]) - energy[g, 0] scatt_eouts.append(scatt_eouts[-1]) # Build the angular distributions associated with each outgoing energy From 130db90c2a74f5a24e4248ef2d0bed8f7b81e4fa Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Wed, 15 Feb 2017 20:27:31 -0500 Subject: [PATCH 08/29] cleaned up comments and removed the conversion to CE --- openmc/mgxs_library.py | 268 ++++++----------------------------------- 1 file changed, 36 insertions(+), 232 deletions(-) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index b4482d60b..e3a8ad4d8 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1721,22 +1721,19 @@ class XSdata(object): check_greater_than('num_azimuthal', num_azimuthal, 0) xsdata = copy.deepcopy(self) + + # First handle the case where the current and requested + # representations are the same if target_representation == self.representation: # Check to make sure the num_polar and num_azimuthal values match if target_representation == 'angle': if num_polar != self.num_polar or num_azimuthal != self.num_azimuthal: - raise NotImplementedError("XSdata.convert_representation " - "cannot translate between " + raise NotImplementedError("XCannot translate between " "`angle` representations with " "different angle bin structures") - # Nothing to do + # Nothing to do as the same structure was requested return xsdata - types = ['total', 'absorption', 'fission', 'nu_fission', - 'scatter_matrix', 'multiplicity_matrix', 'prompt_nu_fission', - 'delayed_nu_fission', 'kappa_fission', 'chi', 'chi_prompt', - 'chi_delayed', 'beta', 'decay_rate', 'inverse_velocity'] - xsdata.representation = target_representation # We have different actions depending on the representation conversion if target_representation == 'isotropic': @@ -1748,17 +1745,21 @@ class XSdata(object): xsdata.num_polar = num_polar xsdata.num_azimuthal = num_azimuthal - # Reset XSdata.xs_shapes to accomodate the new shape + # Reset xs_shapes so it is recalculated the next time it is needed xsdata._xs_shapes = None for i, temp in enumerate(xsdata.temperatures): - for xs in types: + for xs in ['total', 'absorption', 'fission', 'nu_fission', + 'scatter_matrix', 'multiplicity_matrix', + 'prompt_nu_fission', 'delayed_nu_fission', + 'kappa_fission', 'chi', 'chi_prompt', 'chi_delayed', + 'beta', 'decay_rate', 'inverse_velocity']: # Get the original data orig_data = getattr(self, '_' + xs)[i] if orig_data is not None: if target_representation == 'isotropic': # Since we are going from angle to isotropic, the - # current data is just averaged over the angle bins + # current data is just the average over the angle bins new_data = orig_data.mean(axis=(0, 1)) elif target_representation == 'angle': # Since we are going from isotropic to angle, the @@ -1787,7 +1788,7 @@ class XSdata(object): Returns ------- openmc.XSdata - Multi-group cross section data with the same data as self, but + Multi-group cross section data with the same data as in self, but represented as specified in :param:`target_format`. """ @@ -1806,15 +1807,11 @@ class XSdata(object): xsdata = copy.deepcopy(self) xsdata.scatter_format = target_format xsdata.order = target_order - # Reset XSdata.xs_shapes to accomodate the new shape + + # Reset and re-generate XSdata.xs_shapes with the new scattering format xsdata._xs_shapes = None xsdata.xs_shapes - # We have to accomodate the following possibilities: - # histogram -> tabular w/ same or diff order - # histogram -> legendre - # histogram -> histogram w/ same or diff order - for i, temp in enumerate(xsdata.temperatures): orig_data = self._scatter_matrix[i] new_shape = orig_data.shape[:-1] + (xsdata.num_orders,) @@ -1827,21 +1824,23 @@ class XSdata(object): new_data[..., :order] = orig_data[..., :order] elif target_format == 'tabular': mu = np.linspace(-1, 1, xsdata.num_orders) - # Evaluate the legendre on the tabular grid within mu + # Evaluate the legendre on the mu grid for imu in range(len(mu)): new_data[..., imu] = \ np.sum((l + 0.5) * eval_legendre(l, mu[imu]) * orig_data[..., l] for l in range(self.num_orders)) elif target_format == 'histogram': + # This code uses the vectorized integration capabilities + # instead of having an isotropic and angle representation + # path. + # Set the histogram mu grid mu = np.linspace(-1, 1, xsdata.num_orders + 1) - # This code will be written to utilize the vectorized - # integration capabilities instead of having an isotropic - # and angle representation path. + # For every bin perform simpson integration of a finely + # sampled orig_data for h_bin in range(xsdata.num_orders): mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU) - table_shape = new_data.shape[:-1] + (_NMU,) - table_fine = np.zeros(table_shape) + table_fine = np.zeros(new_data.shape[:-1] + (_NMU,)) for imu in range(len(mu_fine)): table_fine[..., imu] = \ np.sum((l + 0.5) * @@ -1850,16 +1849,17 @@ class XSdata(object): for l in range(self.num_orders)) new_data[..., h_bin] = simps(table_fine, mu_fine) - # Remove the very small results from numerical precision - # issues (allowing conversions to be reproduced exactly) + # Remove the very small values resulting from numerical + # precision issues new_data[..., np.abs(new_data) < 1.E-10] = 0. elif self.scatter_format == 'tabular': + # Calculate the mu points of the current data mu_self = np.linspace(-1, 1, self.num_orders) if target_format == 'legendre': - # find the legendre coefficients via integration. To best + # Find the Legendre coefficients via integration. To best # use the vectorized integration capabilities of scipy, - # this will be done with fixed sample integration routines. + # this is done with fixed sample integration routines. mu_fine = np.linspace(-1, 1, _NMU) y = [interp1d(mu_self, orig_data)(mu_fine) * eval_legendre(l, mu_fine) @@ -1867,8 +1867,8 @@ class XSdata(object): for l in range(xsdata.num_orders): new_data[..., l] = simps(y[l], mu_fine) - # Remove the very small results from numerical precision - # issues (allowing conversions to be reproduced exactly) + # Remove the very small values resulting from numerical + # precision issues new_data[..., np.abs(new_data) < 1.E-10] = 0. elif target_format == 'tabular': # Simply use an interpolating function to get the new data @@ -1888,8 +1888,8 @@ class XSdata(object): mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU) new_data[..., h_bin] = simps(interp(mu_fine), mu_fine) - # Remove the very small results from numerical precision - # issues (allowing conversions to be reproduced exactly) + # Remove the very small values resulting from numerical + # precision issues new_data[..., np.abs(new_data) < 1.E-10] = 0. elif self.scatter_format == 'histogram': # The histogram format does not have enough information to @@ -1918,8 +1918,8 @@ class XSdata(object): for l in range(xsdata.num_orders): new_data[..., l] = simps(y[l], mu_fine) - # Remove the very small results from numerical precision - # issues (allowing conversions to be reproduced exactly) + # Remove the very small values resulting from numerical + # precision issues new_data[..., np.abs(new_data) < 1.E-10] = 0. elif target_format == 'tabular': # Simply use an interpolating function to get the new data @@ -1939,209 +1939,13 @@ class XSdata(object): new_data[..., h_bin] = \ norm * simps(interp(mu_fine), mu_fine) - # Remove the very small results from numerical precision - # issues (allowing conversions to be reproduced exactly) + # Remove the very small values resulting from numerical + # precision issues new_data[..., np.abs(new_data) < 1.E-10] = 0. xsdata.set_scatter_matrix(new_data, temp) return xsdata - def convert_to_continuous_energy(self, temperature=294.): - """Converts the XSdata object to an equivalent - openmc.data.IncidentNeutron object. - - Parameters - ---------- - temperature : float - Temperature of dataset to print; defaults to 294K - - Returns - ------- - openmc.data.IncidentNeutron - The continuous-energy IncidentNeutron data library equivalent - to the MGXS data in self - """ - - # Check if this can be performed successfully - if self.representation == 'angle': - raise ValueError("Cannot convert angle-dependent MGXS; convert to " - "an isotropic representation first") - required_types = ['absorption', 'scatter_matrix'] - for type_check in required_types: - if getattr(self, '_' + type_check)[0] is None: - raise ValueError(type_check + ' data is required') - if temperature not in self.temperatures: - raise ValueError("Invalid temperature") - - def convert_xs(group_edges, values): - cexs = openmc.data.Tabulated1D(group_edges, - np.append(values[::-1], - [values[0]]), - breakpoints=[len(group_edges)], - interpolation=[1]) - return cexs - - # Build required metadata - kTs = self.temperatures[:] * openmc.data.K_BOLTZMANN - if self.atomic_weight_ratio: - awr = self.atomic_weight_ratio - else: - awr = 1. - - # Get the temperature index - iT = np.where(self.temperatures == temperature)[0][0] - data = openmc.data.IncidentNeutron(self.name, 1, 1, 0, awr, kTs) - strT = "{}K".format(int(round(self.temperatures[iT]))) - data.energy = {strT: self.energy_groups.group_edges} - energy_midpoints = (self.energy_groups.group_edges[1:] + - self.energy_groups.group_edges[:-1]) / 2. - - # Elastic MGXS: must explicitly be 0 barns to avoid incorrect - # elastic_scatter calculation with data available to us in MG library. - el_rxn = openmc.data.Reaction(2) - el_rxn.xs[strT] = \ - convert_xs(self.energy_groups.group_edges, - np.zeros(self.energy_groups.num_groups)) - data.reactions[2] = el_rxn - - # Absorption MGXS - abs_rxn = openmc.data.Reaction(102) - abs_rxn.xs[strT] = convert_xs(self.energy_groups.group_edges, - np.subtract(self._absorption[iT], - self._fission[iT])) - data.reactions[102] = abs_rxn - - if self._nu_fission[iT] is not None and self._fission[iT] is not None: - fiss_rxn = openmc.data.Reaction(18) - fiss_rxn.xs[strT] = convert_xs(self.energy_groups.group_edges, - self._fission[iT]) - - # Get nu_fission and chi from the presence of both, or the - # nu_fission matrix - if self._nu_fission[iT].shape == self.xs_shapes["[G][G']"]: - nu_fiss = np.sum(self._nu_fission[iT], axis=1)[::-1] - chi = self._nu_fission[iT][::-1] / nu_fiss - else: - nu_fiss = self._nu_fission[iT][::-1] - chi = np.reshape(np.tile(self._chi[iT][::-1], - self.energy_groups.num_groups), - (self.energy_groups.num_groups, - self.energy_groups.num_groups)) - nu = np.divide(nu_fiss, self._fission[iT][::-1]) - - # Build a histogram distribution to represent the yield for the - # incoming groups - prod = openmc.data.Product() - prod.yield_ = \ - openmc.data.Tabulated1D(self.energy_groups.group_edges[:-1], - nu, breakpoints=[len(nu)], - interpolation=[1]) - - # Now build the outgoing energy distribution using discrete energy - # lines within each of the outgoing groups - chi_eouts = [] - for g in range(self.energy_groups.num_groups): - chi_eouts.append(openmc.stats.Discrete(energy_midpoints, - chi[g, :])) - # Ensure the distribution CDF starts with 0 - chi_eouts[-1].c = np.cumsum(chi[g, :]) - chi[g, 0] - - # Create the continuous distribution for the fission - # The angular distribution will remain None (thus isotropic) - chi_distrib = openmc.data.ContinuousTabular( - [len(chi)], [1], self.energy_groups.group_edges[:-1], - chi_eouts) - prod.distribution = \ - [openmc.data.UncorrelatedAngleEnergy(energy=chi_distrib)] - fiss_rxn.products = [prod] - fiss_rxn.center_of_mass = False - data.reactions[18] = fiss_rxn - - # Scattering Data - # First convert the data to a tabular representation - if self.scatter_format == 'tabular': - tabular = self - else: - tabular = self.convert_scatter_format('tabular', 33) - - # Calculate the isotropic scattering matrix and use that to find the - # total scattering x/s and outgoing energy distributions - isotropic_matrix = np.mean(tabular._scatter_matrix[iT], axis=-1)[::-1, - ::-1] - scatt_xs = np.sum(isotropic_matrix, axis=1) - energy = np.zeros((self.energy_groups.num_groups, - self.energy_groups.num_groups)) - for gin in range(self.energy_groups.num_groups): - energy[gin, :] = isotropic_matrix[gin, :] / scatt_xs[gin] - - # Get the anisotropic but normalized angular distribution - distrib = np.zeros((self.energy_groups.num_groups, - self.energy_groups.num_groups, - tabular.num_orders)) - for gin in range(self.energy_groups.num_groups): - for gout in range(self.energy_groups.num_groups): - distrib[gin, gout, :] = \ - np.divide(tabular._scatter_matrix[iT][gin, gout, :], - isotropic_matrix[gin, gout]) - distrib = np.nan_to_num(distrib) - - # Incorporate the scattering multiplication, if required - scatt_prod = openmc.data.Product() - if self._multiplicity_matrix[iT]: - yield_ = self._multiplicity_matrix[iT][::-1, ::-1] - scatt_prod.yield_ = \ - openmc.data.Tabulated1D(self.energy_groups.group_edges[:-1], - yield_, breakpoints=[len(yield_)], - interpolation=[1]) - - # Now build the outgoing energy distribution using discrete energy - # lines within each of the outgoing groups - scatt_eouts = [] - for g in range(self.energy_groups.num_groups): - scatt_eouts.append( - openmc.stats.Tabular(self.energy_groups.group_edges[:-1], - energy[g, :], - interpolation='histogram')) - # Ensure the distribution CDF starts with 0 - scatt_eouts[-1].c = np.cumsum(energy[g, :]) - energy[g, 0] - scatt_eouts.append(scatt_eouts[-1]) - - # Build the angular distributions associated with each outgoing energy - # group - mu = np.linspace(-1., 1., tabular.num_orders) - scatt_angles = [] - for gin in range(self.energy_groups.num_groups): - scatt_angles.append([]) - for gout in range(self.energy_groups.num_groups): - if energy[gin, gout] > 0.: - scatt_angles[gin].append( - openmc.stats.Tabular(mu, distrib[gin, gout, :])) - # Ensure the distribution CDF starts with 0 - scatt_angles[gin][-1].c = \ - np.cumsum(distrib[gin, gout, :]) - \ - distrib[gin, gout, 0] - scatt_angles.append(scatt_angles[-1]) - - # Combine the energy and angle distributions in to a correlated - # angle/energy object - scatt_prod.distribution = \ - [openmc.data.CorrelatedAngleEnergy( - breakpoints=[len(scatt_eouts)], interpolation=[1], - energy=self.energy_groups.group_edges[:-1], - energy_out=scatt_eouts, mu=scatt_angles)] - - # Finally build the reaction with the just-calculated information - # This will be set to the (n,2n) reaction, though any scattering - # reaction would suffice - scatt_rxn = openmc.data.Reaction(16) - scatt_rxn.xs[strT] = \ - convert_xs(self.energy_groups.group_edges, scatt_xs[::-1]) - scatt_rxn.products = [scatt_prod] - scatt_rxn.center_of_mass = False - data.reactions[16] = scatt_rxn - - return data - def to_hdf5(self, file): """Write XSdata to an HDF5 file From 855c1ca19bcecb82e9c7b5c9a497a7236c1f3ca8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 16 Jan 2017 19:58:15 -0600 Subject: [PATCH 09/29] Build OpenMC as a library (all except main.F90) --- CMakeLists.txt | 100 +++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 96 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 64eee2e9f..098de3de2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -266,15 +266,105 @@ add_library(faddeeva STATIC src/Faddeeva.c) #=============================================================================== set(program "openmc") -file(GLOB source src/*.F90 src/xml/openmc_fox.F90) -add_executable(${program} ${source}) +set(LIBOPENMC_FORTRAN_SRC + src/algorithm.F90 + src/angle_distribution.F90 + src/angleenergy_header.F90 + src/bank_header.F90 + src/cmfd_data.F90 + src/cmfd_execute.F90 + src/cmfd_header.F90 + src/cmfd_input.F90 + src/cmfd_loss_operator.F90 + src/cmfd_prod_operator.F90 + src/cmfd_solver.F90 + src/constants.F90 + src/cross_section.F90 + src/dict_header.F90 + src/distribution_multivariate.F90 + src/distribution_univariate.F90 + src/doppler.F90 + src/eigenvalue.F90 + src/endf.F90 + src/endf_header.F90 + src/energy_distribution.F90 + src/energy_grid.F90 + src/error.F90 + src/finalize.F90 + src/geometry.F90 + src/geometry_header.F90 + src/global.F90 + src/hdf5_interface.F90 + src/initialize.F90 + src/input_xml.F90 + src/list_header.F90 + src/material_header.F90 + src/math.F90 + src/matrix_header.F90 + src/mesh.F90 + src/mesh_header.F90 + src/message_passing.F90 + src/mgxs_data.F90 + src/mgxs_header.F90 + src/multipole.F90 + src/multipole_header.F90 + src/nuclide_header.F90 + src/output.F90 + src/particle_header.F90 + src/particle_restart.F90 + src/particle_restart_write.F90 + src/physics_common.F90 + src/physics.F90 + src/physics_mg.F90 + src/plot.F90 + src/plot_header.F90 + src/ppmlib.F90 + src/product_header.F90 + src/progress_header.F90 + src/random_lcg.F90 + src/reaction_header.F90 + src/relaxng + src/sab_header.F90 + src/scattdata_header.F90 + src/secondary_correlated.F90 + src/secondary_kalbach.F90 + src/secondary_nbody.F90 + src/secondary_uncorrelated.F90 + src/set_header.F90 + src/simulation.F90 + src/source.F90 + src/source_header.F90 + src/state_point.F90 + src/stl_vector.F90 + src/string.F90 + src/summary.F90 + src/surface_header.F90 + src/tally.F90 + src/tally_filter.F90 + src/tally_filter_header.F90 + src/tally_header.F90 + src/tally_initialize.F90 + src/timer_header.F90 + src/tracking.F90 + src/track_output.F90 + src/trigger.F90 + src/trigger_header.F90 + src/urr_header.F90 + src/vector_header.F90 + src/volume_calc.F90 + src/volume_header.F90 + src/xml_interface.F90 + src/xml/openmc_fox.F90) +add_library(libopenmc STATIC ${LIBOPENMC_FORTRAN_SRC}) +set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc) +add_executable(${program} src/main.F90) # target_include_directories was added in CMake 2.8.11 and is the recommended # way to set include directories. For lesser versions, we revert to set_property if(CMAKE_VERSION VERSION_LESS 2.8.11) include_directories(${HDF5_INCLUDE_DIRS}) else() - target_include_directories(${program} PUBLIC ${HDF5_INCLUDE_DIRS}) + target_include_directories(libopenmc PUBLIC ${HDF5_INCLUDE_DIRS}) endif() # target_compile_options was added in CMake 2.8.12 and is the recommended way to @@ -288,6 +378,7 @@ if (CMAKE_VERSION VERSION_LESS 2.8.12) set_property(TARGET faddeeva PROPERTY COMPILE_FLAGS "${cflags}") else() target_compile_options(${program} PUBLIC ${f90flags}) + target_compile_options(libopenmc PUBLIC ${f90flags}) target_compile_options(faddeeva PRIVATE ${cflags}) endif() @@ -298,7 +389,8 @@ endforeach() # target_link_libraries treats any arguments starting with - but not -l as # linker flags. Thus, we can pass both linker flags and libraries together. -target_link_libraries(${program} ${ldflags} ${HDF5_LIBRARIES} fox_dom faddeeva) +target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} fox_dom faddeeva) +target_link_libraries(${program} ${ldflags} libopenmc) #=============================================================================== # Install executable, scripts, manpage, license From 50220c0ef1bc4aafc122576ceff2e00b65bb655b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 18 Jan 2017 07:01:30 -0600 Subject: [PATCH 10/29] Pass MPI intracommunicator to initialize subroutine --- src/cmfd_execute.F90 | 22 +++++++------------- src/cmfd_input.F90 | 1 + src/cmfd_solver.F90 | 6 ++++-- src/eigenvalue.F90 | 20 +++++++++--------- src/error.F90 | 4 +--- src/finalize.F90 | 19 ++++++++--------- src/global.F90 | 15 -------------- src/hdf5_interface.F90 | 10 ++++----- src/initialize.F90 | 45 ++++++++++++++++++++++++----------------- src/input_xml.F90 | 1 + src/main.F90 | 15 +++++++++----- src/mesh.F90 | 11 ++++------ src/message_passing.F90 | 17 ++++++++++++++++ src/output.F90 | 1 + src/physics.F90 | 1 + src/physics_mg.F90 | 1 + src/simulation.F90 | 7 ++----- src/source.F90 | 1 + src/state_point.F90 | 30 ++++++++++++--------------- src/summary.F90 | 1 + src/tally.F90 | 17 +++++++--------- src/tracking.F90 | 1 + src/trigger.F90 | 1 + src/volume_calc.F90 | 8 ++++---- 24 files changed, 126 insertions(+), 129 deletions(-) diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index a5d92002b..b1f0041d3 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -22,6 +22,7 @@ contains use cmfd_data, only: set_up_cmfd use cmfd_solver, only: cmfd_solver_execute use error, only: warning, fatal_error + use message_passing, only: master ! CMFD single processor on master if (master) then @@ -90,13 +91,9 @@ contains subroutine calc_fission_source() use constants, only: CMFD_NOACCEL, ZERO, TWO - use global, only: cmfd, cmfd_coremap, master, entropy_on, current_batch - use string, only: to_str - -#ifdef MPI - use global, only: mpi_err + use global, only: cmfd, cmfd_coremap, entropy_on, current_batch use message_passing -#endif + use string, only: to_str integer :: nx ! maximum number of cells in x direction integer :: ny ! maximum number of cells in y direction @@ -202,7 +199,7 @@ contains #ifdef MPI ! Broadcast full source to all procs - call MPI_BCAST(cmfd % cmfd_src, n, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) + call MPI_BCAST(cmfd % cmfd_src, n, MPI_REAL8, 0, mpi_intracomm, mpi_err) #endif end subroutine calc_fission_source @@ -216,16 +213,11 @@ contains use algorithm, only: binary_search use constants, only: ZERO, ONE use error, only: warning, fatal_error - use global, only: meshes, source_bank, work, n_user_meshes, cmfd, & - master + use global, only: meshes, source_bank, work, n_user_meshes, cmfd use mesh_header, only: RegularMesh use mesh, only: count_bank_sites, get_mesh_indices - use string, only: to_str - -#ifdef MPI - use global, only: mpi_err use message_passing -#endif + use string, only: to_str logical, intent(in) :: new_weights ! calcualte new weights @@ -289,7 +281,7 @@ contains ! Broadcast weight factors to all procs #ifdef MPI call MPI_BCAST(cmfd % weightfactors, ng*nx*ny*nz, MPI_REAL8, 0, & - MPI_COMM_WORLD, mpi_err) + mpi_intracomm, mpi_err) #endif end if diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index f10762750..95e1b1811 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -15,6 +15,7 @@ contains subroutine configure_cmfd() use cmfd_header, only: allocate_cmfd + use message_passing, only: master integer :: color ! color group of processor diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 84f0016d3..fb687d10b 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -301,10 +301,12 @@ contains subroutine convergence(iter, innerits, iconv) - use constants, only: ONE, ZERO - use global, only: cmfd_power_monitor, master use, intrinsic :: ISO_FORTRAN_ENV + use constants, only: ONE, ZERO + use global, only: cmfd_power_monitor + use message_passing, only: master + integer, intent(in) :: iter ! outer iteration number integer, intent(in) :: innerits ! inner iteration nubmer logical, intent(out) :: iconv ! convergence logical diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 05c8190a6..a930f4410 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -1,8 +1,5 @@ module eigenvalue -#ifdef MPI - use message_passing -#endif use algorithm, only: binary_search use constants, only: ZERO @@ -11,6 +8,7 @@ module eigenvalue use math, only: t_percentile use mesh, only: count_bank_sites use mesh_header, only: RegularMesh + use message_passing use random_lcg, only: prn, set_particle_seed, advance_prn_seed use string, only: to_str @@ -66,7 +64,7 @@ contains #ifdef MPI start = 0_8 call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, & - MPI_COMM_WORLD, mpi_err) + mpi_intracomm, mpi_err) ! While we would expect the value of start on rank 0 to be 0, the MPI ! standard says that the receive buffer on rank 0 is undefined and not @@ -76,7 +74,7 @@ contains finish = start + n_bank total = finish call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & - MPI_COMM_WORLD, mpi_err) + mpi_intracomm, mpi_err) #else start = 0_8 @@ -150,13 +148,13 @@ contains ! First do an exclusive scan to get the starting indices for start = 0_8 call MPI_EXSCAN(index_temp, start, 1, MPI_INTEGER8, MPI_SUM, & - MPI_COMM_WORLD, mpi_err) + mpi_intracomm, mpi_err) finish = start + index_temp ! Allocate space for bank_position if this hasn't been done yet if (.not. allocated(bank_position)) allocate(bank_position(n_procs)) call MPI_ALLGATHER(start, 1, MPI_INTEGER8, bank_position, 1, & - MPI_INTEGER8, MPI_COMM_WORLD, mpi_err) + MPI_INTEGER8, mpi_intracomm, mpi_err) #else start = 0_8 finish = index_temp @@ -210,7 +208,7 @@ contains if (neighbor /= rank) then n_request = n_request + 1 call MPI_ISEND(temp_sites(index_local), int(n), MPI_BANK, neighbor, & - rank, MPI_COMM_WORLD, request(n_request), mpi_err) + rank, mpi_intracomm, request(n_request), mpi_err) end if ! Increment all indices @@ -254,7 +252,7 @@ contains n_request = n_request + 1 call MPI_IRECV(source_bank(index_local), int(n), MPI_BANK, & - neighbor, neighbor, MPI_COMM_WORLD, request(n_request), mpi_err) + neighbor, neighbor, mpi_intracomm, request(n_request), mpi_err) else ! If the source sites are on this procesor, we can simply copy them @@ -379,7 +377,7 @@ contains #ifdef MPI ! Combine values across all processors call MPI_ALLREDUCE(keff_generation, k_generation(overall_gen), 1, & - MPI_REAL8, MPI_SUM, MPI_COMM_WORLD, mpi_err) + MPI_REAL8, MPI_SUM, mpi_intracomm, mpi_err) #else k_generation(overall_gen) = keff_generation #endif @@ -565,7 +563,7 @@ contains #ifdef MPI ! Send source fraction to all processors n = product(ufs_mesh % dimension) - call MPI_BCAST(source_frac, n, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err) + call MPI_BCAST(source_frac, n, MPI_REAL8, 0, mpi_intracomm, mpi_err) #endif ! Normalize to total weight to get fraction of source in each cell diff --git a/src/error.F90 b/src/error.F90 index 9e5ca8479..1bb6b5257 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -3,9 +3,7 @@ module error use, intrinsic :: ISO_FORTRAN_ENV use constants -#ifdef MPI use message_passing -#endif implicit none @@ -139,7 +137,7 @@ contains #ifdef MPI ! Abort MPI - call MPI_ABORT(MPI_COMM_WORLD, code, mpi_err) + call MPI_ABORT(mpi_intracomm, code, mpi_err) #endif ! Abort program diff --git a/src/finalize.F90 b/src/finalize.F90 index 83a5ac5fa..7ad99d0b5 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -1,17 +1,14 @@ module finalize + use hdf5, only: h5tclose_f, h5close_f + use global + use hdf5_interface, only: hdf5_bank_t + use message_passing use output, only: print_runtime, print_results, & print_overlap_check, write_tallies use tally, only: tally_statistics -#ifdef MPI - use message_passing -#endif - - use hdf5_interface, only: hdf5_bank_t - use hdf5, only: h5tclose_f, h5close_f - implicit none contains @@ -21,7 +18,7 @@ contains ! statistics and writing out tallies !=============================================================================== - subroutine finalize_run() + subroutine openmc_finalize() integer :: hdf5_err @@ -66,7 +63,7 @@ contains call MPI_FINALIZE(mpi_err) #endif - end subroutine finalize_run + end subroutine openmc_finalize !=============================================================================== ! REDUCE_OVERLAP_COUNT accumulates cell overlap check counts to master @@ -77,10 +74,10 @@ contains #ifdef MPI if (master) then call MPI_REDUCE(MPI_IN_PLACE, overlap_check_cnt, n_cells, & - MPI_INTEGER8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_INTEGER8, MPI_SUM, 0, mpi_intracomm, mpi_err) else call MPI_REDUCE(overlap_check_cnt, overlap_check_cnt, n_cells, & - MPI_INTEGER8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_INTEGER8, MPI_SUM, 0, mpi_intracomm, mpi_err) end if #endif diff --git a/src/global.F90 b/src/global.F90 index 02f687018..c48367214 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -267,21 +267,6 @@ module global ! ============================================================================ ! PARALLEL PROCESSING VARIABLES - ! The defaults set here for the number of processors, rank, and master and - ! mpi_enabled flag are for when MPI is not being used at all, i.e. a serial - ! run. In this case, these variables are still used at times. - - integer :: n_procs = 1 ! number of processes - integer :: rank = 0 ! rank of process - logical :: master = .true. ! master process? - logical :: mpi_enabled = .false. ! is MPI in use and initialized? - integer :: mpi_err ! MPI error code -#ifdef MPIF08 - type(MPI_Datatype) :: MPI_BANK -#else - integer :: MPI_BANK ! MPI datatype for fission bank -#endif - #ifdef _OPENMP integer :: n_threads = NONE ! number of OpenMP threads integer :: thread_id ! ID of a given thread diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 5e14605ae..e588cda14 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -17,7 +17,7 @@ module hdf5_interface use error, only: fatal_error #ifdef PHDF5 - use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL + use message_passing, only: mpi_intracomm, MPI_INFO_NULL #endif implicit none @@ -124,10 +124,10 @@ contains call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err) #ifdef PHDF5 #ifdef MPIF08 - call h5pset_fapl_mpio_f(plist, MPI_COMM_WORLD%MPI_VAL, & + call h5pset_fapl_mpio_f(plist, mpi_intracomm%MPI_VAL, & MPI_INFO_NULL%MPI_VAL, hdf5_err) #else - call h5pset_fapl_mpio_f(plist, MPI_COMM_WORLD, MPI_INFO_NULL, hdf5_err) + call h5pset_fapl_mpio_f(plist, mpi_intracomm, MPI_INFO_NULL, hdf5_err) #endif #endif @@ -174,10 +174,10 @@ contains call h5pcreate_f(H5P_FILE_ACCESS_F, plist, hdf5_err) #ifdef PHDF5 #ifdef MPIF08 - call h5pset_fapl_mpio_f(plist, MPI_COMM_WORLD%MPI_VAL, & + call h5pset_fapl_mpio_f(plist, mpi_intracomm%MPI_VAL, & MPI_INFO_NULL%MPI_VAL, hdf5_err) #else - call h5pset_fapl_mpio_f(plist, MPI_COMM_WORLD, MPI_INFO_NULL, hdf5_err) + call h5pset_fapl_mpio_f(plist, mpi_intracomm, MPI_INFO_NULL, hdf5_err) #endif #endif diff --git a/src/initialize.F90 b/src/initialize.F90 index 692319acb..2101333b7 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1,5 +1,12 @@ module initialize + use, intrinsic :: ISO_C_BINDING, only: c_loc + + use hdf5 +#ifdef _OPENMP + use omp_lib +#endif + use bank_header, only: Bank use constants use dict_header, only: DictIntInt, ElemKeyValueII @@ -15,6 +22,7 @@ module initialize hdf5_integer8_t use input_xml, only: read_input_xml, cells_in_univ_dict, read_plots_xml use material_header, only: Material + use message_passing use mgxs_data, only: read_mgxs, create_macro_xs use output, only: title, header, print_version, write_message, & print_usage, print_plot @@ -27,30 +35,23 @@ module initialize use tally_filter use tally, only: init_tally_routines -#ifdef MPI - use message_passing -#endif - -#ifdef _OPENMP - use omp_lib -#endif - - use hdf5 - - use, intrinsic :: ISO_C_BINDING, only: c_loc - implicit none contains !=============================================================================== -! INITIALIZE_RUN takes care of all initialization tasks, i.e. reading +! OPENMC_INIT takes care of all initialization tasks, i.e. reading ! from command line, reading xml input files, initializing random ! number seeds, reading cross sections, initializing starting source, ! setting up timers, etc. !=============================================================================== - subroutine initialize_run() + subroutine openmc_init(intracomm) +#ifdef MPIF08 + type(MPI_Comm), intent(in) :: intracomm +#else + integer, intent(in), optional :: intracomm +#endif ! Start total and initialization timer call time_total%start() @@ -58,7 +59,7 @@ contains #ifdef MPI ! Setup MPI - call initialize_mpi() + call initialize_mpi(intracomm) #endif ! Initialize HDF5 interface @@ -155,7 +156,7 @@ contains ! Stop initialization timer call time_initialize%stop() - end subroutine initialize_run + end subroutine openmc_init #ifdef MPI !=============================================================================== @@ -164,7 +165,12 @@ contains ! each processor. !=============================================================================== - subroutine initialize_mpi() + subroutine initialize_mpi(intracomm) +#ifdef MPIF08 + type(MPI_Comm), intent(in) :: intracomm +#else + integer, intent(in) :: intracomm +#endif integer :: bank_blocks(5) ! Count for each datatype #ifdef MPIF08 @@ -182,8 +188,9 @@ contains call MPI_INIT(mpi_err) ! Determine number of processors and rank of each processor - call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, mpi_err) - call MPI_COMM_RANK(MPI_COMM_WORLD, rank, mpi_err) + mpi_intracomm = intracomm + call MPI_COMM_SIZE(mpi_intracomm, n_procs, mpi_err) + call MPI_COMM_RANK(mpi_intracomm, rank, mpi_err) ! Determine master if (rank == 0) then diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 8f25d6301..d96dd419a 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -15,6 +15,7 @@ module input_xml use hdf5_interface use list_header, only: ListChar, ListInt, ListReal use mesh_header, only: RegularMesh + use message_passing use mgxs_data, only: create_macro_xs, read_mgxs use multipole, only: multipole_read use output, only: write_message diff --git a/src/main.F90 b/src/main.F90 index aff1e2114..8582ef703 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -1,17 +1,22 @@ program main use constants - use finalize, only: finalize_run + use finalize, only: openmc_finalize use global - use initialize, only: initialize_run + use initialize, only: openmc_init + use message_passing use particle_restart, only: run_particle_restart use plot, only: run_plot use simulation, only: run_simulation implicit none - ! set up problem - call initialize_run() + ! Initialize run -- when run with MPI, pass communicator +#ifdef MPI + call openmc_init(MPI_COMM_WORLD) +#else + call openmc_init() +#endif ! start problem based on mode select case (run_mode) @@ -24,6 +29,6 @@ program main end select ! finalize run - call finalize_run() + call openmc_finalize() end program main diff --git a/src/mesh.F90 b/src/mesh.F90 index 0f342c050..b778cbdf0 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -1,13 +1,10 @@ module mesh -#ifdef MPI - use message_passing -#endif - use algorithm, only: binary_search use constants use global use mesh_header + use message_passing implicit none @@ -207,17 +204,17 @@ contains ! collect values from all processors if (master) then call MPI_REDUCE(MPI_IN_PLACE, cnt, n, MPI_REAL8, MPI_SUM, 0, & - MPI_COMM_WORLD, mpi_err) + mpi_intracomm, mpi_err) else ! Receive buffer not significant at other processors call MPI_REDUCE(cnt, dummy, n, MPI_REAL8, MPI_SUM, 0, & - MPI_COMM_WORLD, mpi_err) + mpi_intracomm, mpi_err) end if ! Check if there were sites outside the mesh for any processor if (present(sites_outside)) then call MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, & - MPI_COMM_WORLD, mpi_err) + mpi_intracomm, mpi_err) end if #else diff --git a/src/message_passing.F90 b/src/message_passing.F90 index 75d0bd301..6c13b00e8 100644 --- a/src/message_passing.F90 +++ b/src/message_passing.F90 @@ -8,4 +8,21 @@ module message_passing #endif #endif + ! The defaults set here for the number of processors, rank, and master and + ! mpi_enabled flag are for when MPI is not being used at all, i.e. a serial + ! run. In this case, these variables are still used at times. + + integer :: n_procs = 1 ! number of processes + integer :: rank = 0 ! rank of process + logical :: master = .true. ! master process? + logical :: mpi_enabled = .false. ! is MPI in use and initialized? + integer :: mpi_err ! MPI error code +#ifdef MPIF08 + type(MPI_Datatype) :: MPI_BANK ! MPI datatype for fission bank + type(MPI_Comm) :: mpi_intracomm ! MPI intra-communicator +#else + integer :: MPI_BANK ! MPI datatype for fission bank + integer :: mpi_intracomm ! MPI intra-communicator +#endif + end module message_passing diff --git a/src/output.F90 b/src/output.F90 index a0b4a99b5..c719dbed9 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -11,6 +11,7 @@ module output use math, only: t_percentile use mesh_header, only: RegularMesh use mesh, only: mesh_indices_to_bin, bin_to_mesh_indices + use message_passing, only: master, n_procs use nuclide_header use particle_header, only: LocalCoord, Particle use plot_header diff --git a/src/physics.F90 b/src/physics.F90 index d644b5a36..fceca20a9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -9,6 +9,7 @@ module physics use material_header, only: Material use math use mesh, only: get_mesh_indices + use message_passing use nuclide_header use output, only: write_message use particle_header, only: Particle diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 082b337a4..b5d34838c 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -9,6 +9,7 @@ module physics_mg use math, only: rotate_angle use mgxs_header, only: Mgxs, MgxsContainer use mesh, only: get_mesh_indices + use message_passing use output, only: write_message use particle_header, only: Particle use particle_restart_write, only: write_particle_restart diff --git a/src/simulation.F90 b/src/simulation.F90 index b4a3791f0..187632836 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -1,9 +1,5 @@ module simulation -#ifdef MPI - use message_passing -#endif - use cmfd_execute, only: cmfd_init_batch, execute_cmfd use constants, only: ZERO use eigenvalue, only: count_source_for_ufs, calculate_average_keff, & @@ -13,6 +9,7 @@ module simulation use eigenvalue, only: join_bank_from_threads #endif use global + use message_passing use output, only: write_message, header, print_columns, & print_batch_keff, print_generation use particle_header, only: Particle @@ -322,7 +319,7 @@ contains if (master) call check_triggers() #ifdef MPI call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, & - MPI_COMM_WORLD, mpi_err) + mpi_intracomm, mpi_err) #endif if (satisfy_triggers .or. & (trigger_on .and. current_batch == n_max_batches)) then diff --git a/src/source.F90 b/src/source.F90 index cc85bf6ec..4c0fe9891 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -15,6 +15,7 @@ module source use geometry_header, only: BASE_UNIVERSE use global use hdf5_interface, only: file_create, file_open, file_close, read_dataset + use message_passing, only: rank use output, only: write_message use particle_header, only: Particle use random_lcg, only: prn, set_particle_seed, prn_set_stream diff --git a/src/state_point.F90 b/src/state_point.F90 index bf0947163..257753b5b 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -11,26 +11,22 @@ module state_point ! intervals, using the tag. !=============================================================================== + use, intrinsic :: ISO_C_BINDING, only: c_loc, c_ptr + + use hdf5 use constants + use dict_header, only: ElemKeyValueII, ElemKeyValueCI use endf, only: reaction_name use error, only: fatal_error, warning use global use hdf5_interface + use mesh_header, only: RegularMesh + use message_passing use output, only: write_message, time_stamp + use random_lcg, only: seed use string, only: to_str, count_digits, zero_padded use tally_header, only: TallyObject - use mesh_header, only: RegularMesh - use dict_header, only: ElemKeyValueII, ElemKeyValueCI - use random_lcg, only: seed - -#ifdef MPI - use message_passing -#endif - - use hdf5 - - use, intrinsic :: ISO_C_BINDING, only: c_loc, c_ptr implicit none @@ -553,7 +549,7 @@ contains ! receive buffer without having a temporary variable #ifdef MPI call MPI_REDUCE(MPI_IN_PLACE, global_temp, n_bins, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, mpi_intracomm, mpi_err) #endif ! Transfer values to value on master @@ -567,7 +563,7 @@ contains ! Receive buffer not significant at other processors #ifdef MPI call MPI_REDUCE(global_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, mpi_intracomm, mpi_err) #endif end if @@ -616,7 +612,7 @@ contains ! a receive buffer without having a temporary variable #ifdef MPI call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_SUM, 0, mpi_intracomm, mpi_err) #endif ! At the end of the simulation, store the results back in the @@ -640,7 +636,7 @@ contains ! Receive buffer not significant at other processors #ifdef MPI call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, mpi_intracomm, mpi_err) #endif end if @@ -940,7 +936,7 @@ contains ! Receive source sites from other processes if (i > 0) then call MPI_RECV(source_bank, int(dims(1)), MPI_BANK, i, i, & - MPI_COMM_WORLD, MPI_STATUS_IGNORE, mpi_err) + mpi_intracomm, MPI_STATUS_IGNORE, mpi_err) end if #endif @@ -969,7 +965,7 @@ contains else #ifdef MPI call MPI_SEND(source_bank, int(work), MPI_BANK, 0, rank, & - MPI_COMM_WORLD, mpi_err) + mpi_intracomm, mpi_err) #endif end if diff --git a/src/summary.F90 b/src/summary.F90 index 97abeb206..efe07b1ab 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -8,6 +8,7 @@ module summary use hdf5_interface use material_header, only: Material use mesh_header, only: RegularMesh + use message_passing use nuclide_header use output, only: time_stamp use surface_header diff --git a/src/tally.F90 b/src/tally.F90 index 41987ec35..df1b5d2d6 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -2,10 +2,6 @@ module tally use, intrinsic :: ISO_C_BINDING -#ifdef MPI - use message_passing -#endif - use algorithm, only: binary_search use constants use cross_section, only: multipole_deriv_eval @@ -18,6 +14,7 @@ module tally mesh_intersects_1d, mesh_intersects_2d, & mesh_intersects_3d use mesh_header, only: RegularMesh + use message_passing use output, only: header use particle_header, only: LocalCoord, Particle use string, only: to_str @@ -4073,14 +4070,14 @@ contains ! The MPI_IN_PLACE specifier allows the master to copy values into ! a receive buffer without having a temporary variable call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_SUM, 0, mpi_intracomm, mpi_err) ! Transfer values to value on master t % results(RESULT_VALUE,:,:) = tally_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(tally_temp, dummy, n_bins, MPI_REAL8, & - MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_SUM, 0, mpi_intracomm, mpi_err) ! Reset value on other processors t % results(RESULT_VALUE,:,:) = ZERO @@ -4094,14 +4091,14 @@ contains if (master) then call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err) ! Transfer values back to global_tallies on master global_tallies(RESULT_VALUE, :) = global_temp else ! Receive buffer not significant at other processors call MPI_REDUCE(global_temp, dummy, N_GLOBAL_TALLIES, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err) ! Reset value on other processors global_tallies(RESULT_VALUE, :) = ZERO @@ -4111,11 +4108,11 @@ contains ! last realization if (master) then call MPI_REDUCE(MPI_IN_PLACE, total_weight, 1, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, mpi_intracomm, mpi_err) else ! Receive buffer not significant at other processors call MPI_REDUCE(total_weight, dummy, 1, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + 0, mpi_intracomm, mpi_err) end if end subroutine reduce_tally_results diff --git a/src/tracking.F90 b/src/tracking.F90 index a45f8e448..1b4f6b99b 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -8,6 +8,7 @@ module tracking use geometry_header, only: Universe, BASE_UNIVERSE use global use output, only: write_message + use message_passing use particle_header, only: LocalCoord, Particle use physics, only: collision use physics_mg, only: collision_mg diff --git a/src/trigger.F90 b/src/trigger.F90 index f41cc46cf..bf22327d6 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -10,6 +10,7 @@ module trigger use output, only: warning, write_message use mesh, only: bin_to_mesh_indices use mesh_header, only: RegularMesh + use message_passing, only: master use trigger_header, only: TriggerObject use tally, only: TallyObject use tally_filter, only: MeshFilter diff --git a/src/volume_calc.F90 b/src/volume_calc.F90 index dee646860..a9584d672 100644 --- a/src/volume_calc.F90 +++ b/src/volume_calc.F90 @@ -272,11 +272,11 @@ contains if (master) then #ifdef MPI do j = 1, n_procs - 1 - call MPI_RECV(n, 1, MPI_INTEGER, j, 0, MPI_COMM_WORLD, & + call MPI_RECV(n, 1, MPI_INTEGER, j, 0, mpi_intracomm, & MPI_STATUS_IGNORE, mpi_err) allocate(data(2*n)) - call MPI_RECV(data, 2*n, MPI_INTEGER, j, 1, MPI_COMM_WORLD, & + call MPI_RECV(data, 2*n, MPI_INTEGER, j, 1, mpi_intracomm, & MPI_STATUS_IGNORE, mpi_err) do k = 0, n - 1 do m = 1, master_indices(i_domain) % size() @@ -340,8 +340,8 @@ contains data(2*k + 2) = master_hits(i_domain) % data(k + 1) end do - call MPI_SEND(n, 1, MPI_INTEGER, 0, 0, MPI_COMM_WORLD, mpi_err) - call MPI_SEND(data, 2*n, MPI_INTEGER, 0, 1, MPI_COMM_WORLD, mpi_err) + call MPI_SEND(n, 1, MPI_INTEGER, 0, 0, mpi_intracomm, mpi_err) + call MPI_SEND(data, 2*n, MPI_INTEGER, 0, 1, mpi_intracomm, mpi_err) deallocate(data) #endif end if From c9e6bbd6c68d6d42ffd6b4917e69d11ec36aca9e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Feb 2017 11:47:28 -0600 Subject: [PATCH 11/29] Add comment describing intracomm argument --- src/initialize.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/initialize.F90 b/src/initialize.F90 index 2101333b7..d30f1e014 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -48,9 +48,9 @@ contains subroutine openmc_init(intracomm) #ifdef MPIF08 - type(MPI_Comm), intent(in) :: intracomm + type(MPI_Comm), intent(in) :: intracomm ! MPI intracommunicator #else - integer, intent(in), optional :: intracomm + integer, intent(in), optional :: intracomm ! MPI intracommunicator #endif ! Start total and initialization timer @@ -167,9 +167,9 @@ contains subroutine initialize_mpi(intracomm) #ifdef MPIF08 - type(MPI_Comm), intent(in) :: intracomm + type(MPI_Comm), intent(in) :: intracomm ! MPI intracommunicator #else - integer, intent(in) :: intracomm + integer, intent(in) :: intracomm ! MPI intracommunicator #endif integer :: bank_blocks(5) ! Count for each datatype From 5d2ad92769bdec7615fa0ff0bbc94e36d00953c3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Feb 2017 10:37:35 -0600 Subject: [PATCH 12/29] Add depletable attribute to Material --- openmc/material.py | 119 +++++++++++++++++++++++++-------------------- 1 file changed, 66 insertions(+), 53 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index f276f0f1a..7c271cc30 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -56,6 +56,9 @@ class Material(object): Units used for `density`. Can be one of 'g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum', or 'macro'. The 'macro' unit only applies in the case of a multi-group calculation. + depletable : bool + Indicate whether the material is depletable. This attribute can be used + by downstream depletion applications. elements : list of tuple List in which each item is a 4-tuple consisting of an :class:`openmc.Element` instance, the percent density, the percent @@ -78,6 +81,7 @@ class Material(object): self.temperature = temperature self._density = None self._density_units = '' + self._depletable = False # A list of tuples (nuclide, percent, percent type) self._nuclides = [] @@ -127,37 +131,36 @@ class Material(object): def __repr__(self): string = 'Material\n' - string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) - string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) - string += '{0: <16}{1}{2}\n'.format('\tTemperature', '=\t', - self._temperature) + string += '{: <16}=\t{}\n'.format('\tID', self._id) + string += '{: <16}=\t{}\n'.format('\tName', self._name) + string += '{: <16}=\t{}\n'.format('\tTemperature', self._temperature) - string += '{0: <16}{1}{2}'.format('\tDensity', '=\t', self._density) - string += ' [{0}]\n'.format(self._density_units) + string += '{: <16}=\t{}'.format('\tDensity', self._density) + string += ' [{}]\n'.format(self._density_units) - string += '{0: <16}\n'.format('\tS(a,b) Tables') + string += '{: <16}\n'.format('\tS(a,b) Tables') for sab in self._sab: - string += '{0: <16}{1}{2}\n'.format('\tS(a,b)', '=\t', sab) + string += '{: <16}=\t{}\n'.format('\tS(a,b)', sab) - string += '{0: <16}\n'.format('\tNuclides') + string += '{: <16}\n'.format('\tNuclides') for nuclide, percent, percent_type in self._nuclides: string += '{0: <16}'.format('\t{0.name}'.format(nuclide)) - string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type) + string += '=\t{: <12} [{}]\n'.format(percent, percent_type) if self._macroscopic is not None: - string += '{0: <16}\n'.format('\tMacroscopic Data') - string += '{0: <16}'.format('\t{0}'.format(self._macroscopic)) + string += '{: <16}\n'.format('\tMacroscopic Data') + string += '{: <16}'.format('\t{}'.format(self._macroscopic)) - string += '{0: <16}\n'.format('\tElements') + string += '{: <16}\n'.format('\tElements') for element, percent, percent_type, enr in self._elements: string += '{0: <16}'.format('\t{0.name}'.format(element)) if enr is None: - string += '=\t{0: <12} [{1}]\n'.format(percent, percent_type) + string += '=\t{: <12} [{}]\n'.format(percent, percent_type) else: - string += '=\t{0: <12} [{1}] @ {2} w/o enrichment\n'\ + string += '=\t{: <12} [{}] @ {} w/o enrichment\n'\ .format(percent, percent_type, enr) return string @@ -182,6 +185,10 @@ class Material(object): def density_units(self): return self._density_units + @property + def depletable(self): + return self._depletable + @property def elements(self): return self._elements @@ -234,7 +241,7 @@ class Material(object): @name.setter def name(self, name): if name is not None: - cv.check_type('name for Material ID="{0}"'.format(self._id), + cv.check_type('name for Material ID="{}"'.format(self._id), name, string_types) self._name = name else: @@ -242,10 +249,16 @@ class Material(object): @temperature.setter def temperature(self, temperature): - cv.check_type('Temperature for Material ID="{0}"'.format(self._id), + cv.check_type('Temperature for Material ID="{}"'.format(self._id), temperature, (Real, type(None))) self._temperature = temperature + @depletable.setter + def depletable(self, depletable): + cv.check_type('Depletable flag for Material ID="{}"'.format(self._id), + depletable, bool) + self._depletable = depletable + def set_density(self, units, density=None): """Set the density of the material @@ -264,17 +277,17 @@ class Material(object): if units is 'sum': if density is not None: - msg = 'Density "{0}" for Material ID="{1}" is ignored ' \ + msg = 'Density "{}" for Material ID="{}" is ignored ' \ 'because the unit is "sum"'.format(density, self.id) warnings.warn(msg) else: if density is None: - msg = 'Unable to set the density for Material ID="{0}" ' \ + msg = 'Unable to set the density for Material ID="{}" ' \ 'because a density value must be given when not using ' \ '"sum" unit'.format(self.id) raise ValueError(msg) - cv.check_type('the density for Material ID="{0}"'.format(self.id), + cv.check_type('the density for Material ID="{}"'.format(self.id), density, Real) self._density = density @@ -285,8 +298,8 @@ class Material(object): 'version of openmc') if not isinstance(filename, string_types) and filename is not None: - msg = 'Unable to add OTF material file to Material ID="{0}" with a ' \ - 'non-string name "{1}"'.format(self._id, filename) + msg = 'Unable to add OTF material file to Material ID="{}" with a ' \ + 'non-string name "{}"'.format(self._id, filename) raise ValueError(msg) self._distrib_otf_file = filename @@ -314,23 +327,23 @@ class Material(object): """ if self._macroscopic is not None: - msg = 'Unable to add a Nuclide to Material ID="{0}" as a ' \ + msg = 'Unable to add a Nuclide to Material ID="{}" as a ' \ 'macroscopic data-set has already been added'.format(self._id) raise ValueError(msg) if not isinstance(nuclide, string_types + (openmc.Nuclide,)): - msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \ - 'non-Nuclide value "{1}"'.format(self._id, nuclide) + msg = 'Unable to add a Nuclide to Material ID="{}" with a ' \ + 'non-Nuclide value "{}"'.format(self._id, nuclide) raise ValueError(msg) elif not isinstance(percent, Real): - msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \ - 'non-floating point value "{1}"'.format(self._id, percent) + msg = 'Unable to add a Nuclide to Material ID="{}" with a ' \ + 'non-floating point value "{}"'.format(self._id, percent) raise ValueError(msg) elif percent_type not in ['ao', 'wo', 'at/g-cm']: - msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \ - 'percent type "{1}"'.format(self._id, percent_type) + msg = 'Unable to add a Nuclide to Material ID="{}" with a ' \ + 'percent type "{}"'.format(self._id, percent_type) raise ValueError(msg) if isinstance(nuclide, openmc.Nuclide): @@ -353,7 +366,7 @@ class Material(object): """ if not isinstance(nuclide, openmc.Nuclide): - msg = 'Unable to remove a Nuclide "{0}" in Material ID="{1}" ' \ + msg = 'Unable to remove a Nuclide "{}" in Material ID="{}" ' \ 'since it is not a Nuclide'.format(self._id, nuclide) raise ValueError(msg) @@ -377,15 +390,15 @@ class Material(object): # Ensure no nuclides, elements, or sab are added since these would be # incompatible with macroscopics if self._nuclides or self._elements or self._sab: - msg = 'Unable to add a Macroscopic data set to Material ID="{0}" ' \ - 'with a macroscopic value "{1}" as an incompatible data ' \ + msg = 'Unable to add a Macroscopic data set to Material ID="{}" ' \ + 'with a macroscopic value "{}" as an incompatible data ' \ 'member (i.e., nuclide, element, or S(a,b) table) ' \ 'has already been added'.format(self._id, macroscopic) raise ValueError(msg) if not isinstance(macroscopic, string_types + (openmc.Macroscopic,)): - msg = 'Unable to add a Macroscopic to Material ID="{0}" with a ' \ - 'non-Macroscopic value "{1}"'.format(self._id, macroscopic) + msg = 'Unable to add a Macroscopic to Material ID="{}" with a ' \ + 'non-Macroscopic value "{}"'.format(self._id, macroscopic) raise ValueError(msg) if isinstance(macroscopic, openmc.Macroscopic): @@ -398,7 +411,7 @@ class Material(object): if self._macroscopic is None: self._macroscopic = macroscopic else: - msg = 'Unable to add a Macroscopic to Material ID="{0}". ' \ + msg = 'Unable to add a Macroscopic to Material ID="{}". ' \ 'Only one Macroscopic allowed per ' \ 'Material.'.format(self._id) raise ValueError(msg) @@ -422,7 +435,7 @@ class Material(object): """ if not isinstance(macroscopic, openmc.Macroscopic): - msg = 'Unable to remove a Macroscopic "{0}" in Material ID="{1}" ' \ + msg = 'Unable to remove a Macroscopic "{}" in Material ID="{}" ' \ 'since it is not a Macroscopic'.format(self._id, macroscopic) raise ValueError(msg) @@ -450,23 +463,23 @@ class Material(object): """ if self._macroscopic is not None: - msg = 'Unable to add an Element to Material ID="{0}" as a ' \ + msg = 'Unable to add an Element to Material ID="{}" as a ' \ 'macroscopic data-set has already been added'.format(self._id) raise ValueError(msg) if not isinstance(element, string_types + (openmc.Element,)): - msg = 'Unable to add an Element to Material ID="{0}" with a ' \ - 'non-Element value "{1}"'.format(self._id, element) + msg = 'Unable to add an Element to Material ID="{}" with a ' \ + 'non-Element value "{}"'.format(self._id, element) raise ValueError(msg) if not isinstance(percent, Real): - msg = 'Unable to add an Element to Material ID="{0}" with a ' \ - 'non-floating point value "{1}"'.format(self._id, percent) + msg = 'Unable to add an Element to Material ID="{}" with a ' \ + 'non-floating point value "{}"'.format(self._id, percent) raise ValueError(msg) if percent_type not in ['ao', 'wo']: - msg = 'Unable to add an Element to Material ID="{0}" with a ' \ - 'percent type "{1}"'.format(self._id, percent_type) + msg = 'Unable to add an Element to Material ID="{}" with a ' \ + 'percent type "{}"'.format(self._id, percent_type) raise ValueError(msg) # Copy this Element to separate it from same Element in other Materials @@ -477,14 +490,14 @@ class Material(object): if enrichment is not None: if not isinstance(enrichment, Real): - msg = 'Unable to add an Element to Material ID="{0}" with a ' \ - 'non-floating point enrichment value "{1}"'\ + msg = 'Unable to add an Element to Material ID="{}" with a ' \ + 'non-floating point enrichment value "{}"'\ .format(self._id, enrichment) raise ValueError(msg) elif element.name != 'U': - msg = 'Unable to use enrichment for element {0} which is not ' \ - 'uranium for Material ID="{1}"'.format(element.name, + msg = 'Unable to use enrichment for element {} which is not ' \ + 'uranium for Material ID="{}"'.format(element.name, self._id) raise ValueError(msg) @@ -493,8 +506,8 @@ class Material(object): cv.check_greater_than('enrichment', enrichment, 0., equality=True) if enrichment > 5.0: - msg = 'A uranium enrichment of {0} was given for Material ID='\ - '"{1}". OpenMC assumes the U234/U235 mass ratio is '\ + msg = 'A uranium enrichment of {} was given for Material ID='\ + '"{}". OpenMC assumes the U234/U235 mass ratio is '\ 'constant at 0.008, which is only valid at low ' \ 'enrichments. Consider setting the isotopic ' \ 'composition manually for enrichments over 5%.'.\ @@ -514,7 +527,7 @@ class Material(object): """ if not isinstance(element, openmc.Element): - msg = 'Unable to remove "{0}" in Material ID="{1}" ' \ + msg = 'Unable to remove "{}" in Material ID="{}" ' \ 'since it is not an Element'.format(self.id, element) raise ValueError(msg) @@ -534,13 +547,13 @@ class Material(object): """ if self._macroscopic is not None: - msg = 'Unable to add an S(a,b) table to Material ID="{0}" as a ' \ + msg = 'Unable to add an S(a,b) table to Material ID="{}" as a ' \ 'macroscopic data-set has already been added'.format(self._id) raise ValueError(msg) if not isinstance(name, string_types): - msg = 'Unable to add an S(a,b) table to Material ID="{0}" with a ' \ - 'non-string table name "{1}"'.format(self._id, name) + msg = 'Unable to add an S(a,b) table to Material ID="{}" with a ' \ + 'non-string table name "{}"'.format(self._id, name) raise ValueError(msg) new_name = openmc.data.get_thermal_name(name) From 384a08956af7a92ae3fc70292b27aac698f091ca Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Thu, 16 Feb 2017 21:02:34 -0500 Subject: [PATCH 13/29] Added test of convert_* methods; --- openmc/mgxs_library.py | 92 +++++----- tests/test_mg_convert/inputs_true.dat | 30 ++++ tests/test_mg_convert/results_true.dat | 28 +++ tests/test_mg_convert/test_mg_convert.py | 207 +++++++++++++++++++++++ 4 files changed, 303 insertions(+), 54 deletions(-) create mode 100644 tests/test_mg_convert/inputs_true.dat create mode 100644 tests/test_mg_convert/results_true.dat create mode 100755 tests/test_mg_convert/test_mg_convert.py diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index e3a8ad4d8..d441e1bbb 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -14,11 +14,14 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \ # Supported incoming particle MGXS angular treatment representations _REPRESENTATIONS = ['isotropic', 'angle'] + # Supported scattering angular distribution representations _SCATTER_TYPES = ['tabular', 'legendre', 'histogram'] -# List of MGXS dimension types + +# List of MGXS indexing schemes _XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[DG][G]", "[DG][G']", "[DG][G][G']"] + # Number of mu points for conversion between scattering formats _NMU = 257 @@ -368,15 +371,14 @@ class XSdata(object): @name.setter def name(self, name): + check_type('name for XSdata', name, string_types) self._name = name @energy_groups.setter def energy_groups(self, energy_groups): - # Check validity of energy_groups check_type('energy_groups', energy_groups, openmc.mgxs.EnergyGroups) - if energy_groups.group_edges is None: msg = 'Unable to assign an EnergyGroups object ' \ 'with uninitialized group edges' @@ -387,7 +389,6 @@ class XSdata(object): @num_delayed_groups.setter def num_delayed_groups(self, num_delayed_groups): - # Check validity of num_delayed_groups check_type('num_delayed_groups', num_delayed_groups, Integral) check_less_than('num_delayed_groups', num_delayed_groups, openmc.mgxs.MAX_DELAYED_GROUPS, equality=True) @@ -398,14 +399,12 @@ class XSdata(object): @representation.setter def representation(self, representation): - # Check it is of valid value. check_value('representation', representation, _REPRESENTATIONS) self._representation = representation @atomic_weight_ratio.setter def atomic_weight_ratio(self, atomic_weight_ratio): - # Check validity of type and that the atomic_weight_ratio value is > 0 check_type('atomic_weight_ratio', atomic_weight_ratio, Real) check_greater_than('atomic_weight_ratio', atomic_weight_ratio, 0.0) self._atomic_weight_ratio = atomic_weight_ratio @@ -419,14 +418,12 @@ class XSdata(object): @scatter_format.setter def scatter_format(self, scatter_format): - # check to see it is of a valid type and value check_value('scatter_format', scatter_format, _SCATTER_TYPES) self._scatter_format = scatter_format @order.setter def order(self, order): - # Check type and value check_type('order', order, Integral) check_greater_than('order', order, 0, equality=True) self._order = order @@ -434,7 +431,6 @@ class XSdata(object): @num_polar.setter def num_polar(self, num_polar): - # Make sure we have positive ints check_type('num_polar', num_polar, Integral) check_greater_than('num_polar', num_polar, 0) self._num_polar = num_polar @@ -1688,7 +1684,14 @@ class XSdata(object): def convert_representation(self, target_representation, num_polar=None, num_azimuthal=None): """Produce a new XSdata object with the same data, but converted to the - new representation + new representation (isotropic or angle-dependent). + + This method cannot be used to change the number of polar or + azimuthal bins of an XSdata object that already uses an angular + representation. Finally, this method simply uses an arithmetic mean to + convert from an angular to isotropic representation; no flux-weighting + is applied and therefore the correctness of the solution is not + guaranteed. Parameters ---------- @@ -1728,9 +1731,9 @@ class XSdata(object): # Check to make sure the num_polar and num_azimuthal values match if target_representation == 'angle': if num_polar != self.num_polar or num_azimuthal != self.num_azimuthal: - raise NotImplementedError("XCannot translate between " - "`angle` representations with " - "different angle bin structures") + raise ValueError("Cannot translate between `angle`" + " representations with different angle" + " bin structures") # Nothing to do as the same structure was requested return xsdata @@ -1741,6 +1744,7 @@ class XSdata(object): # values are changed back to None for clarity xsdata._num_polar = None xsdata._num_azimuthal = None + elif target_representation == 'angle': xsdata.num_polar = num_polar xsdata.num_azimuthal = num_azimuthal @@ -1757,16 +1761,19 @@ class XSdata(object): # Get the original data orig_data = getattr(self, '_' + xs)[i] if orig_data is not None: + if target_representation == 'isotropic': # Since we are going from angle to isotropic, the # current data is just the average over the angle bins new_data = orig_data.mean(axis=(0, 1)) + elif target_representation == 'angle': # Since we are going from isotropic to angle, the # current data is just copied for every angle bin new_shape = (num_polar, num_azimuthal) + \ orig_data.shape new_data = np.resize(orig_data, new_shape) + setter = getattr(xsdata, 'set_' + xs) setter(new_data, temp) @@ -1810,18 +1817,19 @@ class XSdata(object): # Reset and re-generate XSdata.xs_shapes with the new scattering format xsdata._xs_shapes = None - xsdata.xs_shapes for i, temp in enumerate(xsdata.temperatures): orig_data = self._scatter_matrix[i] new_shape = orig_data.shape[:-1] + (xsdata.num_orders,) new_data = np.zeros(new_shape) + if self.scatter_format == 'legendre': if target_format == 'legendre': # Then we are changing orders and only need to change # dimensionality of the mu data and pad/truncate as needed order = min(xsdata.num_orders, self.num_orders) new_data[..., :order] = orig_data[..., :order] + elif target_format == 'tabular': mu = np.linspace(-1, 1, xsdata.num_orders) # Evaluate the legendre on the mu grid @@ -1856,6 +1864,7 @@ class XSdata(object): elif self.scatter_format == 'tabular': # Calculate the mu points of the current data mu_self = np.linspace(-1, 1, self.num_orders) + if target_format == 'legendre': # Find the Legendre coefficients via integration. To best # use the vectorized integration capabilities of scipy, @@ -1870,10 +1879,12 @@ class XSdata(object): # Remove the very small values resulting from numerical # precision issues new_data[..., np.abs(new_data) < 1.E-10] = 0. + elif target_format == 'tabular': # Simply use an interpolating function to get the new data mu = np.linspace(-1, 1, xsdata.num_orders) new_data[..., :] = interp1d(mu_self, orig_data)(mu) + elif target_format == 'histogram': # Use an interpolating function to do the bin-wise # integrals @@ -1891,6 +1902,7 @@ class XSdata(object): # Remove the very small values resulting from numerical # precision issues new_data[..., np.abs(new_data) < 1.E-10] = 0. + elif self.scatter_format == 'histogram': # The histogram format does not have enough information to # convert to the other forms without inducing some amount of @@ -1921,10 +1933,12 @@ class XSdata(object): # Remove the very small values resulting from numerical # precision issues new_data[..., np.abs(new_data) < 1.E-10] = 0. + elif target_format == 'tabular': # Simply use an interpolating function to get the new data mu = np.linspace(-1, 1, xsdata.num_orders) new_data[..., :] = interp(mu) * norm + elif target_format == 'histogram': # Use an interpolating function to do the bin-wise # integrals @@ -1942,6 +1956,7 @@ class XSdata(object): # Remove the very small values resulting from numerical # precision issues new_data[..., np.abs(new_data) < 1.E-10] = 0. + xsdata.set_scatter_matrix(new_data, temp) return xsdata @@ -2431,8 +2446,15 @@ class MGXSLibrary(object): def convert_representation(self, target_representation, num_polar=None, num_azimuthal=None): - """Produce a new MGXSLibrary object with the same data, but converted - to the new representation + """Produce a new XSdata object with the same data, but converted to the + new representation (isotropic or angle-dependent). + + This method cannot be used to change the number of polar or + azimuthal bins of an XSdata object that already uses an angular + representation. Finally, this method simply uses an arithmetic mean to + convert from an angular to isotropic representation; no flux-weighting + is applied and therefore the correctness of the solution is not + guaranteed. Parameters ---------- @@ -2456,14 +2478,6 @@ class MGXSLibrary(object): """ - check_value('target_representation', target_representation, - _REPRESENTATIONS) - if target_representation == 'angle': - check_type('num_polar', num_polar, Integral) - check_type('num_azimuthal', num_azimuthal, Integral) - check_greater_than('num_polar', num_polar, 0) - check_greater_than('num_azimuthal', num_azimuthal, 0) - library = copy.deepcopy(self) for i, xsdata in enumerate(self.xsdatas): library.xsdatas[i] = \ @@ -2493,13 +2507,6 @@ class MGXSLibrary(object): """ - check_value('target_format', target_format, _SCATTER_TYPES) - check_type('target_order', target_order, Integral) - if target_format == 'legendre': - check_greater_than('target_order', target_order, 0, equality=True) - else: - check_greater_than('target_order', target_order, 0) - library = copy.deepcopy(self) for i, xsdata in enumerate(self.xsdatas): library.xsdatas[i] = \ @@ -2507,29 +2514,6 @@ class MGXSLibrary(object): return library - def convert_to_continuous_energy(self, h5_filename='ce_mgxs.h5', - library_filename='ce_mgxs.xml'): - """Converts the MGXSLibrary object to an equivalent - library of openmc.data.IncidentNeutron objects - - Parameters - ---------- - h5_filename : str - HDF5 file to write with all the files - library_filename : str - cross_sections.xml file describing the HDF5 file - - """ - - library = openmc.data.DataLibrary() - data = [] - for i, xsdata in enumerate(self.xsdatas): - data.append(xsdata.convert_to_continuous_energy()) - data[-1].export_to_hdf5(h5_filename) - - library.register_file(h5_filename) - library.export_to_xml(library_filename) - def export_to_hdf5(self, filename='mgxs.h5'): """Create an hdf5 file that can be used for a simulation. diff --git a/tests/test_mg_convert/inputs_true.dat b/tests/test_mg_convert/inputs_true.dat new file mode 100644 index 000000000..9b0b0f4b1 --- /dev/null +++ b/tests/test_mg_convert/inputs_true.dat @@ -0,0 +1,30 @@ + + + + + + + + + + + ./mgxs.h5 + + + + + + + + + 1000 + 2000 + 100 + + + + -5 -5 -5 5 5 5 + + + multi-group + diff --git a/tests/test_mg_convert/results_true.dat b/tests/test_mg_convert/results_true.dat new file mode 100644 index 000000000..5cf0eb6d2 --- /dev/null +++ b/tests/test_mg_convert/results_true.dat @@ -0,0 +1,28 @@ +k-combined: +9.957263E-01 5.469041E-05 +k-combined: +9.957263E-01 5.469041E-05 +k-combined: +9.963652E-01 5.444550E-05 +k-combined: +9.957263E-01 5.469041E-05 +k-combined: +9.955157E-01 5.785520E-05 +k-combined: +9.954023E-01 5.694553E-05 +k-combined: +9.955722E-01 6.075634E-05 +k-combined: +9.955692E-01 6.094480E-05 +k-combined: +9.955154E-01 5.775805E-05 +k-combined: +9.956787E-01 5.915112E-05 +k-combined: +9.957022E-01 5.952369E-05 +k-combined: +9.953517E-01 5.920431E-05 +k-combined: +9.957263E-01 5.469041E-05 +k-combined: +9.957263E-01 5.469041E-05 diff --git a/tests/test_mg_convert/test_mg_convert.py b/tests/test_mg_convert/test_mg_convert.py new file mode 100755 index 000000000..c5c0ac16c --- /dev/null +++ b/tests/test_mg_convert/test_mg_convert.py @@ -0,0 +1,207 @@ +#!/usr/bin/env python + +import os +import sys +import hashlib +sys.path.insert(0, os.pardir) + +import numpy as np + +from testing_harness import PyAPITestHarness +import openmc + +# OpenMC simulation parameters +batches = 2000 +inactive = 100 +particles = 1000 + + +def build_mgxs_library(convert): + # Instantiate the energy group data + groups = openmc.mgxs.EnergyGroups(group_edges=[1e-5, 0.625, 20.0e6]) + + # Instantiate the 7-group (C5G7) cross section data + uo2_xsdata = openmc.XSdata('UO2', groups) + uo2_xsdata.order = 2 + uo2_xsdata.set_total([2., 2.]) + uo2_xsdata.set_absorption([1., 1.]) + scatter_matrix = np.array([[[0.75, 0.25], + [0.00, 1.00]], + [[0.75 / 3., 0.25 / 3.], + [0.00 / 3., 1.00 / 3.]], + [[0.75 / 4., 0.25 / 4.], + [0.00 / 4., 1.00 / 4.]]]) + scatter_matrix = np.rollaxis(scatter_matrix, 0, 3) + uo2_xsdata.set_scatter_matrix(scatter_matrix) + uo2_xsdata.set_fission([0.5, 0.5]) + uo2_xsdata.set_nu_fission([1., 1.]) + uo2_xsdata.set_chi([1., 0.]) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.add_xsdatas([uo2_xsdata]) + + if convert is not None: + if isinstance(convert[0], list): + for conv in convert: + if conv[0] in ['legendre', 'tabular', 'histogram']: + mg_cross_sections_file = \ + mg_cross_sections_file.convert_scatter_format( + conv[0], conv[1]) + elif conv[0] in ['angle', 'isotropic']: + mg_cross_sections_file = \ + mg_cross_sections_file.convert_representation( + conv[0], conv[1], conv[1]) + elif convert[0] in ['legendre', 'tabular', 'histogram']: + mg_cross_sections_file = \ + mg_cross_sections_file.convert_scatter_format( + convert[0], convert[1]) + elif convert[0] in ['angle', 'isotropic']: + mg_cross_sections_file = \ + mg_cross_sections_file.convert_representation( + convert[0], convert[1], convert[1]) + + mg_cross_sections_file.export_to_hdf5() + + +class MGXSTestHarness(PyAPITestHarness): + def _build_inputs(self): + # Instantiate some Macroscopic Data + uo2_data = openmc.Macroscopic('UO2') + + # Instantiate some Materials and register the appropriate objects + mat = openmc.Material(material_id=1, name='UO2 fuel') + mat.set_density('macro', 1.0) + mat.add_macroscopic(uo2_data) + + # Instantiate a Materials collection and export to XML + materials_file = openmc.Materials([mat]) + materials_file.cross_sections = "./mgxs.h5" + materials_file.export_to_xml() + + # Instantiate ZCylinder surfaces + left = openmc.XPlane(surface_id=4, x0=-5., name='left') + right = openmc.XPlane(surface_id=5, x0=5., name='right') + bottom = openmc.YPlane(surface_id=6, y0=-5., name='bottom') + top = openmc.YPlane(surface_id=7, y0=5., name='top') + + left.boundary_type = 'reflective' + right.boundary_type = 'vacuum' + top.boundary_type = 'reflective' + bottom.boundary_type = 'reflective' + + # Instantiate Cells + fuel = openmc.Cell(cell_id=1, name='cell 1') + + # Use surface half-spaces to define regions + fuel.region = +left & -right & +bottom & -top + + # Register Materials with Cells + fuel.fill = mat + + # Instantiate Universe + root = openmc.Universe(universe_id=0, name='root universe') + + # Register Cells with Universe + root.add_cells([fuel]) + + # Instantiate a Geometry, register the root Universe, and export to XML + geometry = openmc.Geometry(root) + geometry.export_to_xml() + + settings_file = openmc.Settings() + settings_file.energy_mode = "multi-group" + settings_file.batches = batches + settings_file.inactive = inactive + settings_file.particles = particles + + # Create an initial uniform spatial source distribution + bounds = [-5, -5, -5, 5, 5, 5] + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) + settings_file.source = openmc.source.Source(space=uniform_dist) + + settings_file.export_to_xml() + + def _run_openmc(self): + # Run multiple conversions to compare results + cases = [None, ['legendre', 2], ['legendre', 0], + ['tabular', 33], ['histogram', 32], + [['tabular', 33], ['legendre', 1]], + [['tabular', 33], ['tabular', 3]], + [['tabular', 33], ['histogram', 32]], + [['histogram', 32], ['legendre', 1]], + [['histogram', 32], ['tabular', 3]], + [['histogram', 32], ['histogram', 16]], + [['histogram', 32], ['histogram', 128]], + ['angle', 2], [['angle', 2], ['isotropic', None]]] + + outstr = '' + for case in cases: + build_mgxs_library(case) + + if self._opts.mpi_exec is not None: + returncode = openmc.run(mpi_procs=self._opts.mpi_np, + openmc_exec=self._opts.exe, + mpi_exec=self._opts.mpi_exec) + + else: + returncode = openmc.run(openmc_exec=self._opts.exe) + + assert returncode == 0, 'OpenMC did not exit successfully.' + + sp = openmc.StatePoint('statepoint.' + str(batches) + '.h5') + + # Write out k-combined. + outstr += 'k-combined:\n' + form = '{0:12.6E} {1:12.6E}\n' + outstr += form.format(sp.k_combined[0], sp.k_combined[1]) + sp.close() + + return outstr + + def _get_results(self, outstr, hash_output=False): + # Hash the results if necessary. + if hash_output: + sha512 = hashlib.sha512() + sha512.update(outstr.encode('utf-8')) + outstr = sha512.hexdigest() + + return outstr + + def _cleanup(self): + super(MGXSTestHarness, self)._cleanup() + f = os.path.join(os.getcwd(), 'mgxs.h5') + if os.path.exists(f): + os.remove(f) + + def execute_test(self): + """Build input XMLs, run OpenMC, and verify correct results.""" + try: + self._build_inputs() + inputs = self._get_inputs() + self._write_inputs(inputs) + self._compare_inputs() + outstr = self._run_openmc() + results = self._get_results(outstr) + self._write_results(results) + self._compare_results() + finally: + self._cleanup() + + def update_results(self): + """Update results_true.dat and inputs_true.dat""" + try: + self._build_inputs() + inputs = self._get_inputs() + self._write_inputs(inputs) + self._overwrite_inputs() + outstr = self._run_openmc() + results = self._get_results(outstr) + self._write_results(results) + self._overwrite_results() + finally: + self._cleanup() + + +if __name__ == '__main__': + harness = MGXSTestHarness('statepoint.10.*', False) + harness.main() From 7069c5db315be9f2a756d6e5d8cc295d7be474a0 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 17 Feb 2017 07:48:45 -0500 Subject: [PATCH 14/29] lessened the number of particles used in test_mg_convert and reduced the number of cases by 2 --- tests/test_mg_convert/inputs_true.dat | 2 +- tests/test_mg_convert/results_true.dat | 28 ++++++++++-------------- tests/test_mg_convert/test_mg_convert.py | 5 ++--- 3 files changed, 15 insertions(+), 20 deletions(-) diff --git a/tests/test_mg_convert/inputs_true.dat b/tests/test_mg_convert/inputs_true.dat index 9b0b0f4b1..95ad61586 100644 --- a/tests/test_mg_convert/inputs_true.dat +++ b/tests/test_mg_convert/inputs_true.dat @@ -18,7 +18,7 @@ 1000 - 2000 + 1000 100 diff --git a/tests/test_mg_convert/results_true.dat b/tests/test_mg_convert/results_true.dat index 5cf0eb6d2..cf4777b8e 100644 --- a/tests/test_mg_convert/results_true.dat +++ b/tests/test_mg_convert/results_true.dat @@ -1,28 +1,24 @@ k-combined: -9.957263E-01 5.469041E-05 +9.957366E-01 7.836427E-05 k-combined: -9.957263E-01 5.469041E-05 +9.963227E-01 8.010344E-05 k-combined: -9.963652E-01 5.444550E-05 +9.957366E-01 7.836427E-05 k-combined: -9.957263E-01 5.469041E-05 +9.955152E-01 8.132242E-05 k-combined: -9.955157E-01 5.785520E-05 +9.953937E-01 8.452497E-05 k-combined: -9.954023E-01 5.694553E-05 +9.957744E-01 8.816057E-05 k-combined: -9.955722E-01 6.075634E-05 +9.955332E-01 8.894390E-05 k-combined: -9.955692E-01 6.094480E-05 +9.956463E-01 8.378533E-05 k-combined: -9.955154E-01 5.775805E-05 +9.955875E-01 8.818927E-05 k-combined: -9.956787E-01 5.915112E-05 +9.957017E-01 8.726677E-05 k-combined: -9.957022E-01 5.952369E-05 +9.957366E-01 7.836427E-05 k-combined: -9.953517E-01 5.920431E-05 -k-combined: -9.957263E-01 5.469041E-05 -k-combined: -9.957263E-01 5.469041E-05 +9.957366E-01 7.836427E-05 diff --git a/tests/test_mg_convert/test_mg_convert.py b/tests/test_mg_convert/test_mg_convert.py index c5c0ac16c..ad86e5a50 100755 --- a/tests/test_mg_convert/test_mg_convert.py +++ b/tests/test_mg_convert/test_mg_convert.py @@ -11,7 +11,7 @@ from testing_harness import PyAPITestHarness import openmc # OpenMC simulation parameters -batches = 2000 +batches = 1000 inactive = 100 particles = 1000 @@ -123,7 +123,7 @@ class MGXSTestHarness(PyAPITestHarness): def _run_openmc(self): # Run multiple conversions to compare results - cases = [None, ['legendre', 2], ['legendre', 0], + cases = [['legendre', 2], ['legendre', 0], ['tabular', 33], ['histogram', 32], [['tabular', 33], ['legendre', 1]], [['tabular', 33], ['tabular', 3]], @@ -131,7 +131,6 @@ class MGXSTestHarness(PyAPITestHarness): [['histogram', 32], ['legendre', 1]], [['histogram', 32], ['tabular', 3]], [['histogram', 32], ['histogram', 16]], - [['histogram', 32], ['histogram', 128]], ['angle', 2], [['angle', 2], ['isotropic', None]]] outstr = '' From ffc3bbfd67f0584820d3100801986b7f10d34ab8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Feb 2017 17:10:19 -0600 Subject: [PATCH 15/29] Allow plot, volume, and particle restart run modes directly --- examples/xml/basic/settings.xml | 10 +- examples/xml/boxes/settings.xml | 9 +- examples/xml/lattice/nested/settings.xml | 9 +- examples/xml/lattice/simple/settings.xml | 9 +- examples/xml/pincell/settings.xml | 11 +- examples/xml/pincell_multigroup/settings.xml | 18 +- examples/xml/reflective/settings.xml | 9 +- openmc/settings.py | 32 ++- src/constants.F90 | 3 +- src/input_xml.F90 | 211 ++++++++++-------- tests/test_asymmetric_lattice/inputs_true.dat | 9 +- tests/test_cmfd_feed/settings.xml | 11 +- tests/test_cmfd_nofeed/settings.xml | 11 +- tests/test_complex_cell/settings.xml | 9 +- tests/test_confidence_intervals/settings.xml | 9 +- .../inputs_true.dat | 7 +- tests/test_density/settings.xml | 9 +- tests/test_diff_tally/inputs_true.dat | 9 +- tests/test_distribmat/inputs_true.dat | 9 +- .../test_eigenvalue_genperbatch/settings.xml | 11 +- .../test_eigenvalue_no_inactive/settings.xml | 9 +- tests/test_energy_cutoff/inputs_true.dat | 7 +- tests/test_energy_grid/settings.xml | 9 +- tests/test_energy_laws/settings.xml | 9 +- tests/test_entropy/settings.xml | 9 +- .../case-1/settings.xml | 9 +- .../case-2/settings.xml | 9 +- .../case-3/settings.xml | 9 +- .../case-4/settings.xml | 19 +- tests/test_filter_energyfun/inputs_true.dat | 9 +- tests/test_filter_mesh/inputs_true.dat | 9 +- tests/test_fixed_source/settings.xml | 7 +- tests/test_infinite_cell/settings.xml | 9 +- tests/test_iso_in_lab/inputs_true.dat | 9 +- tests/test_lattice/settings.xml | 9 +- tests/test_lattice_hex/settings.xml | 9 +- tests/test_lattice_mixed/settings.xml | 9 +- tests/test_lattice_multiple/settings.xml | 9 +- tests/test_mg_basic/inputs_true.dat | 9 +- tests/test_mg_max_order/inputs_true.dat | 9 +- tests/test_mg_nuclide/inputs_true.dat | 9 +- tests/test_mg_tallies/inputs_true.dat | 9 +- .../inputs_true.dat | 9 +- .../inputs_true.dat | 9 +- .../inputs_true.dat | 9 +- tests/test_mgxs_library_hdf5/inputs_true.dat | 9 +- tests/test_mgxs_library_mesh/inputs_true.dat | 9 +- .../inputs_true.dat | 9 +- .../inputs_true.dat | 9 +- tests/test_multipole/inputs_true.dat | 9 +- tests/test_output/settings.xml | 9 +- .../test_particle_restart_eigval/settings.xml | 9 +- .../test_particle_restart_fixed/settings.xml | 7 +- tests/test_periodic/inputs_true.dat | 9 +- tests/test_plot/settings.xml | 9 +- tests/test_ptables_off/settings.xml | 9 +- tests/test_quadric_surfaces/settings.xml | 9 +- tests/test_reflective_plane/settings.xml | 9 +- .../test_resonance_scattering/inputs_true.dat | 9 +- tests/test_rotation/settings.xml | 9 +- tests/test_salphabeta/settings.xml | 9 +- tests/test_score_current/settings.xml | 9 +- tests/test_seed/settings.xml | 9 +- tests/test_source/inputs_true.dat | 9 +- tests/test_sourcepoint_batch/settings.xml | 9 +- tests/test_sourcepoint_latest/settings.xml | 9 +- tests/test_sourcepoint_restart/settings.xml | 9 +- tests/test_statepoint_batch/settings.xml | 9 +- tests/test_statepoint_restart/settings.xml | 9 +- tests/test_statepoint_sourcesep/settings.xml | 9 +- tests/test_survival_biasing/settings.xml | 9 +- tests/test_tallies/inputs_true.dat | 9 +- tests/test_tally_aggregation/inputs_true.dat | 9 +- tests/test_tally_arithmetic/inputs_true.dat | 9 +- tests/test_tally_assumesep/settings.xml | 9 +- tests/test_tally_nuclides/settings.xml | 9 +- tests/test_tally_slice_merge/inputs_true.dat | 9 +- tests/test_trace/settings.xml | 9 +- tests/test_track_output/settings.xml | 18 +- tests/test_translation/settings.xml | 9 +- .../test_trigger_batch_interval/settings.xml | 17 +- .../settings.xml | 17 +- tests/test_trigger_no_status/settings.xml | 19 +- tests/test_trigger_tallies/settings.xml | 17 +- tests/test_triso/inputs_true.dat | 9 +- tests/test_uniform_fs/settings.xml | 9 +- tests/test_universe/settings.xml | 9 +- tests/test_void/settings.xml | 9 +- tests/test_volume_calc/inputs_true.dat | 9 +- 89 files changed, 499 insertions(+), 584 deletions(-) diff --git a/examples/xml/basic/settings.xml b/examples/xml/basic/settings.xml index eb9f5bd56..6e622b6ce 100644 --- a/examples/xml/basic/settings.xml +++ b/examples/xml/basic/settings.xml @@ -1,12 +1,10 @@ - - - 15 - 5 - 10000 - + eigenvalue + 15 + 5 + 10000 diff --git a/examples/xml/boxes/settings.xml b/examples/xml/boxes/settings.xml index eff7c1c10..9007a12c5 100644 --- a/examples/xml/boxes/settings.xml +++ b/examples/xml/boxes/settings.xml @@ -2,11 +2,10 @@ - - 15 - 5 - 10000 - + eigenvalue + 15 + 5 + 10000 diff --git a/examples/xml/lattice/nested/settings.xml b/examples/xml/lattice/nested/settings.xml index 2a6aaaf42..879173d1b 100644 --- a/examples/xml/lattice/nested/settings.xml +++ b/examples/xml/lattice/nested/settings.xml @@ -2,11 +2,10 @@ - - 20 - 10 - 10000 - + eigenvalue + 20 + 10 + 10000 diff --git a/examples/xml/lattice/simple/settings.xml b/examples/xml/lattice/simple/settings.xml index 2a6aaaf42..879173d1b 100644 --- a/examples/xml/lattice/simple/settings.xml +++ b/examples/xml/lattice/simple/settings.xml @@ -2,11 +2,10 @@ - - 20 - 10 - 10000 - + eigenvalue + 20 + 10 + 10000 diff --git a/examples/xml/pincell/settings.xml b/examples/xml/pincell/settings.xml index 443af9cde..733de1926 100644 --- a/examples/xml/pincell/settings.xml +++ b/examples/xml/pincell/settings.xml @@ -2,11 +2,10 @@ - - 100 - 10 - 1000 - + eigenvalue + 100 + 10 + 1000 - - 500 - 10 - 10000 - + eigenvalue + 500 + 10 + 10000 diff --git a/openmc/settings.py b/openmc/settings.py index e6a3f9e13..a4dde8da3 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -11,6 +11,9 @@ from openmc.clean_xml import clean_xml_indentation import openmc.checkvalue as cv from openmc import Nuclide, VolumeCalculation, Source, Mesh +_RUN_MODES = ['eigenvalue', 'fixed source', 'plot', 'volume', + 'particle restart'] + class Settings(object): """Settings used for an OpenMC simulation. @@ -81,7 +84,7 @@ class Settings(object): The elastic scattering model to use for resonant isotopes run_cmfd : bool Indicate if coarse mesh finite difference acceleration is to be used - run_mode : {'eigenvalue' or 'fixed source'} + run_mode : {'eigenvalue', 'fixed source', 'plot', 'volume', 'particle restart'} The type of calculation to perform (default is 'eigenvalue') seed : int Seed for the linear congruential pseudorandom number generator @@ -388,10 +391,7 @@ class Settings(object): @run_mode.setter def run_mode(self, run_mode): - if run_mode not in ['eigenvalue', 'fixed source']: - msg = 'Unable to set run mode to "{0}". Only "eigenvalue" ' \ - 'and "fixed source" are supported."'.format(run_mode) - raise ValueError(msg) + cv.check_value('run mode', run_mode, _RUN_MODES) self._run_mode = run_mode @batches.setter @@ -794,18 +794,8 @@ class Settings(object): self._create_fission_neutrons = create_fission_neutrons def _create_run_mode_subelement(self, root): - - if self.run_mode == 'eigenvalue': - elem = ET.SubElement(root, "eigenvalue") - self._create_particles_subelement(elem) - self._create_batches_subelement(elem) - self._create_inactive_subelement(elem) - self._create_generations_per_batch_subelement(elem) - self._create_keff_trigger_subelement(elem) - else: - elem = ET.SubElement(root, "fixed_source") - self._create_particles_subelement(elem) - self._create_batches_subelement(elem) + elem = ET.SubElement(root, "run_mode") + elem.text = self._run_mode def _create_batches_subelement(self, run_mode_element): if self._batches is not None: @@ -814,8 +804,7 @@ class Settings(object): def _create_generations_per_batch_subelement(self, run_mode_element): if self._generations_per_batch is not None: - element = ET.SubElement(run_mode_element, - "generations_per_batch") + element = ET.SubElement(run_mode_element, "generations_per_batch") element.text = str(self._generations_per_batch) def _create_inactive_subelement(self, run_mode_element): @@ -1081,6 +1070,11 @@ class Settings(object): root_element = ET.Element("settings") self._create_run_mode_subelement(root_element) + self._create_particles_subelement(root_element) + self._create_batches_subelement(root_element) + self._create_inactive_subelement(root_element) + self._create_generations_per_batch_subelement(root_element) + self._create_keff_trigger_subelement(root_element) self._create_source_subelement(root_element) self._create_output_subelement(root_element) self._create_statepoint_subelement(root_element) diff --git a/src/constants.F90 b/src/constants.F90 index 8ed657339..43ba2c223 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -424,7 +424,8 @@ module constants MODE_FIXEDSOURCE = 1, & ! Fixed source mode MODE_EIGENVALUE = 2, & ! K eigenvalue mode MODE_PLOTTING = 3, & ! Plotting mode - MODE_PARTICLE = 4 ! Particle restart mode + MODE_PARTICLE = 4, & ! Particle restart mode + MODE_VOLUME = 5 ! Volume calculation mode !============================================================================= ! CMFD CONSTANTS diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 2be8cbd57..6a2efa2e7 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -100,7 +100,6 @@ contains type(Node), pointer :: node_res_scat => null() type(Node), pointer :: node_scatterer => null() type(Node), pointer :: node_trigger => null() - type(Node), pointer :: node_keff_trigger => null() type(Node), pointer :: node_vol => null() type(Node), pointer :: node_tab_leg => null() type(NodeList), pointer :: node_scat_list => null() @@ -224,111 +223,49 @@ contains end if end if - ! Make sure that either eigenvalue or fixed source was specified - if (.not. check_for_node(doc, "eigenvalue") .and. & - .not. check_for_node(doc, "fixed_source")) then - call fatal_error(" or not specified.") - end if + if (check_for_node(doc, "run_mode")) then + call get_node_value(doc, "run_mode", temp_str) + select case (to_lower(temp_str)) + case ("eigenvalue") + run_mode = MODE_EIGENVALUE + case ("fixed source") + run_mode = MODE_FIXEDSOURCE + case ("plot") + run_mode = MODE_PLOTTING + case ("particle restart") + run_mode = MODE_PARTICLE + case ("volume") + run_mode = MODE_VOLUME + end select - ! Eigenvalue information - if (check_for_node(doc, "eigenvalue")) then - ! Set run mode - if (run_mode == NONE) run_mode = MODE_EIGENVALUE + ! Assume XML specifics , , etc. directly + node_mode => doc + else + call warning(" should be specified.") - ! Get pointer to eigenvalue XML block - call get_node_ptr(doc, "eigenvalue", node_mode) - - ! Check number of particles - if (.not. check_for_node(node_mode, "particles")) then - call fatal_error("Need to specify number of particles per generation.") + ! Make sure that either eigenvalue or fixed source was specified + if (.not. check_for_node(doc, "eigenvalue") .and. & + .not. check_for_node(doc, "fixed_source")) then + call fatal_error(" or not specified.") end if - ! Get number of particles - call get_node_value(node_mode, "particles", temp_long) + if (check_for_node(doc, "eigenvalue")) then + ! Set run mode + if (run_mode == NONE) run_mode = MODE_EIGENVALUE - ! If the number of particles was specified as a command-line argument, we - ! don't set it here - if (n_particles == 0) n_particles = temp_long + ! Get pointer to eigenvalue XML block + call get_node_ptr(doc, "eigenvalue", node_mode) + elseif (check_for_node(doc, "fixed_source")) then + ! Set run mode + if (run_mode == NONE) run_mode = MODE_FIXEDSOURCE - ! Get number of basic batches - call get_node_value(node_mode, "batches", n_batches) - if (.not. trigger_on) then - n_max_batches = n_batches - end if - - ! Get number of inactive batches - call get_node_value(node_mode, "inactive", n_inactive) - n_active = n_batches - n_inactive - if (check_for_node(node_mode, "generations_per_batch")) then - call get_node_value(node_mode, "generations_per_batch", gen_per_batch) - end if - - ! Allocate array for batch keff and entropy - allocate(k_generation(n_max_batches*gen_per_batch)) - allocate(entropy(n_max_batches*gen_per_batch)) - entropy = ZERO - - ! Get the trigger information for keff - if (check_for_node(node_mode, "keff_trigger")) then - call get_node_ptr(node_mode, "keff_trigger", node_keff_trigger) - - if (check_for_node(node_keff_trigger, "type")) then - call get_node_value(node_keff_trigger, "type", temp_str) - temp_str = trim(to_lower(temp_str)) - - select case (temp_str) - case ('std_dev') - keff_trigger % trigger_type = STANDARD_DEVIATION - case ('variance') - keff_trigger % trigger_type = VARIANCE - case ('rel_err') - keff_trigger % trigger_type = RELATIVE_ERROR - case default - call fatal_error("Unrecognized keff trigger type " // temp_str) - end select - - else - call fatal_error("Specify keff trigger type in settings XML") - end if - - if (check_for_node(node_keff_trigger, "threshold")) then - call get_node_value(node_keff_trigger, "threshold", & - keff_trigger % threshold) - else - call fatal_error("Specify keff trigger threshold in settings XML") - end if + ! Get pointer to fixed_source XML block + call get_node_ptr(doc, "fixed_source", node_mode) end if end if - ! Fixed source calculation information - if (check_for_node(doc, "fixed_source")) then - ! Set run mode - if (run_mode == NONE) run_mode = MODE_FIXEDSOURCE - - ! Get pointer to fixed_source XML block - call get_node_ptr(doc, "fixed_source", node_mode) - - ! Check number of particles - if (.not. check_for_node(node_mode, "particles")) then - call fatal_error("Need to specify number of particles per batch.") - end if - - ! Get number of particles - call get_node_value(node_mode, "particles", temp_long) - - ! If the number of particles was specified as a command-line argument, we - ! don't set it here - if (n_particles == 0) n_particles = temp_long - - ! Copy batch information - call get_node_value(node_mode, "batches", n_batches) - if (.not. trigger_on) then - n_max_batches = n_batches - end if - n_active = n_batches - n_inactive = 0 - gen_per_batch = 1 - end if + ! Read run parameters + call get_run_parameters(node_mode) ! Check number of active batches, inactive batches, and particles if (n_active <= 0) then @@ -1093,6 +1030,86 @@ contains end subroutine read_settings_xml +!=============================================================================== +! GET_RUN_PARAMETERS +!=============================================================================== + + subroutine get_run_parameters(node_base) + type(Node), pointer :: node_base + + integer(8) :: temp_long + character(MAX_LINE_LEN) :: temp_str + type(Node), pointer :: node_keff_trigger => null() + + ! Check number of particles + if (.not. check_for_node(node_base, "particles")) then + call fatal_error("Need to specify number of particles.") + end if + + ! Get number of particles + call get_node_value(node_base, "particles", temp_long) + + ! If the number of particles was specified as a command-line argument, we + ! don't set it here + if (n_particles == 0) n_particles = temp_long + + ! Get number of basic batches + call get_node_value(node_base, "batches", n_batches) + if (.not. trigger_on) then + n_max_batches = n_batches + end if + n_inactive = 0 + gen_per_batch = 1 + + ! Get number of inactive batches + if (run_mode == MODE_EIGENVALUE) then + call get_node_value(node_base, "inactive", n_inactive) + if (check_for_node(node_base, "generations_per_batch")) then + call get_node_value(node_base, "generations_per_batch", gen_per_batch) + end if + + ! Allocate array for batch keff and entropy + allocate(k_generation(n_max_batches*gen_per_batch)) + allocate(entropy(n_max_batches*gen_per_batch)) + entropy = ZERO + + ! Get the trigger information for keff + if (check_for_node(node_base, "keff_trigger")) then + call get_node_ptr(node_base, "keff_trigger", node_keff_trigger) + + if (check_for_node(node_keff_trigger, "type")) then + call get_node_value(node_keff_trigger, "type", temp_str) + temp_str = trim(to_lower(temp_str)) + + select case (temp_str) + case ('std_dev') + keff_trigger % trigger_type = STANDARD_DEVIATION + case ('variance') + keff_trigger % trigger_type = VARIANCE + case ('rel_err') + keff_trigger % trigger_type = RELATIVE_ERROR + case default + call fatal_error("Unrecognized keff trigger type " // temp_str) + end select + + else + call fatal_error("Specify keff trigger type in settings XML") + end if + + if (check_for_node(node_keff_trigger, "threshold")) then + call get_node_value(node_keff_trigger, "threshold", & + keff_trigger % threshold) + else + call fatal_error("Specify keff trigger threshold in settings XML") + end if + end if + end if + + ! Determine number of active batches + n_active = n_batches - n_inactive + + end subroutine get_run_parameters + !=============================================================================== ! READ_GEOMETRY_XML reads data from a geometry.xml file and parses it, checking ! for errors and placing properly-formatted data in the right data structures diff --git a/tests/test_asymmetric_lattice/inputs_true.dat b/tests/test_asymmetric_lattice/inputs_true.dat index ada4f8f05..32745dad1 100644 --- a/tests/test_asymmetric_lattice/inputs_true.dat +++ b/tests/test_asymmetric_lattice/inputs_true.dat @@ -205,11 +205,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -32 -32 0 32 32 32 diff --git a/tests/test_cmfd_feed/settings.xml b/tests/test_cmfd_feed/settings.xml index 41de07f56..eab90b70f 100644 --- a/tests/test_cmfd_feed/settings.xml +++ b/tests/test_cmfd_feed/settings.xml @@ -2,11 +2,10 @@ - - 20 - 10 - 1000 - + eigenvalue + 20 + 10 + 1000 @@ -19,7 +18,7 @@ - + 10 1 1 -10.0 -1.0 -1.0 diff --git a/tests/test_cmfd_nofeed/settings.xml b/tests/test_cmfd_nofeed/settings.xml index 41de07f56..eab90b70f 100644 --- a/tests/test_cmfd_nofeed/settings.xml +++ b/tests/test_cmfd_nofeed/settings.xml @@ -2,11 +2,10 @@ - - 20 - 10 - 1000 - + eigenvalue + 20 + 10 + 1000 @@ -19,7 +18,7 @@ - + 10 1 1 -10.0 -1.0 -1.0 diff --git a/tests/test_complex_cell/settings.xml b/tests/test_complex_cell/settings.xml index a6fd5da19..70b4e802f 100644 --- a/tests/test_complex_cell/settings.xml +++ b/tests/test_complex_cell/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_confidence_intervals/settings.xml b/tests/test_confidence_intervals/settings.xml index 19a27694e..09e53927d 100644 --- a/tests/test_confidence_intervals/settings.xml +++ b/tests/test_confidence_intervals/settings.xml @@ -3,11 +3,10 @@ true - - 10 - 2 - 100 - + eigenvalue + 10 + 2 + 100 diff --git a/tests/test_create_fission_neutrons/inputs_true.dat b/tests/test_create_fission_neutrons/inputs_true.dat index ad2e47f9b..b4245d2c8 100644 --- a/tests/test_create_fission_neutrons/inputs_true.dat +++ b/tests/test_create_fission_neutrons/inputs_true.dat @@ -18,10 +18,9 @@ - - 100 - 10 - + fixed source + 100 + 10 -1 -1 -1 1 1 1 diff --git a/tests/test_density/settings.xml b/tests/test_density/settings.xml index a6fd5da19..70b4e802f 100644 --- a/tests/test_density/settings.xml +++ b/tests/test_density/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_diff_tally/inputs_true.dat b/tests/test_diff_tally/inputs_true.dat index 0010c6d97..0d75e48c6 100644 --- a/tests/test_diff_tally/inputs_true.dat +++ b/tests/test_diff_tally/inputs_true.dat @@ -297,11 +297,10 @@ - - 100 - 3 - 0 - + eigenvalue + 100 + 3 + 0 -160 -160 -183 160 160 183 diff --git a/tests/test_distribmat/inputs_true.dat b/tests/test_distribmat/inputs_true.dat index d9b32ca84..746bb8bef 100644 --- a/tests/test_distribmat/inputs_true.dat +++ b/tests/test_distribmat/inputs_true.dat @@ -37,11 +37,10 @@ - - 1000 - 5 - 0 - + eigenvalue + 1000 + 5 + 0 -1 -1 -1 1 1 1 diff --git a/tests/test_eigenvalue_genperbatch/settings.xml b/tests/test_eigenvalue_genperbatch/settings.xml index a64b477b8..fefc2d059 100644 --- a/tests/test_eigenvalue_genperbatch/settings.xml +++ b/tests/test_eigenvalue_genperbatch/settings.xml @@ -1,12 +1,11 @@ - - 7 - 3 - 1000 - 3 - + eigenvalue + 7 + 3 + 1000 + 3 diff --git a/tests/test_eigenvalue_no_inactive/settings.xml b/tests/test_eigenvalue_no_inactive/settings.xml index 36d323ac0..f13a1665e 100644 --- a/tests/test_eigenvalue_no_inactive/settings.xml +++ b/tests/test_eigenvalue_no_inactive/settings.xml @@ -1,11 +1,10 @@ - - 10 - 0 - 1000 - + eigenvalue + 10 + 0 + 1000 diff --git a/tests/test_energy_cutoff/inputs_true.dat b/tests/test_energy_cutoff/inputs_true.dat index eafb2e389..7f67288d1 100644 --- a/tests/test_energy_cutoff/inputs_true.dat +++ b/tests/test_energy_cutoff/inputs_true.dat @@ -17,10 +17,9 @@ - - 100 - 10 - + fixed source + 100 + 10 -1 -1 -1 1 1 1 diff --git a/tests/test_energy_grid/settings.xml b/tests/test_energy_grid/settings.xml index 1e4b5937b..4b8d4fceb 100644 --- a/tests/test_energy_grid/settings.xml +++ b/tests/test_energy_grid/settings.xml @@ -3,11 +3,10 @@ 20000 - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_energy_laws/settings.xml b/tests/test_energy_laws/settings.xml index 1c3f444f0..946345eeb 100644 --- a/tests/test_energy_laws/settings.xml +++ b/tests/test_energy_laws/settings.xml @@ -1,10 +1,9 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_entropy/settings.xml b/tests/test_entropy/settings.xml index d6c6c4a47..df6a851ef 100644 --- a/tests/test_entropy/settings.xml +++ b/tests/test_entropy/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_filter_distribcell/case-1/settings.xml b/tests/test_filter_distribcell/case-1/settings.xml index 3ca5c1a32..14c6f3020 100644 --- a/tests/test_filter_distribcell/case-1/settings.xml +++ b/tests/test_filter_distribcell/case-1/settings.xml @@ -1,11 +1,10 @@ - - 1 - 0 - 1000 - + eigenvalue + 1 + 0 + 1000 diff --git a/tests/test_filter_distribcell/case-2/settings.xml b/tests/test_filter_distribcell/case-2/settings.xml index 3ca5c1a32..14c6f3020 100644 --- a/tests/test_filter_distribcell/case-2/settings.xml +++ b/tests/test_filter_distribcell/case-2/settings.xml @@ -1,11 +1,10 @@ - - 1 - 0 - 1000 - + eigenvalue + 1 + 0 + 1000 diff --git a/tests/test_filter_distribcell/case-3/settings.xml b/tests/test_filter_distribcell/case-3/settings.xml index 98c02fc30..ac716f320 100644 --- a/tests/test_filter_distribcell/case-3/settings.xml +++ b/tests/test_filter_distribcell/case-3/settings.xml @@ -1,11 +1,10 @@ - - 3 - 0 - 100 - + eigenvalue + 3 + 0 + 100 diff --git a/tests/test_filter_distribcell/case-4/settings.xml b/tests/test_filter_distribcell/case-4/settings.xml index 90e48b9c1..f3f0779bc 100644 --- a/tests/test_filter_distribcell/case-4/settings.xml +++ b/tests/test_filter_distribcell/case-4/settings.xml @@ -1,13 +1,12 @@ - - 1000 - 1 - 0 - - - - -1 -1 -1 1 1 1 - - + eigenvalue + 1000 + 1 + 0 + + + -1 -1 -1 1 1 1 + + diff --git a/tests/test_filter_energyfun/inputs_true.dat b/tests/test_filter_energyfun/inputs_true.dat index df627fa04..48f6c5033 100644 --- a/tests/test_filter_energyfun/inputs_true.dat +++ b/tests/test_filter_energyfun/inputs_true.dat @@ -298,11 +298,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -160 -160 -183 160 160 183 diff --git a/tests/test_filter_mesh/inputs_true.dat b/tests/test_filter_mesh/inputs_true.dat index 9f1951b56..1d9d6ac5c 100644 --- a/tests/test_filter_mesh/inputs_true.dat +++ b/tests/test_filter_mesh/inputs_true.dat @@ -297,11 +297,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -160 -160 -183 160 160 183 diff --git a/tests/test_fixed_source/settings.xml b/tests/test_fixed_source/settings.xml index 1e9b85d5a..e7e7f2f5b 100644 --- a/tests/test_fixed_source/settings.xml +++ b/tests/test_fixed_source/settings.xml @@ -1,10 +1,9 @@ - - 10 - 100 - + fixed source + 10 + 100 294 diff --git a/tests/test_infinite_cell/settings.xml b/tests/test_infinite_cell/settings.xml index a6fd5da19..70b4e802f 100644 --- a/tests/test_infinite_cell/settings.xml +++ b/tests/test_infinite_cell/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_iso_in_lab/inputs_true.dat b/tests/test_iso_in_lab/inputs_true.dat index 6b4b14c2c..682f8020f 100644 --- a/tests/test_iso_in_lab/inputs_true.dat +++ b/tests/test_iso_in_lab/inputs_true.dat @@ -297,11 +297,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -160 -160 -183 160 160 183 diff --git a/tests/test_lattice/settings.xml b/tests/test_lattice/settings.xml index d152920c1..ebe98a283 100644 --- a/tests/test_lattice/settings.xml +++ b/tests/test_lattice/settings.xml @@ -10,11 +10,10 @@ =============================================================== --> - - 10 - 5 - 100 - + eigenvalue + 10 + 5 + 100 diff --git a/tests/test_lattice_hex/settings.xml b/tests/test_lattice_hex/settings.xml index 810529cb1..0d87c3b86 100644 --- a/tests/test_lattice_hex/settings.xml +++ b/tests/test_lattice_hex/settings.xml @@ -1,10 +1,9 @@ - - 10 - 5 - 500 - + eigenvalue + 10 + 5 + 500 diff --git a/tests/test_lattice_mixed/settings.xml b/tests/test_lattice_mixed/settings.xml index b67824d03..69d50204e 100644 --- a/tests/test_lattice_mixed/settings.xml +++ b/tests/test_lattice_mixed/settings.xml @@ -1,10 +1,9 @@ - - 10 - 5 - 500 - + eigenvalue + 10 + 5 + 500 diff --git a/tests/test_lattice_multiple/settings.xml b/tests/test_lattice_multiple/settings.xml index 517637a59..569c80982 100644 --- a/tests/test_lattice_multiple/settings.xml +++ b/tests/test_lattice_multiple/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 100 - + eigenvalue + 10 + 5 + 100 diff --git a/tests/test_mg_basic/inputs_true.dat b/tests/test_mg_basic/inputs_true.dat index 182e5f46d..7141e57dd 100644 --- a/tests/test_mg_basic/inputs_true.dat +++ b/tests/test_mg_basic/inputs_true.dat @@ -84,11 +84,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 0.0 0.0 0.0 10.0 10.0 5.0 diff --git a/tests/test_mg_max_order/inputs_true.dat b/tests/test_mg_max_order/inputs_true.dat index 857e95bf5..a5feed722 100644 --- a/tests/test_mg_max_order/inputs_true.dat +++ b/tests/test_mg_max_order/inputs_true.dat @@ -30,11 +30,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 0.0 0.0 0.0 10.0 10.0 5.0 diff --git a/tests/test_mg_nuclide/inputs_true.dat b/tests/test_mg_nuclide/inputs_true.dat index 829b908ff..b5cb31b59 100644 --- a/tests/test_mg_nuclide/inputs_true.dat +++ b/tests/test_mg_nuclide/inputs_true.dat @@ -84,11 +84,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 0.0 0.0 0.0 10.0 10.0 5.0 diff --git a/tests/test_mg_tallies/inputs_true.dat b/tests/test_mg_tallies/inputs_true.dat index 56d97a394..244d92b06 100644 --- a/tests/test_mg_tallies/inputs_true.dat +++ b/tests/test_mg_tallies/inputs_true.dat @@ -84,11 +84,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 0.0 0.0 0.0 10.0 10.0 5.0 diff --git a/tests/test_mgxs_library_ce_to_mg/inputs_true.dat b/tests/test_mgxs_library_ce_to_mg/inputs_true.dat index 920dc4b5e..e31ee4d04 100644 --- a/tests/test_mgxs_library_ce_to_mg/inputs_true.dat +++ b/tests/test_mgxs_library_ce_to_mg/inputs_true.dat @@ -38,11 +38,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -0.63 -0.63 -1 0.63 0.63 1 diff --git a/tests/test_mgxs_library_condense/inputs_true.dat b/tests/test_mgxs_library_condense/inputs_true.dat index ac3f32c4f..ad6526fc1 100644 --- a/tests/test_mgxs_library_condense/inputs_true.dat +++ b/tests/test_mgxs_library_condense/inputs_true.dat @@ -38,11 +38,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -0.63 -0.63 -1 0.63 0.63 1 diff --git a/tests/test_mgxs_library_distribcell/inputs_true.dat b/tests/test_mgxs_library_distribcell/inputs_true.dat index fd0401527..9103a4b96 100644 --- a/tests/test_mgxs_library_distribcell/inputs_true.dat +++ b/tests/test_mgxs_library_distribcell/inputs_true.dat @@ -65,11 +65,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -10.71 -10.71 -1 10.71 10.71 1 diff --git a/tests/test_mgxs_library_hdf5/inputs_true.dat b/tests/test_mgxs_library_hdf5/inputs_true.dat index ac3f32c4f..ad6526fc1 100644 --- a/tests/test_mgxs_library_hdf5/inputs_true.dat +++ b/tests/test_mgxs_library_hdf5/inputs_true.dat @@ -38,11 +38,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -0.63 -0.63 -1 0.63 0.63 1 diff --git a/tests/test_mgxs_library_mesh/inputs_true.dat b/tests/test_mgxs_library_mesh/inputs_true.dat index 07b2076d7..421f4d8c7 100644 --- a/tests/test_mgxs_library_mesh/inputs_true.dat +++ b/tests/test_mgxs_library_mesh/inputs_true.dat @@ -297,11 +297,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -160 -160 -183 160 160 183 diff --git a/tests/test_mgxs_library_no_nuclides/inputs_true.dat b/tests/test_mgxs_library_no_nuclides/inputs_true.dat index ac3f32c4f..ad6526fc1 100644 --- a/tests/test_mgxs_library_no_nuclides/inputs_true.dat +++ b/tests/test_mgxs_library_no_nuclides/inputs_true.dat @@ -38,11 +38,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -0.63 -0.63 -1 0.63 0.63 1 diff --git a/tests/test_mgxs_library_nuclides/inputs_true.dat b/tests/test_mgxs_library_nuclides/inputs_true.dat index 05227ca8c..74dfc7b2e 100644 --- a/tests/test_mgxs_library_nuclides/inputs_true.dat +++ b/tests/test_mgxs_library_nuclides/inputs_true.dat @@ -38,11 +38,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -0.63 -0.63 -1 0.63 0.63 1 diff --git a/tests/test_multipole/inputs_true.dat b/tests/test_multipole/inputs_true.dat index 4b158b07c..2a1e865ad 100644 --- a/tests/test_multipole/inputs_true.dat +++ b/tests/test_multipole/inputs_true.dat @@ -34,11 +34,10 @@ - - 1000 - 5 - 0 - + eigenvalue + 1000 + 5 + 0 -1 -1 -1 1 1 1 diff --git a/tests/test_output/settings.xml b/tests/test_output/settings.xml index e2f3fc018..22b6267a3 100644 --- a/tests/test_output/settings.xml +++ b/tests/test_output/settings.xml @@ -3,11 +3,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_particle_restart_eigval/settings.xml b/tests/test_particle_restart_eigval/settings.xml index c01617e22..da37b4fc2 100644 --- a/tests/test_particle_restart_eigval/settings.xml +++ b/tests/test_particle_restart_eigval/settings.xml @@ -1,11 +1,10 @@ - - 12 - 5 - 1200 - + eigenvalue + 12 + 5 + 1200 diff --git a/tests/test_particle_restart_fixed/settings.xml b/tests/test_particle_restart_fixed/settings.xml index 1d8193f8f..9c731fde8 100644 --- a/tests/test_particle_restart_fixed/settings.xml +++ b/tests/test_particle_restart_fixed/settings.xml @@ -1,10 +1,9 @@ - - 12 - 1000 - + fixed source + 12 + 1000 diff --git a/tests/test_periodic/inputs_true.dat b/tests/test_periodic/inputs_true.dat index 764fec50d..5f7da5c47 100644 --- a/tests/test_periodic/inputs_true.dat +++ b/tests/test_periodic/inputs_true.dat @@ -25,11 +25,10 @@ - - 1000 - 4 - 0 - + eigenvalue + 1000 + 4 + 0 -5.0 -5.0 -5.0 5.0 5.0 5.0 diff --git a/tests/test_plot/settings.xml b/tests/test_plot/settings.xml index 03985b3ae..37623cb1f 100644 --- a/tests/test_plot/settings.xml +++ b/tests/test_plot/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_ptables_off/settings.xml b/tests/test_ptables_off/settings.xml index ab5404b20..5ae20fd38 100644 --- a/tests/test_ptables_off/settings.xml +++ b/tests/test_ptables_off/settings.xml @@ -3,11 +3,10 @@ false - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_quadric_surfaces/settings.xml b/tests/test_quadric_surfaces/settings.xml index 9f0e8ed05..81e5ad185 100644 --- a/tests/test_quadric_surfaces/settings.xml +++ b/tests/test_quadric_surfaces/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_reflective_plane/settings.xml b/tests/test_reflective_plane/settings.xml index a6fd5da19..70b4e802f 100644 --- a/tests/test_reflective_plane/settings.xml +++ b/tests/test_reflective_plane/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_resonance_scattering/inputs_true.dat b/tests/test_resonance_scattering/inputs_true.dat index 804c6c4f6..2801f115d 100644 --- a/tests/test_resonance_scattering/inputs_true.dat +++ b/tests/test_resonance_scattering/inputs_true.dat @@ -15,11 +15,10 @@ - - 1000 - 10 - 5 - + eigenvalue + 1000 + 10 + 5 -4 -4 -4 4 4 4 diff --git a/tests/test_rotation/settings.xml b/tests/test_rotation/settings.xml index a6fd5da19..70b4e802f 100644 --- a/tests/test_rotation/settings.xml +++ b/tests/test_rotation/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_salphabeta/settings.xml b/tests/test_salphabeta/settings.xml index a6fd5da19..70b4e802f 100644 --- a/tests/test_salphabeta/settings.xml +++ b/tests/test_salphabeta/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_score_current/settings.xml b/tests/test_score_current/settings.xml index 517637a59..569c80982 100644 --- a/tests/test_score_current/settings.xml +++ b/tests/test_score_current/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 100 - + eigenvalue + 10 + 5 + 100 diff --git a/tests/test_seed/settings.xml b/tests/test_seed/settings.xml index 0514e8a0a..11da44588 100644 --- a/tests/test_seed/settings.xml +++ b/tests/test_seed/settings.xml @@ -3,11 +3,10 @@ 239407351 - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_source/inputs_true.dat b/tests/test_source/inputs_true.dat index 727636c51..bb80609af 100644 --- a/tests/test_source/inputs_true.dat +++ b/tests/test_source/inputs_true.dat @@ -13,11 +13,10 @@ - - 1000 - 10 - 5 - + eigenvalue + 1000 + 10 + 5 diff --git a/tests/test_sourcepoint_batch/settings.xml b/tests/test_sourcepoint_batch/settings.xml index c816fe700..13096d551 100644 --- a/tests/test_sourcepoint_batch/settings.xml +++ b/tests/test_sourcepoint_batch/settings.xml @@ -4,11 +4,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_sourcepoint_latest/settings.xml b/tests/test_sourcepoint_latest/settings.xml index fa7f07dfa..58dfa671d 100644 --- a/tests/test_sourcepoint_latest/settings.xml +++ b/tests/test_sourcepoint_latest/settings.xml @@ -3,11 +3,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_sourcepoint_restart/settings.xml b/tests/test_sourcepoint_restart/settings.xml index 4e9b12d24..1a7a21357 100644 --- a/tests/test_sourcepoint_restart/settings.xml +++ b/tests/test_sourcepoint_restart/settings.xml @@ -4,11 +4,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_statepoint_batch/settings.xml b/tests/test_statepoint_batch/settings.xml index 0d8c3e88f..e2f8dad47 100644 --- a/tests/test_statepoint_batch/settings.xml +++ b/tests/test_statepoint_batch/settings.xml @@ -3,11 +3,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_statepoint_restart/settings.xml b/tests/test_statepoint_restart/settings.xml index ec9adfb1d..88382f74b 100644 --- a/tests/test_statepoint_restart/settings.xml +++ b/tests/test_statepoint_restart/settings.xml @@ -3,11 +3,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_statepoint_sourcesep/settings.xml b/tests/test_statepoint_sourcesep/settings.xml index 17d4ee2e2..86489bf33 100644 --- a/tests/test_statepoint_sourcesep/settings.xml +++ b/tests/test_statepoint_sourcesep/settings.xml @@ -4,11 +4,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_survival_biasing/settings.xml b/tests/test_survival_biasing/settings.xml index b0ff3fafc..6d5b66789 100644 --- a/tests/test_survival_biasing/settings.xml +++ b/tests/test_survival_biasing/settings.xml @@ -8,11 +8,10 @@ 1.2 - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_tallies/inputs_true.dat b/tests/test_tallies/inputs_true.dat index c349c3b94..88ff72f9b 100644 --- a/tests/test_tallies/inputs_true.dat +++ b/tests/test_tallies/inputs_true.dat @@ -297,11 +297,10 @@ - - 400 - 5 - 0 - + eigenvalue + 400 + 5 + 0 -160 -160 -183 160 160 183 diff --git a/tests/test_tally_aggregation/inputs_true.dat b/tests/test_tally_aggregation/inputs_true.dat index 23a35ad79..8223b8d6a 100644 --- a/tests/test_tally_aggregation/inputs_true.dat +++ b/tests/test_tally_aggregation/inputs_true.dat @@ -297,11 +297,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -160 -160 -183 160 160 183 diff --git a/tests/test_tally_arithmetic/inputs_true.dat b/tests/test_tally_arithmetic/inputs_true.dat index 323b56e49..4ade94e93 100644 --- a/tests/test_tally_arithmetic/inputs_true.dat +++ b/tests/test_tally_arithmetic/inputs_true.dat @@ -297,11 +297,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -160 -160 -183 160 160 183 diff --git a/tests/test_tally_assumesep/settings.xml b/tests/test_tally_assumesep/settings.xml index 517637a59..569c80982 100644 --- a/tests/test_tally_assumesep/settings.xml +++ b/tests/test_tally_assumesep/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 100 - + eigenvalue + 10 + 5 + 100 diff --git a/tests/test_tally_nuclides/settings.xml b/tests/test_tally_nuclides/settings.xml index b2ddb4248..32afc717a 100644 --- a/tests/test_tally_nuclides/settings.xml +++ b/tests/test_tally_nuclides/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 100 - + eigenvalue + 10 + 5 + 100 diff --git a/tests/test_tally_slice_merge/inputs_true.dat b/tests/test_tally_slice_merge/inputs_true.dat index 894173760..93b4af5cd 100644 --- a/tests/test_tally_slice_merge/inputs_true.dat +++ b/tests/test_tally_slice_merge/inputs_true.dat @@ -297,11 +297,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -160 -160 -183 160 160 183 diff --git a/tests/test_trace/settings.xml b/tests/test_trace/settings.xml index e50684eef..ce614711b 100644 --- a/tests/test_trace/settings.xml +++ b/tests/test_trace/settings.xml @@ -3,11 +3,10 @@ 5 1 453 - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_track_output/settings.xml b/tests/test_track_output/settings.xml index ef74341db..299ee72c5 100644 --- a/tests/test_track_output/settings.xml +++ b/tests/test_track_output/settings.xml @@ -1,19 +1,11 @@ - - - - - - 2 - 0 - 100 - - - - + eigenvalue + 2 + 0 + 100 @@ -26,5 +18,5 @@ 1 1 1 1 1 2 - + diff --git a/tests/test_translation/settings.xml b/tests/test_translation/settings.xml index a6fd5da19..70b4e802f 100644 --- a/tests/test_translation/settings.xml +++ b/tests/test_translation/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_trigger_batch_interval/settings.xml b/tests/test_trigger_batch_interval/settings.xml index b8e1e9c96..ecba8d347 100644 --- a/tests/test_trigger_batch_interval/settings.xml +++ b/tests/test_trigger_batch_interval/settings.xml @@ -1,14 +1,13 @@ - - 15 - 5 - 1000 - - std_dev - 0.004 - - + eigenvalue + 15 + 5 + 1000 + + std_dev + 0.004 + true diff --git a/tests/test_trigger_no_batch_interval/settings.xml b/tests/test_trigger_no_batch_interval/settings.xml index 5412af35f..9a4c94226 100644 --- a/tests/test_trigger_no_batch_interval/settings.xml +++ b/tests/test_trigger_no_batch_interval/settings.xml @@ -1,14 +1,13 @@ - - 15 - 5 - 1000 - - std_dev - 0.004 - - + eigenvalue + 15 + 5 + 1000 + + std_dev + 0.004 + true diff --git a/tests/test_trigger_no_status/settings.xml b/tests/test_trigger_no_status/settings.xml index 3d215a432..b85816240 100644 --- a/tests/test_trigger_no_status/settings.xml +++ b/tests/test_trigger_no_status/settings.xml @@ -1,20 +1,19 @@ - - 10 - 5 - 1000 - - std_dev - 0.009 - - + eigenvalue + 10 + 5 + 1000 + + std_dev + 0.009 + false 15 1 - + diff --git a/tests/test_trigger_tallies/settings.xml b/tests/test_trigger_tallies/settings.xml index 895c0ee72..d8f814bf3 100644 --- a/tests/test_trigger_tallies/settings.xml +++ b/tests/test_trigger_tallies/settings.xml @@ -1,14 +1,13 @@ - - 10 - 5 - 1000 - - std_dev - 0.001 - - + eigenvalue + 10 + 5 + 1000 + + std_dev + 0.001 + true diff --git a/tests/test_triso/inputs_true.dat b/tests/test_triso/inputs_true.dat index be6d2d0ea..43b15b037 100644 --- a/tests/test_triso/inputs_true.dat +++ b/tests/test_triso/inputs_true.dat @@ -430,11 +430,10 @@ - - 100 - 5 - 0 - + eigenvalue + 100 + 5 + 0 0.0 0.0 0.0 diff --git a/tests/test_uniform_fs/settings.xml b/tests/test_uniform_fs/settings.xml index d7956e743..3e6ef672f 100644 --- a/tests/test_uniform_fs/settings.xml +++ b/tests/test_uniform_fs/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_universe/settings.xml b/tests/test_universe/settings.xml index a6fd5da19..70b4e802f 100644 --- a/tests/test_universe/settings.xml +++ b/tests/test_universe/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 1000 - + eigenvalue + 10 + 5 + 1000 diff --git a/tests/test_void/settings.xml b/tests/test_void/settings.xml index b2ddb4248..32afc717a 100644 --- a/tests/test_void/settings.xml +++ b/tests/test_void/settings.xml @@ -1,11 +1,10 @@ - - 10 - 5 - 100 - + eigenvalue + 10 + 5 + 100 diff --git a/tests/test_volume_calc/inputs_true.dat b/tests/test_volume_calc/inputs_true.dat index af26cf8e7..28b4928cd 100644 --- a/tests/test_volume_calc/inputs_true.dat +++ b/tests/test_volume_calc/inputs_true.dat @@ -26,11 +26,10 @@ - - 1000 - 4 - 0 - + eigenvalue + 1000 + 4 + 0 -1.0 -1.0 -5.0 1.0 1.0 5.0 From 698efb04e21eb685e73b636d2369ac1c2437d061 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 16 Feb 2017 17:26:23 -0600 Subject: [PATCH 16/29] Allow volume calculations to be run alone --- src/input_xml.F90 | 25 +++++++++++++--------- src/main.F90 | 3 +++ src/simulation.F90 | 3 --- tests/test_volume_calc/inputs_true.dat | 10 +-------- tests/test_volume_calc/results_true.dat | 1 - tests/test_volume_calc/test_volume_calc.py | 19 ++++++---------- 6 files changed, 25 insertions(+), 36 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 6a2efa2e7..93bb1e733 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -264,16 +264,18 @@ contains end if end if - ! Read run parameters - call get_run_parameters(node_mode) + if (run_mode == MODE_EIGENVALUE .or. run_mode == MODE_FIXEDSOURCE) then + ! Read run parameters + call get_run_parameters(node_mode) - ! Check number of active batches, inactive batches, and particles - if (n_active <= 0) then - call fatal_error("Number of active batches must be greater than zero.") - elseif (n_inactive < 0) then - call fatal_error("Number of inactive batches must be non-negative.") - elseif (n_particles <= 0) then - call fatal_error("Number of particles must be greater than zero.") + ! Check number of active batches, inactive batches, and particles + if (n_active <= 0) then + call fatal_error("Number of active batches must be greater than zero.") + elseif (n_inactive < 0) then + call fatal_error("Number of inactive batches must be non-negative.") + elseif (n_particles <= 0) then + call fatal_error("Number of particles must be greater than zero.") + end if end if ! Copy random number seed if specified @@ -317,7 +319,10 @@ contains ! Get point to list of elements and make sure there is at least one call get_node_list(doc, "source", node_source_list) n = get_list_size(node_source_list) - if (n == 0) call fatal_error("No source specified in settings XML file.") + + if (run_mode == MODE_EIGENVALUE .or. run_mode == MODE_FIXEDSOURCE) then + if (n == 0) call fatal_error("No source specified in settings XML file.") + end if ! Allocate array for sources allocate(external_source(n)) diff --git a/src/main.F90 b/src/main.F90 index 8582ef703..1cd3a9e54 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -8,6 +8,7 @@ program main use particle_restart, only: run_particle_restart use plot, only: run_plot use simulation, only: run_simulation + use volume_calc, only: run_volume_calculations implicit none @@ -26,6 +27,8 @@ program main call run_plot() case (MODE_PARTICLE) if (master) call run_particle_restart() + case (MODE_VOLUME) + call run_volume_calculations() end select ! finalize run diff --git a/src/simulation.F90 b/src/simulation.F90 index 187632836..707f42737 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -39,9 +39,6 @@ contains type(Particle) :: p integer(8) :: i_work - ! Volume calculations - if (size(volume_calcs) > 0) call run_volume_calculations() - if (.not. restart_run) call initialize_source() ! Display header diff --git a/tests/test_volume_calc/inputs_true.dat b/tests/test_volume_calc/inputs_true.dat index 28b4928cd..627d30a40 100644 --- a/tests/test_volume_calc/inputs_true.dat +++ b/tests/test_volume_calc/inputs_true.dat @@ -26,15 +26,7 @@ - eigenvalue - 1000 - 4 - 0 - - - -1.0 -1.0 -5.0 1.0 1.0 5.0 - - + volume cell 1 2 3 diff --git a/tests/test_volume_calc/results_true.dat b/tests/test_volume_calc/results_true.dat index a36e61fb5..8eb2a61ac 100644 --- a/tests/test_volume_calc/results_true.dat +++ b/tests/test_volume_calc/results_true.dat @@ -1,4 +1,3 @@ -k-combined: 4.165450e-02 3.582533e-04 Volume calculation 0 Domain 1: 31.4693 +/- 0.0721 cm^3 Domain 2: 2.0933 +/- 0.0310 cm^3 diff --git a/tests/test_volume_calc/test_volume_calc.py b/tests/test_volume_calc/test_volume_calc.py index fa267efb6..cb4ecc2d7 100644 --- a/tests/test_volume_calc/test_volume_calc.py +++ b/tests/test_volume_calc/test_volume_calc.py @@ -52,22 +52,12 @@ class VolumeTest(PyAPITestHarness): # Define settings settings = openmc.Settings() - settings.particles = 1000 - settings.batches = 4 - settings.inactive = 0 - settings.source = openmc.Source(space=openmc.stats.Box( - [-1., -1., -5.], [1., 1., 5.])) + settings.run_mode = 'volume' settings.volume_calculations = vol_calcs settings.export_to_xml() def _get_results(self): - # Read the statepoint file. - statepoint = os.path.join(os.getcwd(), self._sp_name) - sp = openmc.StatePoint(statepoint) - - # Write out k-combined. - outstr = 'k-combined: {:12.6e} {:12.6e}\n'.format(*sp.k_combined) - + outstr = '' for i, filename in enumerate(sorted(glob.glob(os.path.join( os.getcwd(), 'volume_*.h5')))): outstr += 'Volume calculation {}\n'.format(i) @@ -83,6 +73,9 @@ class VolumeTest(PyAPITestHarness): return outstr + def _test_output_created(self): + pass + if __name__ == '__main__': - harness = VolumeTest('statepoint.4.h5') + harness = VolumeTest('') harness.main() From aeffb739ec82bed69f5a7efda44a4805c9bfc4bc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 17 Feb 2017 06:53:32 -0600 Subject: [PATCH 17/29] Make sure prism functions are in documentation --- docs/source/pythonapi/index.rst | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 1125a4b32..0f22d0102 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -115,15 +115,17 @@ Many of the above classes are derived from several abstract classes: openmc.Region openmc.Lattice -One function is also available to create a hexagonal region defined by the -intersection of six surface half-spaces. +Two helpers function are also available to create rectangular and hexagonal +prisms defined by the intersection of four and six surface half-spaces, +respectively. .. autosummary:: :toctree: generated :nosignatures: :template: myfunction.rst - openmc.make_hexagon_region + openmc.get_hexagonal_prism + openmc.get_rectangular_prism Constructing Tallies -------------------- From f28a0fcddaf0876b6e3b9c1cf1616795e8633db3 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 17 Feb 2017 09:29:14 -0600 Subject: [PATCH 18/29] Allow arbitrary MPI arguments in openmc.run() --- openmc/executor.py | 22 +++++++++---------- .../test_mgxs_library_ce_to_mg.py | 14 +++++------- .../test_statepoint_restart.py | 7 +++--- tests/testing_harness.py | 11 +++++----- 4 files changed, 25 insertions(+), 29 deletions(-) diff --git a/openmc/executor.py b/openmc/executor.py index a20437e1c..3b9b8abd8 100644 --- a/openmc/executor.py +++ b/openmc/executor.py @@ -44,8 +44,8 @@ def plot_geometry(output=True, openmc_exec='openmc', cwd='.'): def run(particles=None, threads=None, geometry_debug=False, - restart_file=None, tracks=False, mpi_procs=1, output=True, - openmc_exec='openmc', mpi_exec='mpiexec', cwd='.'): + restart_file=None, tracks=False, output=True, cwd='.', + openmc_exec='openmc', mpi_args=None): """Run an OpenMC simulation. Parameters @@ -56,23 +56,23 @@ def run(particles=None, threads=None, geometry_debug=False, Number of OpenMP threads. If OpenMC is compiled with OpenMP threading enabled, the default is implementation-dependent but is usually equal to the number of hardware threads available (or a value set by the - OMP_NUM_THREADS environment variable). + :envvar:`OMP_NUM_THREADS` environment variable). geometry_debug : bool, optional Turn on geometry debugging during simulation. Defaults to False. restart_file : str, optional Path to restart file to use tracks : bool, optional Write tracks for all particles. Defaults to False. - mpi_procs : int, optional - Number of MPI processes. output : bool, optional Capture OpenMC output from standard out. Defaults to True. + cwd : str, optional + Path to working directory to run in. Defaults to the current working + directory. openmc_exec : str, optional Path to OpenMC executable. Defaults to 'openmc'. - mpi_exec : str, optional - MPI execute command. Defaults to 'mpiexec'. - cwd : str, optional - Path to working directory to run in. Defaults to the current working directory. + mpi_args : list of str, optional + MPI execute command and any additional MPI arguments to pass, + e.g. ['mpiexec', '-n', '8']. """ @@ -94,8 +94,8 @@ def run(particles=None, threads=None, geometry_debug=False, if tracks: post_args += '-t' - if isinstance(mpi_procs, Integral) and mpi_procs > 1: - pre_args += '{} -n {} '.format(mpi_exec, mpi_procs) + if mpi_args is not None: + pre_args = ' '.join(mpi_args) + ' ' command = pre_args + openmc_exec + ' ' + post_args diff --git a/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py b/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py index 7dc679b1f..f9f9a7509 100644 --- a/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py +++ b/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py @@ -41,10 +41,9 @@ class MGXSTestHarness(PyAPITestHarness): def _run_openmc(self): # Initial run if self._opts.mpi_exec is not None: - returncode = openmc.run(mpi_procs=self._opts.mpi_np, - openmc_exec=self._opts.exe, - mpi_exec=self._opts.mpi_exec) - + mpi_args = [self._opts.mpi_exec, '-n', self._opts.mpi_np] + returncode = openmc.run(openmc_exec=self._opts.exe, + mpi_args=mpi_args) else: returncode = openmc.run(openmc_exec=self._opts.exe) @@ -74,10 +73,9 @@ class MGXSTestHarness(PyAPITestHarness): # Re-run MG mode. if self._opts.mpi_exec is not None: - returncode = openmc.run(mpi_procs=self._opts.mpi_np, - openmc_exec=self._opts.exe, - mpi_exec=self._opts.mpi_exec) - + mpi_args = [self._opts.mpi_exec, '-n', self._opts.mpi_np] + returncode = openmc.run(openmc_exec=self._opts.exe, + mpi_args=mpi_args) else: returncode = openmc.run(openmc_exec=self._opts.exe) diff --git a/tests/test_statepoint_restart/test_statepoint_restart.py b/tests/test_statepoint_restart/test_statepoint_restart.py index d39bf7cd5..61522ee05 100644 --- a/tests/test_statepoint_restart/test_statepoint_restart.py +++ b/tests/test_statepoint_restart/test_statepoint_restart.py @@ -50,11 +50,10 @@ class StatepointRestartTestHarness(TestHarness): # Run OpenMC if self._opts.mpi_exec is not None: - returncode = openmc.run(mpi_procs=self._opts.mpi_np, - restart_file=statepoint, + mpi_args = [self._opts.mpi_exec, '-n', self._opts.mpi_np] + returncode = openmc.run(restart_file=statepoint, openmc_exec=self._opts.exe, - mpi_exec=self._opts.mpi_exec) - + mpi_args=mpi_args) else: returncode = openmc.run(openmc_exec=self._opts.exe, restart_file=statepoint) diff --git a/tests/testing_harness.py b/tests/testing_harness.py index b833615c7..fc3b4f934 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -24,7 +24,7 @@ class TestHarness(object): self.parser = OptionParser() self.parser.add_option('--exe', dest='exe', default='openmc') self.parser.add_option('--mpi_exec', dest='mpi_exec', default=None) - self.parser.add_option('--mpi_np', dest='mpi_np', type=int, default=2) + self.parser.add_option('--mpi_np', dest='mpi_np', default='2') self.parser.add_option('--update', dest='update', action='store_true', default=False) self._opts = None @@ -62,9 +62,9 @@ class TestHarness(object): def _run_openmc(self): if self._opts.mpi_exec is not None: - returncode = openmc.run(mpi_procs=self._opts.mpi_np, - openmc_exec=self._opts.exe, - mpi_exec=self._opts.mpi_exec) + returncode = openmc.run( + openmc_exec=self._opts.exe, + mpi_args=[self._opts.mpi_exec, '-n', self._opts.mpi_np]) else: returncode = openmc.run(openmc_exec=self._opts.exe) @@ -190,8 +190,7 @@ class ParticleRestartTestHarness(TestHarness): # Set arguments args = {'openmc_exec': self._opts.exe} if self._opts.mpi_exec is not None: - args.update({'mpi_procs': self._opts.mpi_np, - 'mpi_exec': self._opts.mpi_exec}) + args['mpi_args'] = [self._opts.mpi_exec, '-n', self._opts.mpi_np] # Initial run returncode = openmc.run(**args) From 1fc810003433385fd9e534d0d669536644c905f1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 17 Feb 2017 09:43:40 -0600 Subject: [PATCH 19/29] Don't show run results for run modes where it doesn't make sense --- src/finalize.F90 | 45 ----------------------------- src/input_xml.F90 | 70 ++++++++++++++++++++++++---------------------- src/simulation.F90 | 56 +++++++++++++++++++++++++++++++++++-- 3 files changed, 90 insertions(+), 81 deletions(-) diff --git a/src/finalize.F90 b/src/finalize.F90 index 7ad99d0b5..0f1af2356 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -5,9 +5,6 @@ module finalize use global use hdf5_interface, only: hdf5_bank_t use message_passing - use output, only: print_runtime, print_results, & - print_overlap_check, write_tallies - use tally, only: tally_statistics implicit none @@ -22,30 +19,6 @@ contains integer :: hdf5_err - ! Start finalization timer - call time_finalize%start() - - if (run_mode /= MODE_PLOTTING .and. run_mode /= MODE_PARTICLE) then - ! Calculate statistics for tallies and write to tallies.out - if (master) then - if (n_realizations > 1) call tally_statistics() - end if - if (output_tallies) then - if (master) call write_tallies() - end if - if (check_overlaps) call reduce_overlap_count() - end if - - ! Stop timers and show timing statistics - call time_finalize%stop() - call time_total%stop() - if (master .and. (run_mode /= MODE_PLOTTING .and. & - run_mode /= MODE_PARTICLE)) then - call print_runtime() - call print_results() - if (check_overlaps) call print_overlap_check() - end if - ! Deallocate arrays call free_memory() @@ -65,22 +38,4 @@ contains end subroutine openmc_finalize -!=============================================================================== -! REDUCE_OVERLAP_COUNT accumulates cell overlap check counts to master -!=============================================================================== - - subroutine reduce_overlap_count() - -#ifdef MPI - if (master) then - call MPI_REDUCE(MPI_IN_PLACE, overlap_check_cnt, n_cells, & - MPI_INTEGER8, MPI_SUM, 0, mpi_intracomm, mpi_err) - else - call MPI_REDUCE(overlap_check_cnt, overlap_check_cnt, n_cells, & - MPI_INTEGER8, MPI_SUM, 0, mpi_intracomm, mpi_err) - end if -#endif - - end subroutine reduce_overlap_count - end module finalize diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 93bb1e733..46eea2c51 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -79,7 +79,6 @@ contains integer :: temp_int integer :: temp_int_array3(3) integer, allocatable :: temp_int_array(:) - integer(8) :: temp_long real(8), allocatable :: temp_real(:) integer :: n_tracks logical :: file_exists @@ -223,44 +222,47 @@ contains end if end if - if (check_for_node(doc, "run_mode")) then - call get_node_value(doc, "run_mode", temp_str) - select case (to_lower(temp_str)) - case ("eigenvalue") - run_mode = MODE_EIGENVALUE - case ("fixed source") - run_mode = MODE_FIXEDSOURCE - case ("plot") - run_mode = MODE_PLOTTING - case ("particle restart") - run_mode = MODE_PARTICLE - case ("volume") - run_mode = MODE_VOLUME - end select + ! Check run mode if it hasn't been set from the command line + if (run_mode == NONE) then + if (check_for_node(doc, "run_mode")) then + call get_node_value(doc, "run_mode", temp_str) + select case (to_lower(temp_str)) + case ("eigenvalue") + run_mode = MODE_EIGENVALUE + case ("fixed source") + run_mode = MODE_FIXEDSOURCE + case ("plot") + run_mode = MODE_PLOTTING + case ("particle restart") + run_mode = MODE_PARTICLE + case ("volume") + run_mode = MODE_VOLUME + end select - ! Assume XML specifics , , etc. directly - node_mode => doc - else - call warning(" should be specified.") + ! Assume XML specifics , , etc. directly + node_mode => doc + else + call warning(" should be specified.") - ! Make sure that either eigenvalue or fixed source was specified - if (.not. check_for_node(doc, "eigenvalue") .and. & - .not. check_for_node(doc, "fixed_source")) then - call fatal_error(" or not specified.") - end if + ! Make sure that either eigenvalue or fixed source was specified + if (.not. check_for_node(doc, "eigenvalue") .and. & + .not. check_for_node(doc, "fixed_source")) then + call fatal_error(" or not specified.") + end if - if (check_for_node(doc, "eigenvalue")) then - ! Set run mode - if (run_mode == NONE) run_mode = MODE_EIGENVALUE + if (check_for_node(doc, "eigenvalue")) then + ! Set run mode + if (run_mode == NONE) run_mode = MODE_EIGENVALUE - ! Get pointer to eigenvalue XML block - call get_node_ptr(doc, "eigenvalue", node_mode) - elseif (check_for_node(doc, "fixed_source")) then - ! Set run mode - if (run_mode == NONE) run_mode = MODE_FIXEDSOURCE + ! Get pointer to eigenvalue XML block + call get_node_ptr(doc, "eigenvalue", node_mode) + elseif (check_for_node(doc, "fixed_source")) then + ! Set run mode + if (run_mode == NONE) run_mode = MODE_FIXEDSOURCE - ! Get pointer to fixed_source XML block - call get_node_ptr(doc, "fixed_source", node_mode) + ! Get pointer to fixed_source XML block + call get_node_ptr(doc, "fixed_source", node_mode) + end if end if end if diff --git a/src/simulation.F90 b/src/simulation.F90 index 707f42737..1274c2b00 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -11,13 +11,15 @@ module simulation use global use message_passing use output, only: write_message, header, print_columns, & - print_batch_keff, print_generation + print_batch_keff, print_generation, print_runtime, & + print_results, print_overlap_check, write_tallies use particle_header, only: Particle use random_lcg, only: set_particle_seed use source, only: initialize_source, sample_external_source use state_point, only: write_state_point, write_source_point use string, only: to_str - use tally, only: synchronize_tallies, setup_active_usertallies + use tally, only: synchronize_tallies, setup_active_usertallies, & + tally_statistics use trigger, only: check_triggers use tracking, only: transport use volume_calc, only: run_volume_calculations @@ -110,6 +112,8 @@ contains if (master) call header("SIMULATION FINISHED", level=1) + call finalize_simulation() + ! Clear particle call p % clear() @@ -376,4 +380,52 @@ contains end subroutine replay_batch_history +!=============================================================================== +! FINALIZE_SIMULATION calculates tally statistics, writes tallies, and displays +! execution time and results +!=============================================================================== + + subroutine finalize_simulation + + ! Start finalization timer + call time_finalize%start() + + ! Calculate statistics for tallies and write to tallies.out + if (master) then + if (n_realizations > 1) call tally_statistics() + end if + if (output_tallies) then + if (master) call write_tallies() + end if + if (check_overlaps) call reduce_overlap_count() + + ! Stop timers and show timing statistics + call time_finalize%stop() + call time_total%stop() + if (master) then + call print_runtime() + call print_results() + if (check_overlaps) call print_overlap_check() + end if + + end subroutine finalize_simulation + +!=============================================================================== +! REDUCE_OVERLAP_COUNT accumulates cell overlap check counts to master +!=============================================================================== + + subroutine reduce_overlap_count() + +#ifdef MPI + if (master) then + call MPI_REDUCE(MPI_IN_PLACE, overlap_check_cnt, n_cells, & + MPI_INTEGER8, MPI_SUM, 0, mpi_intracomm, mpi_err) + else + call MPI_REDUCE(overlap_check_cnt, overlap_check_cnt, n_cells, & + MPI_INTEGER8, MPI_SUM, 0, mpi_intracomm, mpi_err) + end if +#endif + + end subroutine reduce_overlap_count + end module simulation From 46643382f477624e864ea0686fe480d794608490 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 17 Feb 2017 10:57:40 -0600 Subject: [PATCH 20/29] Update documentation --- docs/source/conf.py | 2 +- docs/source/usersguide/input.rst | 159 ++++++++++++++++--------------- 2 files changed, 82 insertions(+), 79 deletions(-) diff --git a/docs/source/conf.py b/docs/source/conf.py index d7fc23c43..22fc13a3f 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -248,7 +248,7 @@ napoleon_use_ivar = True intersphinx_mapping = { 'python': ('https://docs.python.org/3', None), - 'numpy': ('http://docs.scipy.org/doc/numpy/', None), + 'numpy': ('https://docs.scipy.org/doc/numpy/', None), 'pandas': ('http://pandas.pydata.org/pandas-docs/stable/', None), 'matplotlib': ('http://matplotlib.org/', None) } diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 809e02b91..fff479cb9 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -94,6 +94,18 @@ Settings Specification -- settings.xml All simulation parameters and miscellaneous options are specified in the settings.xml file. +```` Element +--------------------- + +The ```` element indicates the total number of batches to execute, +where each batch corresponds to a tally realization. In a fixed source +calculation, each batch consists of a number of source particles. In an +eigenvalue calculation, each batch consists of one or many fission source +iterations (generations), where each generation itself consists of a number of +source neutrons. + + *Default*: None + ```` Element ---------------------------------- @@ -132,67 +144,6 @@ you care. This element has the following attributes/sub-elements: *Default*: 0.0 -.. _eigenvalue: - -```` Element ------------------------- - -The ```` element indicates that a :math:`k`-eigenvalue calculation -should be performed. It has the following attributes/sub-elements: - - :batches: - The total number of batches, where each batch corresponds to multiple - fission source iterations. Batching is done to eliminate correlation between - realizations of random variables. - - *Default*: None - - :generations_per_batch: - The number of total fission source iterations per batch. - - *Default*: 1 - - :inactive: - The number of inactive batches. In general, the starting cycles in a - criticality calculation can not be used to contribute to tallies since the - fission source distribution and eigenvalue are generally not converged - immediately. - - *Default*: None - - :particles: - The number of neutrons to simulate per fission source iteration. - - *Default*: None - - :keff_trigger: - This tag specifies a precision trigger on the combined :math:`k_{eff}`. The - trigger is a convergence criterion on the uncertainty of the estimated - eigenvalue. It has the following attributes/sub-elements: - - :type: - The type of precision trigger. Accepted options are "variance", "std_dev", - and "rel_err". - - :variance: - Variance of the batch mean :math:`\sigma^2` - - :std_dev: - Standard deviation of the batch mean :math:`\sigma` - - :rel_err: - Relative error of the batch mean :math:`\frac{\sigma}{\mu}` - - *Default*: None - - :threshold: - The precision trigger's convergence criterion for the - combined :math:`k_{eff}`. - - *Default*: None - - .. note:: See section on the :ref:`trigger` for more information. - ```` Element ------------------------- @@ -247,23 +198,57 @@ problem. It has the following attributes/sub-elements: *Default*: None -```` Element +```` Element +----------------------------------- + +The ```` element indicates the number of total fission +source iterations per batch for an eigenvalue calculation. This element is +ignored for all run modes other than "eigenvalue". + + *Default*: 1 + +```` Element +---------------------- + +The ```` element indicates the number of inactive batches used in a +k-eigenvalue calculation. In general, the starting fission source iterations in +an eigenvalue calculation can not be used to contribute to tallies since the +fission source distribution and eigenvalue are generally not converged +immediately. + + *Default*: 0 + +```` Element -------------------------- -The ```` element indicates that a fixed source calculation should -be performed. It has the following attributes/sub-elements: +The ```` element specifies a precision trigger on the combined +:math:`k_{eff}`. The trigger is a convergence criterion on the uncertainty of +the estimated eigenvalue. It has the following attributes/sub-elements: - :batches: - The total number of batches. For fixed source calculations, each batch - represents a realization of random variables for tallies. + :type: + The type of precision trigger. Accepted options are "variance", "std_dev", + and "rel_err". + + :variance: + Variance of the batch mean :math:`\sigma^2` + + :std_dev: + Standard deviation of the batch mean :math:`\sigma` + + :rel_err: + Relative error of the batch mean :math:`\frac{\sigma}{\mu}` *Default*: None - :particles: - The number of particles to simulate per batch. + :threshold: + The precision trigger's convergence criterion for the + combined :math:`k_{eff}`. *Default*: None +.. note:: See section on the :ref:`trigger` for more information. + + ```` Element --------------------------- @@ -336,6 +321,15 @@ will abort. *Default*: Current working directory +```` Element +----------------------- + +This element indicates the number of neutrons to simulate per fission source +iteration when a k-eigenvalue calculation is performed or the number of neutrons +per batch for a fixed source simulation. + + *Default*: None + ```` Element --------------------- @@ -408,7 +402,16 @@ The ```` element indicates whether or not CMFD acceleration should be turned on or off. This element has no attributes or sub-elements and can be set to either "false" or "true". - *Defualt*: false + *Default*: false + +```` Element +---------------------- + +The ```` element indicates which run mode should be used when OpenMC +is executed. This element has no attributes or sub-elements and can be set to +"eigenvalue", "fixed source", "plot", "volume", or "particle restart". + + *Default*: None ```` Element ------------------ @@ -774,13 +777,13 @@ number, and particle number, respectively. ------------------------- OpenMC includes tally precision triggers which allow the user to define -uncertainty thresholds on :math:`k_{eff}` in the ```` subelement of -``settings.xml``, and/or tallies in ``tallies.xml``. When using triggers, +uncertainty thresholds on :math:`k_{eff}` in the ```` subelement +of ``settings.xml``, and/or tallies in ``tallies.xml``. When using triggers, OpenMC will run until it completes as many batches as defined by ````. -At this point, the uncertainties on all tallied values are computed and -compared with their corresponding trigger thresholds. If any triggers have not -been met, OpenMC will continue until either all trigger thresholds have been -satisfied or ```` has been reached. +At this point, the uncertainties on all tallied values are computed and compared +with their corresponding trigger thresholds. If any triggers have not been met, +OpenMC will continue until either all trigger thresholds have been satisfied or +```` has been reached. The ```` element provides an active "toggle switch" for tally precision trigger(s), the maximum number of batches and the batch interval. It @@ -793,8 +796,8 @@ has the following attributes/sub-elements: :max_batches: This describes the maximum number of batches allowed when using trigger(s). - .. note:: When max_batches is set, the number of ``batches`` shown in - ```` element represents minimum number of batches to + .. note:: When max_batches is set, the number of ``batches`` shown in the + ```` element represents minimum number of batches to simulate when using the trigger(s). :batch_interval: From 3298c247900011af8d453e28d8e1f8bf2661c490 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 17 Feb 2017 11:13:00 -0600 Subject: [PATCH 21/29] Fix variable declaration in timer_header --- src/timer_header.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/timer_header.F90 b/src/timer_header.F90 index 3a633f700..6f9785d56 100644 --- a/src/timer_header.F90 +++ b/src/timer_header.F90 @@ -42,11 +42,11 @@ contains function timer_get_value(self) result(elapsed) class(Timer), intent(in) :: self ! the timer - real(8) :: elapsed ! total elapsed time + real(8) :: elapsed ! total elapsed time integer(8) :: end_counts ! current number of counts integer(8) :: count_rate ! system-dependent counting rate - real :: elapsed_time ! elapsed time since last start + real(8) :: elapsed_time ! elapsed time since last start if (self % running) then call system_clock(end_counts, count_rate) From 606253e11cf795acc851bdabff8be04ea26206ea Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 17 Feb 2017 14:53:34 -0600 Subject: [PATCH 22/29] Write depletable attribute to XML and summary files --- openmc/material.py | 5 ++++- openmc/summary.py | 18 ++++++++++-------- src/input_xml.F90 | 7 +++++++ src/material_header.F90 | 3 ++- src/summary.F90 | 6 ++++++ 5 files changed, 29 insertions(+), 10 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 7c271cc30..7d59c52ee 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -255,7 +255,7 @@ class Material(object): @depletable.setter def depletable(self, depletable): - cv.check_type('Depletable flag for Material ID="{}"'.format(self._id), + cv.check_type('Depletable flag for Material ID="{}"'.format(self.id), depletable, bool) self._depletable = depletable @@ -771,6 +771,9 @@ class Material(object): if len(self._name) > 0: element.set("name", str(self._name)) + if self._depletable: + element.set("depletable", "true") + # Create temperature XML subelement if self.temperature is not None: subelement = ET.SubElement(element, "temperature") diff --git a/openmc/summary.py b/openmc/summary.py index b283ed502..5bc4c56a2 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -97,23 +97,25 @@ class Summary(object): # Values - Material objects self.materials = {} - for key in self._f['materials'].keys(): + for key, group in self._f['materials'].items(): if key == 'n_materials': continue material_id = int(key.lstrip('material ')) - index = self._f['materials'][key]['index'].value - name = self._f['materials'][key]['name'].value.decode() - density = self._f['materials'][key]['atom_density'].value - nuc_densities = self._f['materials'][key]['nuclide_densities'][...] - nuclides = self._f['materials'][key]['nuclides'].value + + index = group['index'].value + name = group['name'].value.decode() + density = group['atom_density'].value + nuc_densities = group['nuclide_densities'][...] + nuclides = group['nuclides'].value # Create the Material material = openmc.Material(material_id=material_id, name=name) + material.depletable = bool(group.attrs['depletable']) # Read the names of the S(a,b) tables for this Material and add them - if 'sab_names' in self._f['materials'][key]: - sab_tables = self._f['materials'][key]['sab_names'].value + if 'sab_names' in group: + sab_tables = group['sab_names'].value for sab_table in sab_tables: name = sab_table.decode() material.add_s_alpha_beta(name) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 46eea2c51..fa2e9881c 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2251,6 +2251,13 @@ contains call fatal_error("Must specify id of material in materials XML file") end if + ! Check if material is depletable + if (check_for_node(node_mat, "depletable")) then + call get_node_value(node_mat, "depletable", temp_str) + if (to_lower(temp_str) == "true" .or. temp_str == "1") & + mat % depletable = .true. + end if + ! Check to make sure 'id' hasn't been used if (material_dict % has_key(mat % id)) then call fatal_error("Two or more materials use the same unique ID: " & diff --git a/src/material_header.F90 b/src/material_header.F90 index f7c1b5db0..fdc36548f 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -31,8 +31,9 @@ module material_header character(20), allocatable :: names(:) ! isotope names character(20), allocatable :: sab_names(:) ! name of S(a,b) table - ! Does this material contain fissionable nuclides? + ! Does this material contain fissionable nuclides? Is it depletable? logical :: fissionable = .false. + logical :: depletable = .false. ! enforce isotropic scattering in lab logical, allocatable :: p0(:) diff --git a/src/summary.F90 b/src/summary.F90 index efe07b1ab..78a8276f1 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -536,6 +536,12 @@ contains material_group = create_group(materials_group, "material " // & trim(to_str(m%id))) + if (m % depletable) then + call write_attribute(material_group, "depletable", 1) + else + call write_attribute(material_group, "depletable", 0) + end if + ! Write internal OpenMC index for this material call write_dataset(material_group, "index", i) From e6a3c04817126a524549b82ed95dffdc2a6bdb8a Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 18 Feb 2017 05:20:30 -0500 Subject: [PATCH 23/29] docstring modifications and reduced histories in test_mg_convert in an attempt to bring its runtime down --- openmc/mgxs_library.py | 6 ++---- tests/test_mg_convert/inputs_true.dat | 6 +++--- tests/test_mg_convert/results_true.dat | 24 ++++++++++++------------ tests/test_mg_convert/test_mg_convert.py | 6 +++--- 4 files changed, 20 insertions(+), 22 deletions(-) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index d441e1bbb..deef7e82d 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1690,8 +1690,7 @@ class XSdata(object): azimuthal bins of an XSdata object that already uses an angular representation. Finally, this method simply uses an arithmetic mean to convert from an angular to isotropic representation; no flux-weighting - is applied and therefore the correctness of the solution is not - guaranteed. + is applied and therefore reaction rates will not be preserved. Parameters ---------- @@ -2453,8 +2452,7 @@ class MGXSLibrary(object): azimuthal bins of an XSdata object that already uses an angular representation. Finally, this method simply uses an arithmetic mean to convert from an angular to isotropic representation; no flux-weighting - is applied and therefore the correctness of the solution is not - guaranteed. + is applied and therefore the reaction rates will not be preserved. Parameters ---------- diff --git a/tests/test_mg_convert/inputs_true.dat b/tests/test_mg_convert/inputs_true.dat index 95ad61586..9bfb1dfba 100644 --- a/tests/test_mg_convert/inputs_true.dat +++ b/tests/test_mg_convert/inputs_true.dat @@ -17,9 +17,9 @@ - 1000 - 1000 - 100 + 100 + 10 + 5 diff --git a/tests/test_mg_convert/results_true.dat b/tests/test_mg_convert/results_true.dat index cf4777b8e..a52b929dc 100644 --- a/tests/test_mg_convert/results_true.dat +++ b/tests/test_mg_convert/results_true.dat @@ -1,24 +1,24 @@ k-combined: -9.957366E-01 7.836427E-05 +9.930873E-01 2.221904E-03 k-combined: -9.963227E-01 8.010344E-05 +9.948148E-01 1.216270E-03 k-combined: -9.957366E-01 7.836427E-05 +9.930873E-01 2.221904E-03 k-combined: -9.955152E-01 8.132242E-05 +9.755034E-01 6.178296E-03 k-combined: -9.953937E-01 8.452497E-05 +9.738059E-01 4.529068E-03 k-combined: -9.957744E-01 8.816057E-05 +9.866847E-01 9.485912E-03 k-combined: -9.955332E-01 8.894390E-05 +9.755024E-01 6.179047E-03 k-combined: -9.956463E-01 8.378533E-05 +9.738061E-01 4.529462E-03 k-combined: -9.955875E-01 8.818927E-05 +9.866835E-01 9.485832E-03 k-combined: -9.957017E-01 8.726677E-05 +9.719024E-01 4.213166E-03 k-combined: -9.957366E-01 7.836427E-05 +9.930873E-01 2.221904E-03 k-combined: -9.957366E-01 7.836427E-05 +9.930873E-01 2.221904E-03 diff --git a/tests/test_mg_convert/test_mg_convert.py b/tests/test_mg_convert/test_mg_convert.py index ad86e5a50..607195a82 100755 --- a/tests/test_mg_convert/test_mg_convert.py +++ b/tests/test_mg_convert/test_mg_convert.py @@ -11,9 +11,9 @@ from testing_harness import PyAPITestHarness import openmc # OpenMC simulation parameters -batches = 1000 -inactive = 100 -particles = 1000 +batches = 10 +inactive = 5 +particles = 100 def build_mgxs_library(convert): From ca746424e1646d7ef39d119c3127a2fa4583a866 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 18 Feb 2017 12:15:26 -0500 Subject: [PATCH 24/29] Closed StatePoint inbetween reads and writes in the test_mgxs_library_ce_to_mg; this was failing on my machine, but apparently not elsewhere. I thought I fixed it with a recompile of HDF5 which stopped the non-sensical error messages; turns out it came back and this is actually the culprit --- .../test_mgxs_library_ce_to_mg.py | 3 ++ tests/test_salphabeta/geometry.xml | 15 ------- tests/test_salphabeta/materials.xml | 41 ------------------- tests/test_salphabeta/settings.xml | 16 -------- 4 files changed, 3 insertions(+), 72 deletions(-) delete mode 100644 tests/test_salphabeta/geometry.xml delete mode 100644 tests/test_salphabeta/materials.xml delete mode 100644 tests/test_salphabeta/settings.xml diff --git a/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py b/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py index 7dc679b1f..e6511fd62 100644 --- a/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py +++ b/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py @@ -72,6 +72,9 @@ class MGXSTestHarness(PyAPITestHarness): if os.path.exists('./tallies.xml'): os.remove('./tallies.xml') + # Close the statepoint to allow writing + sp.close() + # Re-run MG mode. if self._opts.mpi_exec is not None: returncode = openmc.run(mpi_procs=self._opts.mpi_np, diff --git a/tests/test_salphabeta/geometry.xml b/tests/test_salphabeta/geometry.xml deleted file mode 100644 index 2b978b914..000000000 --- a/tests/test_salphabeta/geometry.xml +++ /dev/null @@ -1,15 +0,0 @@ - - - - - - - - - - - - - - - diff --git a/tests/test_salphabeta/materials.xml b/tests/test_salphabeta/materials.xml deleted file mode 100644 index bfe0a6224..000000000 --- a/tests/test_salphabeta/materials.xml +++ /dev/null @@ -1,41 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/tests/test_salphabeta/settings.xml b/tests/test_salphabeta/settings.xml deleted file mode 100644 index a6fd5da19..000000000 --- a/tests/test_salphabeta/settings.xml +++ /dev/null @@ -1,16 +0,0 @@ - - - - - 10 - 5 - 1000 - - - - - -4 -4 -4 4 4 4 - - - - From a7b174aff623825c1e07b5c0ddc54692f5fcffce Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 18 Feb 2017 12:16:58 -0500 Subject: [PATCH 25/29] Ahh!! I deleted the XML files! forgot that this was a legacy test --- tests/test_salphabeta/geometry.xml | 15 +++++++++++ tests/test_salphabeta/materials.xml | 41 +++++++++++++++++++++++++++++ tests/test_salphabeta/settings.xml | 16 +++++++++++ 3 files changed, 72 insertions(+) create mode 100644 tests/test_salphabeta/geometry.xml create mode 100644 tests/test_salphabeta/materials.xml create mode 100644 tests/test_salphabeta/settings.xml diff --git a/tests/test_salphabeta/geometry.xml b/tests/test_salphabeta/geometry.xml new file mode 100644 index 000000000..2b978b914 --- /dev/null +++ b/tests/test_salphabeta/geometry.xml @@ -0,0 +1,15 @@ + + + + + + + + + + + + + + + diff --git a/tests/test_salphabeta/materials.xml b/tests/test_salphabeta/materials.xml new file mode 100644 index 000000000..bfe0a6224 --- /dev/null +++ b/tests/test_salphabeta/materials.xml @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_salphabeta/settings.xml b/tests/test_salphabeta/settings.xml new file mode 100644 index 000000000..a6fd5da19 --- /dev/null +++ b/tests/test_salphabeta/settings.xml @@ -0,0 +1,16 @@ + + + + + 10 + 5 + 1000 + + + + + -4 -4 -4 4 4 4 + + + + From db33004a0f2c7db8eb8a67d92c09f25cc3a1ef93 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 18 Feb 2017 14:51:12 -0500 Subject: [PATCH 26/29] Added test for Legendre scattering - current tests converted all legenres to tabular distributions, this avoids that --- tests/test_mg_legendre/inputs_true.dat | 47 ++++++++++++++++++++++ tests/test_mg_legendre/results_true.dat | 2 + tests/test_mg_legendre/test_mg_legendre.py | 28 +++++++++++++ 3 files changed, 77 insertions(+) create mode 100644 tests/test_mg_legendre/inputs_true.dat create mode 100644 tests/test_mg_legendre/results_true.dat create mode 100644 tests/test_mg_legendre/test_mg_legendre.py diff --git a/tests/test_mg_legendre/inputs_true.dat b/tests/test_mg_legendre/inputs_true.dat new file mode 100644 index 000000000..fe0b5cf16 --- /dev/null +++ b/tests/test_mg_legendre/inputs_true.dat @@ -0,0 +1,47 @@ + + + + + + + + + + + + + + + + + ../1d_mgxs.h5 + + + + + + + + + + + + + + + + + 100 + 10 + 5 + + + + 0.0 0.0 0.0 10.0 10.0 5.0 + + + multi-group + + false + + diff --git a/tests/test_mg_legendre/results_true.dat b/tests/test_mg_legendre/results_true.dat new file mode 100644 index 000000000..d0f8c319e --- /dev/null +++ b/tests/test_mg_legendre/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +1.110122E+00 2.549637E-02 diff --git a/tests/test_mg_legendre/test_mg_legendre.py b/tests/test_mg_legendre/test_mg_legendre.py new file mode 100644 index 000000000..a16912c5d --- /dev/null +++ b/tests/test_mg_legendre/test_mg_legendre.py @@ -0,0 +1,28 @@ +#!/usr/bin/env python + +import os +import sys + +sys.path.insert(0, os.pardir) +from testing_harness import PyAPITestHarness +from input_set import MGInputSet + + +class MGMaxOrderTestHarness(PyAPITestHarness): + def __init__(self, statepoint_name, tallies_present, mg=False): + PyAPITestHarness.__init__(self, statepoint_name, tallies_present) + self._input_set = MGInputSet() + + def _build_inputs(self): + """Write input XML files.""" + reps = ['iso'] + self._input_set.build_default_materials_and_geometry(reps=reps) + self._input_set.build_default_settings() + # Enforce Legendre scattering + self._input_set.settings.tabular_legendre = {'enable': False} + self._input_set.export() + + +if __name__ == '__main__': + harness = MGMaxOrderTestHarness('statepoint.10.*', False, mg=True) + harness.main() From db744637536fb41830f74b1fb3db66272507f159 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 18 Feb 2017 14:39:56 -0600 Subject: [PATCH 27/29] Address @nelsonag comments on documentation --- docs/source/pythonapi/index.rst | 2 +- docs/source/usersguide/input.rst | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 0f22d0102..fb6cfe927 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -115,7 +115,7 @@ Many of the above classes are derived from several abstract classes: openmc.Region openmc.Lattice -Two helpers function are also available to create rectangular and hexagonal +Two helper function are also available to create rectangular and hexagonal prisms defined by the intersection of four and six surface half-spaces, respectively. diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index fff479cb9..a64454e92 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -214,14 +214,15 @@ The ```` element indicates the number of inactive batches used in a k-eigenvalue calculation. In general, the starting fission source iterations in an eigenvalue calculation can not be used to contribute to tallies since the fission source distribution and eigenvalue are generally not converged -immediately. +immediately. This element is ignored for all run modes other than "eigenvalue". *Default*: 0 ```` Element -------------------------- -The ```` element specifies a precision trigger on the combined +The ```` element (ignored for all run modes other than +"eigenvalue".) specifies a precision trigger on the combined :math:`k_{eff}`. The trigger is a convergence criterion on the uncertainty of the estimated eigenvalue. It has the following attributes/sub-elements: From 471ebbfa891fc7df8721969f0c16839315294c04 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 19 Feb 2017 06:10:20 -0500 Subject: [PATCH 28/29] Moved check for small values to end of convert_scatter_format and added blank line --- openmc/mgxs_library.py | 22 +++------------------- 1 file changed, 3 insertions(+), 19 deletions(-) diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index deef7e82d..af6dc8f98 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1837,6 +1837,7 @@ class XSdata(object): np.sum((l + 0.5) * eval_legendre(l, mu[imu]) * orig_data[..., l] for l in range(self.num_orders)) + elif target_format == 'histogram': # This code uses the vectorized integration capabilities # instead of having an isotropic and angle representation @@ -1856,10 +1857,6 @@ class XSdata(object): for l in range(self.num_orders)) new_data[..., h_bin] = simps(table_fine, mu_fine) - # Remove the very small values resulting from numerical - # precision issues - new_data[..., np.abs(new_data) < 1.E-10] = 0. - elif self.scatter_format == 'tabular': # Calculate the mu points of the current data mu_self = np.linspace(-1, 1, self.num_orders) @@ -1875,10 +1872,6 @@ class XSdata(object): for l in range(xsdata.num_orders): new_data[..., l] = simps(y[l], mu_fine) - # Remove the very small values resulting from numerical - # precision issues - new_data[..., np.abs(new_data) < 1.E-10] = 0. - elif target_format == 'tabular': # Simply use an interpolating function to get the new data mu = np.linspace(-1, 1, xsdata.num_orders) @@ -1898,10 +1891,6 @@ class XSdata(object): mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU) new_data[..., h_bin] = simps(interp(mu_fine), mu_fine) - # Remove the very small values resulting from numerical - # precision issues - new_data[..., np.abs(new_data) < 1.E-10] = 0. - elif self.scatter_format == 'histogram': # The histogram format does not have enough information to # convert to the other forms without inducing some amount of @@ -1929,10 +1918,6 @@ class XSdata(object): for l in range(xsdata.num_orders): new_data[..., l] = simps(y[l], mu_fine) - # Remove the very small values resulting from numerical - # precision issues - new_data[..., np.abs(new_data) < 1.E-10] = 0. - elif target_format == 'tabular': # Simply use an interpolating function to get the new data mu = np.linspace(-1, 1, xsdata.num_orders) @@ -1952,9 +1937,8 @@ class XSdata(object): new_data[..., h_bin] = \ norm * simps(interp(mu_fine), mu_fine) - # Remove the very small values resulting from numerical - # precision issues - new_data[..., np.abs(new_data) < 1.E-10] = 0. + # Remove small values resulting from numerical precision issues + new_data[..., np.abs(new_data) < 1.E-10] = 0. xsdata.set_scatter_matrix(new_data, temp) From 38ce9acdfbc13b5aac8b9b1370d71eb12525d4a8 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 19 Feb 2017 09:43:57 -0500 Subject: [PATCH 29/29] Updated openmc.run commands and input files to reflect the new run modes --- tests/test_mg_convert/inputs_true.dat | 9 ++++----- tests/test_mg_convert/test_mg_convert.py | 6 +++--- tests/test_mg_legendre/inputs_true.dat | 9 ++++----- 3 files changed, 11 insertions(+), 13 deletions(-) diff --git a/tests/test_mg_convert/inputs_true.dat b/tests/test_mg_convert/inputs_true.dat index 9bfb1dfba..85eeb3a24 100644 --- a/tests/test_mg_convert/inputs_true.dat +++ b/tests/test_mg_convert/inputs_true.dat @@ -16,11 +16,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 -5 -5 -5 5 5 5 diff --git a/tests/test_mg_convert/test_mg_convert.py b/tests/test_mg_convert/test_mg_convert.py index 607195a82..f4cd90e71 100755 --- a/tests/test_mg_convert/test_mg_convert.py +++ b/tests/test_mg_convert/test_mg_convert.py @@ -138,9 +138,9 @@ class MGXSTestHarness(PyAPITestHarness): build_mgxs_library(case) if self._opts.mpi_exec is not None: - returncode = openmc.run(mpi_procs=self._opts.mpi_np, - openmc_exec=self._opts.exe, - mpi_exec=self._opts.mpi_exec) + mpi_args = [self._opts.mpi_exec, '-n', self._opts.mpi_np] + returncode = openmc.run(openmc_exec=self._opts.exe, + mpi_args=mpi_args) else: returncode = openmc.run(openmc_exec=self._opts.exe) diff --git a/tests/test_mg_legendre/inputs_true.dat b/tests/test_mg_legendre/inputs_true.dat index fe0b5cf16..ca24a0dcf 100644 --- a/tests/test_mg_legendre/inputs_true.dat +++ b/tests/test_mg_legendre/inputs_true.dat @@ -30,11 +30,10 @@ - - 100 - 10 - 5 - + eigenvalue + 100 + 10 + 5 0.0 0.0 0.0 10.0 10.0 5.0