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

Fixed array ordering of scattering matrix

Browse source
This commit is contained in:
Adam Nelson 2016-11-13 13:29:45 -05:00
parent 50a4079597
commit 2455388315
6 changed files with 84 additions and 51 deletions
Download patch file Download diff file Expand all files Collapse all files

25
openmc/mgxs/library.py
View file

@ -287,7 +287,7 @@ class Library(object):
def by_nuclide(self, by_nuclide):
cv.check_type('by_nuclide', by_nuclide, bool)
if by_nuclide == True and self.domain_type == 'mesh':
if by_nuclide and self.domain_type == 'mesh':
raise ValueError('Unable to create MGXS library by nuclide with '
'mesh domain')
@ -297,7 +297,7 @@ class Library(object):
def domain_type(self, domain_type):
cv.check_value('domain type', domain_type, openmc.mgxs.DOMAIN_TYPES)
if self.by_nuclide == True and domain_type == 'mesh':
if self.by_nuclide and domain_type == 'mesh':
raise ValueError('Unable to create MGXS library by nuclide with '
'mesh domain')
@ -372,7 +372,15 @@ class Library(object):
@scatter_format.setter
def scatter_format(self, scatter_format):
cv.check_value('scatter_format', scatter_format, openmc.mgxs.MU_TREATMENTS)
cv.check_value('scatter_format', scatter_format,
openmc.mgxs.MU_TREATMENTS)
if scatter_format == 'histogram' and self.correction == 'P0':
msg = 'The P0 correction will be ignored since the ' \
'scatter format is set to histogram'
warn(msg)
self.correction = None
self._scatter_format = scatter_format
@legendre_order.setter
@ -405,6 +413,13 @@ class Library(object):
msg = 'The histogram bins will be ignored since the ' \
'scatter format is set to legendre'
warn(msg)
elif self.scatter_format == 'histogram':
if self.correction == 'P0':
msg = 'The P0 correction will be ignored since ' \
'a histogram representation of the scattering '\
'kernel is requested'
warn(msg, RuntimeWarning)
self.correction = None
self._histogram_bins = histogram_bins
@ -471,7 +486,7 @@ class Library(object):
mgxs.delayed_groups = None
else:
delayed_groups \
= list(range(1,self.num_delayed_groups+1))
= list(range(1, self.num_delayed_groups + 1))
mgxs.delayed_groups = delayed_groups
# If a tally trigger was specified, add it to the MGXS
@ -1351,7 +1366,7 @@ class Library(object):
'scattering matrix is not provided.')
# Total or transport can be present, but if using
# self.correction=="P0", then we should use transport.
if (((self.correction is "P0") and
if (((self.correction == "P0") and
('nu-transport' not in self.mgxs_types))):
error_flag = True
warn('A "nu-transport" MGXS type is required since a "P0" '

17
openmc/mgxs/mgxs.py
View file

@ -3775,15 +3775,16 @@ class ScatterMatrixXS(MatrixMGXS):
# tally data is stored in order of increasing energies
if order_groups == 'increasing':
xs = xs[:, ::-1, ::-1, ...]
# Place the histogram bins before the incoming groups
if self.scatter_format == 'histogram':
xs = np.swapaxes(xs, 3, 2)
xs = np.swapaxes(xs, 2, 1)
# import pdb; pdb.set_trace()
if squeeze:
xs = np.squeeze(xs)
xs = np.atleast_2d(xs)
# We want to squeeze out everything but the in_groups, out_groups,
# and, if needed, num_mu_bins dimension. These must not be squeezed
# so 1-group problems have the correct shape.
if self.scatter_format == 'histogram':
axes = (5, 4, 0)
else:
axes = (4, 3, 0)
xs = np.squeeze(xs, axis=axes)
return xs

43
openmc/mgxs_library.py
View file

@ -15,7 +15,7 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \
# Supported incoming particle MGXS angular treatment representations
_REPRESENTATIONS = ['isotropic', 'angle']
_SCATTER_TYPES = ['tabular', 'legendre', 'histogram']
_XS_SHAPES = ["[Order][G][G']", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]",
_XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]",
"[G'][DG]"]
@ -134,7 +134,7 @@ class XSdata(object):
Note that some cross sections can be input in more than one shape so they
are listed multiple times:
[Order][G][G']: scatter_matrix
[G][G'][Order]: scatter_matrix
[G]: total, absorption, fission, kappa_fission, nu_fission,
prompt_nu_fission, inverse_velocity
@ -298,9 +298,9 @@ class XSdata(object):
self.num_delayed_groups)
self._xs_shapes["[G'][DG]"] = (self.energy_groups.num_groups,
self.num_delayed_groups)
self._xs_shapes["[Order][G][G']"] \
= (self.num_orders, self.energy_groups.num_groups,
self.energy_groups.num_groups)
self._xs_shapes["[G][G'][Order]"] \
= (self.energy_groups.num_groups,
self.energy_groups.num_groups, self.num_orders)
# If representation is by angle prepend num polar and num azim
if self.representation == 'angle':
@ -718,7 +718,7 @@ class XSdata(object):
"""
# Get the accepted shapes for this xs
shapes = [self.xs_shapes["[Order][G][G']"]]
shapes = [self.xs_shapes["[G][G'][Order]"]]
# Convert to a numpy array so we can easily get the shape for checking
scatter = np.asarray(scatter)
@ -1534,11 +1534,10 @@ class XSdata(object):
if self.representation == 'isotropic':
if self.scatter_format == 'legendre':
# Get the scattering orders in the outermost dimension
self._scatter_matrix[i] = np.zeros((self.num_orders,
self.energy_groups.num_groups,
self.energy_groups.num_groups))
self._scatter_matrix[i] = \
np.zeros(self.xs_shapes["[G][G'][Order]"])
for moment in range(self.num_orders):
self._scatter_matrix[i][moment, :, :] = \
self._scatter_matrix[i][:, :, moment] = \
scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
moment=moment, subdomains=subdomain)
else:
@ -1657,7 +1656,7 @@ class XSdata(object):
if self.num_polar is not None:
grp.attrs['num_polar'] = self.num_polar
grp.attrs['scatter_shape'] = np.string_("[Order][G][G']")
grp.attrs['scatter_shape'] = np.string_("[G][G'][Order]")
if self.scatter_format is not None:
grp.attrs['scatter_format'] = np.string_(self.scatter_format)
if self.order is not None:
@ -1736,14 +1735,13 @@ class XSdata(object):
# Get the sparse scattering data to print to the library
G = self.energy_groups.num_groups
if self.representation == 'isotropic':
g_out_bounds = np.zeros((G, 2), dtype=np.int)
for g_in in range(G):
if self.scatter_format == 'legendre':
nz = np.nonzero(self._scatter_matrix[i][0, g_in, :])
nz = np.nonzero(self._scatter_matrix[i][g_in, :, 0])
elif self.scatter_format == 'histogram':
nz = np.nonzero(np.sum(
self._scatter_matrix[i][:, g_in, :], axis=0))
self._scatter_matrix[i][g_in, :, :], axis=1))
if len(nz[0]) == 0:
g_out_bounds[g_in, 0] = 0
g_out_bounds[g_in, 1] = 0
@ -1757,8 +1755,8 @@ class XSdata(object):
for g_in in range(G):
for g_out in range(g_out_bounds[g_in, 0],
g_out_bounds[g_in, 1] + 1):
for l in range(len(matrix[:, g_in, g_out])):
flat_scatt.append(matrix[l, g_in, g_out])
for l in range(len(matrix[g_in, g_out, :])):
flat_scatt.append(matrix[g_in, g_out, l])
# And write it.
scatt_grp = xs_grp.create_group('scatter_data')
@ -1767,7 +1765,6 @@ class XSdata(object):
# Repeat for multiplicity
if self._multiplicity_matrix[i] is not None:
# Now create the flattened scatter matrix array
matrix = self._multiplicity_matrix[i][:, :]
flat_mult = []
@ -1793,7 +1790,13 @@ class XSdata(object):
for p in range(Np):
for a in range(Na):
for g_in in range(G):
matrix = self._scatter_matrix[i][p, a, 0, g_in, :]
if self.scatter_format == 'legendre':
matrix = \
self._scatter_matrix[i][p, a, g_in, :, 0]
elif self.scatter_format == 'histogram':
matrix = \
np.sum(self._scatter_matrix[i][p, a, g_in, :, :],
axis=1)
nz = np.nonzero(matrix)
g_out_bounds[p, a, g_in, 0] = nz[0][0]
g_out_bounds[p, a, g_in, 1] = nz[0][-1]
@ -1806,8 +1809,8 @@ class XSdata(object):
for g_in in range(G):
for g_out in range(g_out_bounds[p, a, g_in, 0],
g_out_bounds[p, a, g_in, 1] + 1):
for l in range(len(matrix[:, g_in, g_out])):
flat_scatt.append(matrix[l, g_in, g_out])
for l in range(len(matrix[g_in, g_out, :])):
flat_scatt.append(matrix[g_in, g_out, l])
# And write it.
scatt_grp = xs_grp.create_group('scatter_data')

40
scripts/openmc-update-mgxs
View file

@ -23,12 +23,10 @@ optional arguments:
from __future__ import print_function
import os
from shutil import move
import warnings
import xml.etree.ElementTree as ET
import argparse
import h5py
import numpy as np
import openmc.mgxs_library
@ -51,7 +49,7 @@ def parse_args():
if args['output'] == '':
filename = args['input'].name
extension = filenameos.path.splitext()
extension = os.path.splitext(filename)
if extension == '.xml':
filename = filename[:filename.rfind('.')] + '.h5'
args['output'] = filename
@ -84,6 +82,8 @@ if __name__ == '__main__':
temp = tree.find('group_structure').text.strip()
temp = np.array(temp.split())
group_structure = temp.astype(np.float)
# Convert from MeV to eV
group_structure *= 1.e6
energy_groups = openmc.mgxs.EnergyGroups(group_structure)
temp = tree.find('inverse-velocity')
if temp is not None:
@ -103,7 +103,7 @@ if __name__ == '__main__':
temperature = get_data(xsdata_elem, 'kT')
if temperature is not None:
temperature = \
float(temperature) / openmc.data.K_BOLTZMANN
float(temperature) / openmc.data.K_BOLTZMANN * 1.E6
else:
temperature = 294.
temperatures = [temperature]
@ -163,7 +163,7 @@ if __name__ == '__main__':
if temp is not None:
temp = np.array(temp.split(), dtype=float)
total = temp.astype(np.float)
total.shape = xsd[i].vector_shape
total.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_total(total, temperature)
if inverse_velocity is not None:
@ -172,41 +172,55 @@ if __name__ == '__main__':
temp = get_data(xsdata_elem, 'absorption')
temp = np.array(temp.split())
absorption = temp.astype(np.float)
absorption.shape = xsd[i].vector_shape
absorption.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_absorption(absorption, temperature)
temp = get_data(xsdata_elem, 'scatter')
temp = np.array(temp.split())
scatter = temp.astype(np.float)
scatter.shape = xsd[i].pn_matrix_shape
# This is now a flattened-array of something that started with a
# shape of [Order][G][G']; we need to unflatten and then switch the
# ordering
in_shape = (order_dim, energy_groups.num_groups,
energy_groups.num_groups)
if representation == 'angle':
in_shape = (n_pol, n_azi) + in_shape
scatter.shape = in_shape
scatter = np.swapaxes(scatter, 2, 3)
scatter = np.swapaxes(scatter, 3, 4)
else:
scatter.shape = in_shape
scatter = np.swapaxes(scatter, 0, 1)
scatter = np.swapaxes(scatter, 1, 2)
xsd[i].set_scatter_matrix(scatter, temperature)
temp = get_data(xsdata_elem, 'multiplicity')
if temp is not None:
temp = np.array(temp.split())
multiplicity = temp.astype(np.float)
multiplicity.shape = xsd[i].matrix_shape
multiplicity.shape = xsd[i].xs_shapes["[G][G']"]
xsd[i].set_multiplicity_matrix(multiplicity, temperature)
temp = get_data(xsdata_elem, 'fission')
if temp is not None:
temp = np.array(temp.split())
fission = temp.astype(np.float)
fission.shape = xsd[i].vector_shape
fission.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_fission(fission, temperature)
temp = get_data(xsdata_elem, 'kappa_fission')
if temp is not None:
temp = np.array(temp.split())
kappa_fission = temp.astype(np.float)
kappa_fission.shape = xsd[i].vector_shape
kappa_fission.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_kappa_fission(kappa_fission, temperature)
temp = get_data(xsdata_elem, 'chi')
if temp is not None:
temp = np.array(temp.split())
chi = temp.astype(np.float)
chi.shape = xsd[i].vector_shape
chi.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_chi(chi, temperature)
else:
chi = None
@ -216,9 +230,9 @@ if __name__ == '__main__':
temp = np.array(temp.split())
nu_fission = temp.astype(np.float)
if chi is not None:
nu_fission.shape = xsd[i].vector_shape
nu_fission.shape = xsd[i].xs_shapes['[G]']
else:
nu_fission.shape = xsd[i].matrix_shape
nu_fission.shape = xsd[i].xs_shapes["[G][G']"]
xsd[i].set_nu_fission(nu_fission, temperature)
# Build library as we go, but first we have enough to initialize it

6
src/mgxs_header.F90
View file

@ -355,7 +355,7 @@ module mgxs_header
if (attribute_exists(xs_id, "scatter_shape")) then
call read_attribute(temp_str, xs_id, "scatter_shape")
temp_str = trim(temp_str)
if (to_lower(temp_str) /= "[order][g][g']") then
if (to_lower(temp_str) /= "[g][g'][order]") then
call fatal_error("Invalid scatter_shape option!")
end if
end if
@ -1794,8 +1794,8 @@ module mgxs_header
order_dim = order + 1
end if
! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for passing
! to ScattData
! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for
! passing to ScattData
allocate(input_scatt(energy_groups, this % n_azi, this % n_pol))
index = 1

4
src/scattdata_header.F90
View file

@ -274,7 +274,7 @@ contains
! Get this by summing the un-normalized P0 coefficient in matrix
! over all outgoing groups
do gin = 1, groups
this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1)
this % scattxs(gin) = sum(matrix(gin) % data(:, :))
end do
allocate(energy(groups))
@ -282,7 +282,7 @@ contains
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
do gin = 1, groups
allocate(energy(gin) % data(gmin(gin):gmax(gin)))
do gout = 1, groups
do gout = gmin(gin), gmax(gin)
norm = sum(matrix(gin) % data(:, gout))
energy(gin) % data(gout) = norm
if (norm /= ZERO) then