crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Replace Chain.index_reaction with Chain.reactions.

Browse source 
Now the ReactionRates controls all the indexing needed to build a depletion
matrix. Got all tests fixed here too. Decided to add get/set methods on
ReactionRates which take strings.
This commit is contained in:
Paul Romano 2018-02-17 14:31:04 -06:00
parent c9abcfc0a1
commit 5c4ea0d640
7 changed files with 147 additions and 113 deletions
Download patch file Download diff file Expand all files Collapse all files

27
openmc/deplete/chain.py
View file

@ -113,19 +113,19 @@ class Chain(object):
Attributes
----------
nuclides : list of Nuclide
List of nuclides in chain.
nuclides : list of openmc.deplete.Nuclide
Nuclides present in the chain.
reactions : list of str
Reactions that are tracked in the depletion chain
nuclide_dict : OrderedDict of str to int
Maps a nuclide name to an index in nuclides.
index_reaction : OrderedDict of str to int
Dictionary mapping a reaction name to an index in ReactionRates.
"""
def __init__(self):
self.nuclides = []
self.reactions = []
self.nuclide_dict = OrderedDict()
self.index_reaction = OrderedDict()
def __contains__(self, nuclide):
return nuclide in self.nuclide_dict
@ -190,7 +190,6 @@ class Chain(object):
missing_fpy = []
missing_fp = []
reaction_index = 0
for idx, parent in enumerate(sorted(decay_data, key=_get_zai)):
data = decay_data[parent]
@ -232,9 +231,8 @@ class Chain(object):
Z = data.nuclide['atomic_number'] + delta_Z
daughter = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
if name not in chain.index_reaction:
chain.index_reaction[name] = reaction_index
reaction_index += 1
if name not in chain.reactions:
chain.reactions.append(name)
if daughter not in decay_data:
missing_rx_product.append((parent, name, daughter))
@ -256,9 +254,8 @@ class Chain(object):
nuclide.reactions.append(
ReactionTuple('fission', 0, q_value, 1.0))
if 'fission' not in chain.index_reaction:
chain.index_reaction['fission'] = reaction_index
reaction_index += 1
if 'fission' not in chain.reactions:
chain.reactions.append('fission')
else:
missing_fpy.append(parent)
@ -340,16 +337,14 @@ class Chain(object):
print('Decay chain "', filename, '" is invalid.')
raise
reaction_index = 0
for i, nuclide_elem in enumerate(root.findall('nuclide_table')):
nuc = Nuclide.from_xml(nuclide_elem)
chain.nuclide_dict[nuc.name] = i
# Check for reaction paths
for rx in nuc.reactions:
if rx.type not in chain.index_reaction:
chain.index_reaction[rx.type] = reaction_index
reaction_index += 1
if rx.type not in chain.reactions:
chain.reactions.append(rx.type)
chain.nuclides.append(nuc)

12
openmc/deplete/openmc_wrapper.py
View file

@ -366,10 +366,11 @@ class OpenMCOperator(Operator):
def initialize_reaction_rates(self):
"""Create reaction rates object. """
# Create dictionary to map reactions to indices
index_rx = {rx: i for i, rx in enumerate(self.chain.reactions)}
self.reaction_rates = ReactionRates(
self.burn_mat_to_ind,
self.burn_nuc_to_ind,
self.chain.index_reaction)
self.burn_mat_to_ind, self.burn_nuc_to_ind, index_rx)
def form_matrix(self, y, mat):
"""Forms the depletion matrix.
@ -520,7 +521,7 @@ class OpenMCOperator(Operator):
# transmutation. The nuclides for the tally are set later when eval() is
# called.
tally_dep = openmc.capi.Tally(1)
tally_dep.scores = self.chain.index_reaction.keys()
tally_dep.scores = self.chain.reactions
tally_dep.filters = [mat_filter]
def total_density_list(self):
@ -577,11 +578,10 @@ class OpenMCOperator(Operator):
# Extract tally bins
materials = list(self.mat_tally_ind.keys())
nuclides = openmc.capi.tallies[1].nuclides
reactions = list(self.chain.index_reaction.keys())
# Form fast map
nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides]
react_ind = [rates.index_rx[react] for react in reactions]
react_ind = [rates.index_rx[react] for react in self.chain.reactions]
# Compute fission power
# TODO : improve this calculation

59
openmc/deplete/reaction_rates.py
View file

