diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index 2187b04860..2af64e172a 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -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) diff --git a/openmc/deplete/openmc_wrapper.py b/openmc/deplete/openmc_wrapper.py index 46d07a23f7..c95170a461 100644 --- a/openmc/deplete/openmc_wrapper.py +++ b/openmc/deplete/openmc_wrapper.py @@ -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 diff --git a/openmc/deplete/reaction_rates.py b/openmc/deplete/reaction_rates.py index b437773826..a479085174 100644 --- a/openmc/deplete/reaction_rates.py +++ b/openmc/deplete/reaction_rates.py @@ -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 diff --git a/openmc/deplete/utilities.py b/openmc/deplete/utilities.py index 5433edce44..d155f321d9 100644 --- a/openmc/deplete/utilities.py +++ b/openmc/deplete/utilities.py @@ -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) diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index 3a065fb46f..ae900e62d7 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -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"] diff --git a/tests/unit_tests/test_deplete_integrator.py b/tests/unit_tests/test_deplete_integrator.py index 59d08b4842..78dd6bcef2 100644 --- a/tests/unit_tests/test_deplete_integrator.py +++ b/tests/unit_tests/test_deplete_integrator.py @@ -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) diff --git a/tests/unit_tests/test_deplete_reaction.py b/tests/unit_tests/test_deplete_reaction.py index a98535e1d8..de628f8c66 100644 --- a/tests/unit_tests/test_deplete_reaction.py +++ b/tests/unit_tests/test_deplete_reaction.py @@ -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