From 384a08956af7a92ae3fc70292b27aac698f091ca Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Thu, 16 Feb 2017 21:02:34 -0500 Subject: [PATCH] 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 e3a8ad4d85..d441e1bbbf 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 0000000000..9b0b0f4b17 --- /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 0000000000..5cf0eb6d26 --- /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 0000000000..c5c0ac16c7 --- /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()