@ -13,20 +13,20 @@ class ReactionRates(np.ndarray):
Parameters
----------
mat_to_ind : OrderedDict of str to int
index_mat : OrderedDict of str to int
A dictionary mapping material ID as string to index.
nuc_to_ind : OrderedDict of str to int
index_nuc : OrderedDict of str to int
A dictionary mapping nuclide name as string to index.
react_to_ind : OrderedDict of str to int
index_rx : OrderedDict of str to int
A dictionary mapping reaction name as string to index.
Attributes
----------
mat_to_ind : OrderedDict of str to int
index_mat : OrderedDict of str to int
A dictionary mapping material ID as string to index.
nuc_to_ind : OrderedDict of str to int
index_nuc : OrderedDict of str to int
A dictionary mapping nuclide name as string to index.
react_to_ind : OrderedDict of str to int
index_rx : OrderedDict of str to int
A dictionary mapping reaction name as string to index.
n_mat : int
Number of materials.
@ -34,10 +34,8 @@ class ReactionRates(np.ndarray):
Number of nucs.
n_react : int
Number of reactions.
rates : numpy.array
Array storing rates indexed by the above dictionaries.
"""
"""
def __new__(cls, index_mat, index_nuc, index_rx):
# Create appropriately-sized zeroed-out ndarray
shape = (len(index_mat), len(index_nuc), len(index_rx))
@ -83,3 +81,46 @@ class ReactionRates(np.ndarray):
def n_react(self):
"""Number of reactions."""
return len(self.index_rx)
def get(self, mat, nuc, rx):
"""Get reaction rate by material/nuclide/reaction
Parameters
----------
mat : str
Material ID as a string
nuc : str
Nuclide name
rx : str
Name of the reaction
Returns
-------
float
Reaction rate corresponding to given material, nuclide, and reaction
"""
mat = self.index_mat[mat]
nuc = self.index_nuc[nuc]
rx = self.index_rx[rx]
return self[mat, nuc, rx]
def set(self, mat, nuc, rx, value):
"""Set reaction rate by material/nuclide/reaction
Parameters
----------
mat : str
Material ID as a string
nuc : str
Nuclide name
rx : str
Name of the reaction
value : float
Corresponding reaction rate to set
"""
mat = self.index_mat[mat]
nuc = self.index_nuc[nuc]
rx = self.index_rx[rx]
self[mat, nuc, rx] = value

48
openmc/deplete/utilities.py
View file

@ -7,26 +7,26 @@ the results module.
import numpy as np
def evaluate_single_nuclide(results, cell, nuc):
"""Evaluates a single nuclide in a single cell from a results list.
def evaluate_single_nuclide(results, mat, nuc):
"""Evaluates a single nuclide in a single material from a results list.
Parameters
----------
results : list of results
The results to extract data from. Must be sorted and continuous.
cell : str
Cell name to evaluate
mat : str
Material name to evaluate
nuc : str
Nuclide name to evaluate
Returns
-------
time : numpy.array
Time vector.
concentration : numpy.array
Total number of atoms in the cell.
"""
time : numpy.ndarray
Time vector
concentration : numpy.ndarray
Total number of atoms in the material
"""
n_points = len(results)
time = np.zeros(n_points)
concentration = np.zeros(n_points)
@ -34,39 +34,39 @@ def evaluate_single_nuclide(results, cell, nuc):
# Evaluate value in each region
for i, result in enumerate(results):
time[i] = result.time[0]
concentration[i] = result[0, cell, nuc]
concentration[i] = result[0, mat, nuc]
return time, concentration
def evaluate_reaction_rate(results, cell, nuc, rxn):
"""Evaluates a single nuclide reaction rate in a single cell from a results list.
def evaluate_reaction_rate(results, mat, nuc, rx):
"""Return reaction rate in a single material/nuclide from a results list.
Parameters
----------
results : list of Results
results : list of openmc.deplete.Results
The results to extract data from. Must be sorted and continuous.
cell : str
Cell name to evaluate
mat : str
Material name to evaluate
nuc : str
Nuclide name to evaluate
rxn : str
rx : str
Reaction rate to evaluate
Returns
-------
time : numpy.array
time : numpy.ndarray
Time vector.
rate : numpy.array
rate : numpy.ndarray
Reaction rate.
"""
"""
n_points = len(results)
time = np.zeros(n_points)
rate = np.zeros(n_points)
# Evaluate value in each region
for i, result in enumerate(results):
time[i] = result.time[0]
rate[i] = result.rates[0][cell, nuc, rxn] * result[0, cell, nuc]
rate[i] = result.rates[0].get(mat, nuc, rx) * result[0, mat, nuc]
return time, rate
@ -76,17 +76,17 @@ def evaluate_eigenvalue(results):
Parameters
----------
results : list of Results
results : list of openmc.deplete.Results
The results to extract data from. Must be sorted and continuous.
Returns
-------
time : numpy.array
time : numpy.ndarray
Time vector.
eigenvalue : numpy.array
eigenvalue : numpy.ndarray
Eigenvalue.
"""
"""
n_points = len(results)
time = np.zeros(n_points)
eigenvalue = np.zeros(n_points)

