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()