crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #787 from nelsonag/convert

Browse source 
Modified openmc.MGXS to allow users to convert between representations and scattering formats
This commit is contained in:
Will Boyd 2017-02-19 11:30:21 -05:00 • committed by GitHub
commit bd5d00c9e3
10 changed files with 750 additions and 15 deletions
Download patch file Download diff file Expand all files Collapse all files

2
openmc/mgxs/mgxs.py
View file

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

419
openmc/mgxs_library.py
View file

@ -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.

6
src/scattdata_header.F90
View file

@ -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

29
tests/test_mg_convert/inputs_true.dat Normal file
View file

@ -0,0 +1,29 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="1" material="1" name="cell 1" region="4 -5 6 -7" universe="0" />
<surface boundary="reflective" coeffs="-5.0" id="4" name="left" type="x-plane" />
<surface boundary="vacuum" coeffs="5.0" id="5" name="right" type="x-plane" />
<surface boundary="reflective" coeffs="-5.0" id="6" name="bottom" type="y-plane" />
<surface boundary="reflective" coeffs="5.0" id="7" name="top" type="y-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>./mgxs.h5</cross_sections>
<material id="1" name="UO2 fuel">
<density units="macro" value="1.0" />
<macroscopic name="UO2" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>-5 -5 -5 5 5 5</parameters>
</space>
</source>
<energy_mode>multi-group</energy_mode>
</settings>

24
tests/test_mg_convert/results_true.dat Normal file
View file

@ -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

206
tests/test_mg_convert/test_mg_convert.py Executable file
View file

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

46
tests/test_mg_legendre/inputs_true.dat Normal file
View file

@ -0,0 +1,46 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<cell id="10000" material="10000" region="10000 -10001 10002 -10003 10004 -10005" universe="0" />
<cell id="10001" material="10001" region="10000 -10001 10002 -10003 10005 -10006" universe="0" />
<cell id="10002" material="10002" region="10000 -10001 10002 -10003 10006 -10007" universe="0" />
<surface boundary="reflective" coeffs="0.0" id="10000" type="x-plane" />
<surface boundary="reflective" coeffs="10.0" id="10001" type="x-plane" />
<surface boundary="reflective" coeffs="0.0" id="10002" type="y-plane" />
<surface boundary="reflective" coeffs="10.0" id="10003" type="y-plane" />
<surface boundary="reflective" coeffs="0.0" id="10004" type="z-plane" />
<surface coeffs="1.6667" id="10005" type="z-plane" />
<surface coeffs="3.3334" id="10006" type="z-plane" />
<surface boundary="reflective" coeffs="5.0" id="10007" type="z-plane" />
</geometry>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<cross_sections>../1d_mgxs.h5</cross_sections>
<material id="10000" name="1">
<density units="macro" value="1.0" />
<macroscopic name="uo2_iso" />
</material>
<material id="10001" name="2">
<density units="macro" value="1.0" />
<macroscopic name="clad_iso" />
</material>
<material id="10002" name="3">
<density units="macro" value="1.0" />
<macroscopic name="lwtr_iso" />
</material>
</materials>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<run_mode>eigenvalue</run_mode>
<particles>100</particles>
<batches>10</batches>
<inactive>5</inactive>
<source strength="1.0">
<space type="box">
<parameters>0.0 0.0 0.0 10.0 10.0 5.0</parameters>
</space>
</source>
<energy_mode>multi-group</energy_mode>
<tabular_legendre>
<enable>false</enable>
</tabular_legendre>
</settings>

2
tests/test_mg_legendre/results_true.dat Normal file
View file

@ -0,0 +1,2 @@
k-combined:
1.110122E+00 2.549637E-02

28
tests/test_mg_legendre/test_mg_legendre.py Normal file
View file

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

3
tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py
View file

@ -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]