57
tests/unit_tests/test_deplete_chain.py
View file

@ -13,19 +13,18 @@ _test_filename = str(Path(__file__).parents[2] / 'chains' / 'chain_test.xml')
def test_init():
"""Test depletion chain initialization."""
dep = Chain()
chain = Chain()
assert isinstance(dep.nuclides, list)
assert isinstance(dep.nuclide_dict, Mapping)
assert isinstance(dep.react_to_ind, Mapping)
assert isinstance(chain.nuclides, list)
assert isinstance(chain.nuclide_dict, Mapping)
def test_len():
"""Test depletion chain length."""
dep = Chain()
dep.nuclides = ["NucA", "NucB", "NucC"]
chain = Chain()
chain.nuclides = ["NucA", "NucB", "NucC"]
assert len(dep) == 3
assert len(chain) == 3
def test_from_endf():
@ -40,13 +39,13 @@ def test_from_xml():
# the components external to depletion_chain.py are simple storage
# types.
dep = Chain.from_xml(_test_filename)
chain = Chain.from_xml(_test_filename)
# Basic checks
assert len(dep) == 3
assert len(chain) == 3
# A tests
nuc = dep["A"]
nuc = chain["A"]
assert nuc.name == "A"
assert nuc.half_life == 2.36520E+04
@ -61,7 +60,7 @@ def test_from_xml():
assert [r.branching_ratio for r in nuc.reactions] == [1.0]
# B tests
nuc = dep["B"]
nuc = chain["B"]
assert nuc.name == "B"
assert nuc.half_life == 3.29040E+04
@ -76,7 +75,7 @@ def test_from_xml():
assert [r.branching_ratio for r in nuc.reactions] == [1.0]
# C tests
nuc = dep["C"]
nuc = chain["C"]
assert nuc.name == "C"
assert nuc.n_decay_modes == 0
@ -135,22 +134,20 @@ def test_form_matrix():
""" Using chain_test, and a dummy reaction rate, compute the matrix. """
# Relies on test_from_xml passing.
dep = Chain.from_xml(_test_filename)
chain = Chain.from_xml(_test_filename)
cell_ind = {"10000": 0, "10001": 1}
mat_ind = {"10000": 0, "10001": 1}
nuc_ind = {"A": 0, "B": 1, "C": 2}
react_ind = dep.react_to_ind
react_ind = {rx: i for i, rx in enumerate(chain.reactions)}
react = reaction_rates.ReactionRates(cell_ind, nuc_ind, react_ind)
react = reaction_rates.ReactionRates(mat_ind, nuc_ind, react_ind)
dep.nuc_to_react_ind = nuc_ind
react.set("10000", "C", "fission", 1.0)
react.set("10000", "A", "(n,gamma)", 2.0)
react.set("10000", "B", "(n,gamma)", 3.0)
react.set("10000", "C", "(n,gamma)", 4.0)
react["10000", "C", "fission"] = 1.0
react["10000", "A", "(n,gamma)"] = 2.0
react["10000", "B", "(n,gamma)"] = 3.0
react["10000", "C", "(n,gamma)"] = 4.0
mat = dep.form_matrix(react[0, :, :])
mat = chain.form_matrix(react[0, :, :])
# Loss A, decay, (n, gamma)
mat00 = -np.log(2) / 2.36520E+04 - 2
# A -> B, decay, 0.6 branching ratio
@ -185,11 +182,11 @@ def test_form_matrix():
def test_getitem():
""" Test nuc_by_ind converter function. """
dep = Chain()
dep.nuclides = ["NucA", "NucB", "NucC"]
dep.nuclide_dict = {nuc: dep.nuclides.index(nuc)
for nuc in dep.nuclides}
chain = Chain()
chain.nuclides = ["NucA", "NucB", "NucC"]
chain.nuclide_dict = {nuc: chain.nuclides.index(nuc)
for nuc in chain.nuclides}
assert "NucA" == dep["NucA"]
assert "NucB" == dep["NucB"]
assert "NucC" == dep["NucC"]
assert "NucA" == chain["NucA"]
assert "NucB" == chain["NucB"]
assert "NucC" == chain["NucC"]

6
tests/unit_tests/test_deplete_integrator.py
View file

