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

Added test of convert_* methods;

Browse source
This commit is contained in:
Adam Nelson 2017-02-16 21:02:34 -05:00
parent 130db90c2a
commit 384a08956a
4 changed files with 303 additions and 54 deletions
Download patch file Download diff file Expand all files Collapse all files

92
openmc/mgxs_library.py
View file

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

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

@ -0,0 +1,30 @@
<?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>
<eigenvalue>
<particles>1000</particles>
<batches>2000</batches>
<inactive>100</inactive>
</eigenvalue>
<source strength="1.0">
<space type="box">
<parameters>-5 -5 -5 5 5 5</parameters>
</space>
</source>
<energy_mode>multi-group</energy_mode>
</settings>

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

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

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

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