diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index 9eb642e322..7bbd257ca0 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -252,7 +252,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 4fb347b73c..af6dc8f98c 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
@@ -14,10 +14,17 @@ 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 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
+
class XSdata(object):
"""A multi-group cross section data set providing all the
@@ -185,6 +192,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._representation = 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
@@ -318,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'
@@ -337,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)
@@ -348,14 +399,12 @@ class XSdata(object):
@representation.setter
def representation(self, representation):
- # Check it is of valid type.
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
@@ -369,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
@@ -384,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
@@ -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,269 @@ 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 (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 reaction rates will not be preserved.
+
+ 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)
+
+ # 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 ValueError("Cannot translate between `angle`"
+ " representations with different angle"
+ " bin structures")
+ # Nothing to do as the same structure was requested
+ return xsdata
+
+ xsdata.representation = target_representation
+ # We have different actions depending on the representation conversion
+ if target_representation == 'isotropic':
+ # 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
+
+ # 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 ['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 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)
+
+ 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 in self, but
+ represented as specified in :param:`target_format`.
+
+ """
+
+ from scipy.interpolate import interp1d
+ from scipy.integrate import 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 and re-generate XSdata.xs_shapes with the new scattering format
+ xsdata._xs_shapes = None
+
+ 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
+ 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)
+ # 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_fine = np.zeros(new_data.shape[:-1] + (_NMU,))
+ 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)
+
+ 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,
+ # 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)
+ for l in range(xsdata.num_orders)]
+ for l in range(xsdata.num_orders):
+ new_data[..., l] = simps(y[l], mu_fine)
+
+ 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)
+
+ elif self.scatter_format == 'histogram':
+ # The histogram format does not have enough information to
+ # 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.
+ # 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 = [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] = simps(y[l], mu_fine)
+
+ 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
+ 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.
+ for h_bin in range(xsdata.num_orders):
+ mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU)
+ new_data[..., h_bin] = \
+ norm * simps(interp(mu_fine), mu_fine)
+
+ # 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)
+
+ return xsdata
+
def to_hdf5(self, file):
"""Write XSdata to an HDF5 file
@@ -1757,7 +2067,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, :, :],
@@ -1995,6 +2305,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 +2427,75 @@ class MGXSLibrary(object):
result = xsdata
return result
+ 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 (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 reaction rates will not be preserved.
+
+ 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`.
+
+ """
+
+ 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 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 : {'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.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`.
+
+ """
+
+ 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.
diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90
index 3d9df0bbf0..511e2a2376 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
diff --git a/tests/test_mg_convert/inputs_true.dat b/tests/test_mg_convert/inputs_true.dat
new file mode 100644
index 0000000000..85eeb3a24a
--- /dev/null
+++ b/tests/test_mg_convert/inputs_true.dat
@@ -0,0 +1,29 @@
+
+
+
+
+
+
+
+
+
+
+ ./mgxs.h5
+
+
+
+
+
+
+
+ eigenvalue
+ 100
+ 10
+ 5
+
+
+ -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..a52b929dcf
--- /dev/null
+++ b/tests/test_mg_convert/results_true.dat
@@ -0,0 +1,24 @@
+k-combined:
+9.930873E-01 2.221904E-03
+k-combined:
+9.948148E-01 1.216270E-03
+k-combined:
+9.930873E-01 2.221904E-03
+k-combined:
+9.755034E-01 6.178296E-03
+k-combined:
+9.738059E-01 4.529068E-03
+k-combined:
+9.866847E-01 9.485912E-03
+k-combined:
+9.755024E-01 6.179047E-03
+k-combined:
+9.738061E-01 4.529462E-03
+k-combined:
+9.866835E-01 9.485832E-03
+k-combined:
+9.719024E-01 4.213166E-03
+k-combined:
+9.930873E-01 2.221904E-03
+k-combined:
+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
new file mode 100755
index 0000000000..f4cd90e713
--- /dev/null
+++ b/tests/test_mg_convert/test_mg_convert.py
@@ -0,0 +1,206 @@
+#!/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 = 10
+inactive = 5
+particles = 100
+
+
+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 = [['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]],
+ ['angle', 2], [['angle', 2], ['isotropic', None]]]
+
+ outstr = ''
+ for case in cases:
+ build_mgxs_library(case)
+
+ if self._opts.mpi_exec is not None:
+ 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)
+
+ 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()
diff --git a/tests/test_mg_legendre/inputs_true.dat b/tests/test_mg_legendre/inputs_true.dat
new file mode 100644
index 0000000000..ca24a0dcf7
--- /dev/null
+++ b/tests/test_mg_legendre/inputs_true.dat
@@ -0,0 +1,46 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ ../1d_mgxs.h5
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ eigenvalue
+ 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 0000000000..d0f8c319e7
--- /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 0000000000..a16912c5db
--- /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()
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 f9f9a75097..dbc2e0916c 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
@@ -71,6 +71,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:
mpi_args = [self._opts.mpi_exec, '-n', self._opts.mpi_np]