@ -86,10 +86,8 @@ def test_save_results(run_in_tmpdir):
for nuc_i, nuc in enumerate(nuc_list):
assert res[0][i, mat, nuc] == x1[i][mat_i][nuc_i]
assert res[1][i, mat, nuc] == x2[i][mat_i][nuc_i]
np.testing.assert_array_equal(res[0].rates[i][mat, nuc, :],
rate1[i][mat, nuc, :])
np.testing.assert_array_equal(res[1].rates[i][mat, nuc, :],
rate2[i][mat, nuc, :])
np.testing.assert_array_equal(res[0].rates[i], rate1[i])
np.testing.assert_array_equal(res[1].rates[i], rate2[i])
np.testing.assert_array_equal(res[0].k, eigvl1)
np.testing.assert_array_equal(res[0].time, t1)

51
tests/unit_tests/test_deplete_reaction.py
View file

@ -1,35 +1,38 @@
"""Tests for the openmc.deplete.ReactionRates class."""
from openmc.deplete import reaction_rates
import numpy as np
from openmc.deplete import ReactionRates
def test_indexing():
"""Tests the __getitem__ and __setitem__ routines simultaneously."""
def test_get_set():
"""Tests the get/set methods."""
mat_to_ind = {"10000" : 0, "10001" : 1}
nuc_to_ind = {"U238" : 0, "U235" : 1}
react_to_ind = {"fission" : 0, "(n,gamma)" : 1}
rates = reaction_rates.ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
rates = ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
assert rates.shape == (2, 2, 2)
assert np.all(rates == 0.0)
rates["10000", "U238", "fission"] = 1.0
rates["10001", "U238", "fission"] = 2.0
rates["10000", "U235", "fission"] = 3.0
rates["10001", "U235", "fission"] = 4.0
rates["10000", "U238", "(n,gamma)"] = 5.0
rates["10001", "U238", "(n,gamma)"] = 6.0
rates["10000", "U235", "(n,gamma)"] = 7.0
rates["10001", "U235", "(n,gamma)"] = 8.0
rates.set("10000", "U238", "fission", 1.0)
rates.set("10001", "U238", "fission", 2.0)
rates.set("10000", "U235", "fission", 3.0)
rates.set("10001", "U235", "fission", 4.0)
rates.set("10000", "U238", "(n,gamma)", 5.0)
rates.set("10001", "U238", "(n,gamma)", 6.0)
rates.set("10000", "U235", "(n,gamma)", 7.0)
rates.set("10001", "U235", "(n,gamma)", 8.0)
# String indexing
assert rates["10000", "U238", "fission"] == 1.0
assert rates["10001", "U238", "fission"] == 2.0
assert rates["10000", "U235", "fission"] == 3.0
assert rates["10001", "U235", "fission"] == 4.0
assert rates["10000", "U238", "(n,gamma)"] == 5.0
assert rates["10001", "U238", "(n,gamma)"] == 6.0
assert rates["10000", "U235", "(n,gamma)"] == 7.0
assert rates["10001", "U235", "(n,gamma)"] == 8.0
assert rates.get("10000", "U238", "fission") == 1.0
assert rates.get("10001", "U238", "fission") == 2.0
assert rates.get("10000", "U235", "fission") == 3.0
assert rates.get("10001", "U235", "fission") == 4.0
assert rates.get("10000", "U238", "(n,gamma)") == 5.0
assert rates.get("10001", "U238", "(n,gamma)") == 6.0
assert rates.get("10000", "U235", "(n,gamma)") == 7.0
assert rates.get("10001", "U235", "(n,gamma)") == 8.0
# Int indexing
assert rates[0, 0, 0] == 1.0
@ -44,7 +47,7 @@ def test_indexing():
rates[0, 0, 0] = 5.0
assert rates[0, 0, 0] == 5.0
assert rates["10000", "U238", "fission"] == 5.0
assert rates.get("10000", "U238", "fission") == 5.0
def test_n_mat():
@ -53,7 +56,7 @@ def test_n_mat():
nuc_to_ind = {"U238" : 0, "U235" : 1, "Gd157" : 2}
react_to_ind = {"fission" : 0, "(n,gamma)" : 1, "(n,2n)" : 2, "(n,3n)" : 3}
rates = reaction_rates.ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
rates = ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
assert rates.n_mat == 2
@ -64,7 +67,7 @@ def test_n_nuc():
nuc_to_ind = {"U238" : 0, "U235" : 1, "Gd157" : 2}
react_to_ind = {"fission" : 0, "(n,gamma)" : 1, "(n,2n)" : 2, "(n,3n)" : 3}
rates = reaction_rates.ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
rates = ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
assert rates.n_nuc == 3
@ -75,6 +78,6 @@ def test_n_react():
nuc_to_ind = {"U238" : 0, "U235" : 1, "Gd157" : 2}
react_to_ind = {"fission" : 0, "(n,gamma)" : 1, "(n,2n)" : 2, "(n,3n)" : 3}
rates = reaction_rates.ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
rates = ReactionRates(mat_to_ind, nuc_to_ind, react_to_ind)
assert rates.n_react == 4