From 476515e9fde83a22f117b364f1422abeec51a2c1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 30 Apr 2020 06:34:19 -0500 Subject: [PATCH] Make sure (n,p), (n,d), etc. reactions produce light nuclides in depletion --- openmc/deplete/chain.py | 10 +++ tests/unit_tests/test_deplete_chain.py | 96 ++++++++++++++------------ 2 files changed, 61 insertions(+), 45 deletions(-) diff --git a/openmc/deplete/chain.py b/openmc/deplete/chain.py index b78590bc25..bc42b13ecb 100644 --- a/openmc/deplete/chain.py +++ b/openmc/deplete/chain.py @@ -531,6 +531,16 @@ class Chain: if target is not None and path_rate != 0.0: k = self.nuclide_dict[target] matrix[k, i] += path_rate * br + + # Determine light nuclide production, e.g., (n,d) should + # produce H2 + light_nucs = _SECONDARY_PARTICLES.get(r_type) + if light_nucs is not None: + for light_nuc in light_nucs: + k = self.nuclide_dict.get(light_nuc) + if k is not None: + matrix[k, i] += path_rate * br + else: for product, y in fission_yields[nuc.name].items(): yield_val = y * path_rate diff --git a/tests/unit_tests/test_deplete_chain.py b/tests/unit_tests/test_deplete_chain.py index ba3be4e37a..97053a4911 100644 --- a/tests/unit_tests/test_deplete_chain.py +++ b/tests/unit_tests/test_deplete_chain.py @@ -1,9 +1,10 @@ """Tests for openmc.deplete.Chain class.""" from collections.abc import Mapping +from itertools import product +from math import log import os from pathlib import Path -from itertools import product import numpy as np from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram @@ -13,14 +14,17 @@ from tests import cdtemp _TEST_CHAIN = """\ - + + + - + + @@ -84,7 +88,13 @@ def test_from_xml(simple_chain): chain = simple_chain # Basic checks - assert len(chain) == 3 + assert len(chain) == 4 + + # H1 tests + nuc = chain["H1"] + assert nuc.name == "H1" + assert nuc.n_decay_modes == 0 + assert nuc.n_reaction_paths == 0 # A tests nuc = chain["A"] @@ -96,10 +106,10 @@ def test_from_xml(simple_chain): assert [m.target for m in modes] == ["B", "C"] assert [m.type for m in modes] == ["beta1", "beta2"] assert [m.branching_ratio for m in modes] == [0.6, 0.4] - assert nuc.n_reaction_paths == 1 - assert [r.target for r in nuc.reactions] == ["C"] - assert [r.type for r in nuc.reactions] == ["(n,gamma)"] - assert [r.branching_ratio for r in nuc.reactions] == [1.0] + assert nuc.n_reaction_paths == 2 + assert [r.target for r in nuc.reactions] == ["C", "B"] + assert [r.type for r in nuc.reactions] == ["(n,gamma)", "(n,p)"] + assert [r.branching_ratio for r in nuc.reactions] == [1.0, 1.0] # B tests nuc = chain["B"] @@ -111,10 +121,10 @@ def test_from_xml(simple_chain): assert [m.target for m in modes] == ["A"] assert [m.type for m in modes] == ["beta"] assert [m.branching_ratio for m in modes] == [1.0] - assert nuc.n_reaction_paths == 1 - assert [r.target for r in nuc.reactions] == ["C"] - assert [r.type for r in nuc.reactions] == ["(n,gamma)"] - assert [r.branching_ratio for r in nuc.reactions] == [1.0] + assert nuc.n_reaction_paths == 2 + assert [r.target for r in nuc.reactions] == ["C", "A"] + assert [r.type for r in nuc.reactions] == ["(n,gamma)", "(n,d)"] + assert [r.branching_ratio for r in nuc.reactions] == [1.0, 1.0] # C tests nuc = chain["C"] @@ -139,18 +149,26 @@ def test_export_to_xml(run_in_tmpdir): # Prevent different MPI ranks from conflicting filename = 'test{}.xml'.format(comm.rank) + H1 = nuclide.Nuclide("H1") + A = nuclide.Nuclide("A") A.half_life = 2.36520e4 A.decay_modes = [ nuclide.DecayTuple("beta1", "B", 0.6), nuclide.DecayTuple("beta2", "C", 0.4) ] - A.reactions = [nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0)] + A.reactions = [ + nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0), + nuclide.ReactionTuple("(n,p)", "B", 0.0, 1.0) + ] B = nuclide.Nuclide("B") B.half_life = 3.29040e4 B.decay_modes = [nuclide.DecayTuple("beta", "A", 1.0)] - B.reactions = [nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0)] + B.reactions = [ + nuclide.ReactionTuple("(n,gamma)", "C", 0.0, 1.0), + nuclide.ReactionTuple("(n,d)", "A", 0.0, 1.0) + ] C = nuclide.Nuclide("C") C.reactions = [ @@ -162,7 +180,7 @@ def test_export_to_xml(run_in_tmpdir): 0.0253: {"A": 0.0292737, "B": 0.002566345}}) chain = Chain() - chain.nuclides = [A, B, C] + chain.nuclides = [H1, A, B, C] chain.export_to_xml(filename) chain_xml = open(filename, 'r').read() @@ -179,43 +197,31 @@ def test_form_matrix(simple_chain): nuclides = ["A", "B", "C"] react = reaction_rates.ReactionRates(mats, nuclides, chain.reactions) - react.set("10000", "C", "fission", 1.0) react.set("10000", "A", "(n,gamma)", 2.0) + react.set("10000", "A", "(n,p)", 0.1) react.set("10000", "B", "(n,gamma)", 3.0) + react.set("10000", "B", "(n,d)", 0.2) react.set("10000", "C", "(n,gamma)", 4.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 - mat10 = np.log(2) / 2.36520E+04 * 0.6 - # A -> C, decay, 0.4 branching ratio + (n,gamma) - mat20 = np.log(2) / 2.36520E+04 * 0.4 + 2 + decay_constant = log(2) / 2.36520E+04 + mat = np.zeros((4, 4)) + mat[0, 1] = 0.1 # A -> B, (n,p) + mat[1, 1] = -decay_constant - 2.1 # Loss A, decay, (n,gamma), (n,p) + mat[2, 1] = decay_constant*0.6 + 0.1 # A -> B, decay, 0.6 branching ratio + (n,p) + mat[3, 1] = decay_constant*0.4 + 2 # A -> C, decay, 0.4 branching ratio + (n,gamma) - # B -> A, decay, 1.0 branching ratio - mat01 = np.log(2)/3.29040E+04 - # Loss B, decay, (n, gamma) - mat11 = -np.log(2)/3.29040E+04 - 3 - # B -> C, (n, gamma) - mat21 = 3 + decay_constant = log(2.0) / 3.29040E+04 + mat[1, 2] = decay_constant + 0.2 # B -> A, decay, 1.0 branching ratio + (n,d) + mat[2, 2] = -decay_constant - 3.2 # Loss B, decay, (n,gamma), (n,d) + mat[3, 2] = 3 # B -> C, (n,gamma) - # C -> A fission, (n, gamma) - mat02 = 0.0292737 * 1.0 + 4.0 * 0.7 - # C -> B fission, (n, gamma) - mat12 = 0.002566345 * 1.0 + 4.0 * 0.3 - # Loss C, fission, (n, gamma) - mat22 = -1.0 - 4.0 + mat[1, 3] = 0.0292737 * 1.0 + 4.0 * 0.7 # C -> A fission, (n,gamma) + mat[2, 3] = 0.002566345 * 1.0 + 4.0 * 0.3 # C -> B fission, (n,gamma) + mat[3, 3] = -1.0 - 4.0 # Loss C, fission, (n,gamma) - assert mat[0, 0] == mat00 - assert mat[1, 0] == mat10 - assert mat[2, 0] == mat20 - assert mat[0, 1] == mat01 - assert mat[1, 1] == mat11 - assert mat[2, 1] == mat21 - assert mat[0, 2] == mat02 - assert mat[1, 2] == mat12 - assert mat[2, 2] == mat22 + sp_matrix = chain.form_matrix(react[0]) + assert np.allclose(mat, sp_matrix.toarray()) # Pass equivalent fission yields directly # Ensure identical matrix is formed @@ -226,7 +232,7 @@ def test_form_matrix(simple_chain): def test_getitem(): - """ Test nuc_by_ind converter function. """ + """Test nuc_by_ind converter function.""" chain = Chain() chain.nuclides = ["NucA", "NucB", "NucC"] chain.nuclide_dict = {nuc: chain.nuclides.index(nuc)