mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Merge pull request #1528 from drewejohnson/condense-chain
Use Chain.reduce to follow paths and reduce the size
This commit is contained in:
Paul Romano
GitHub
commit
72a639593d
5 changed files with 354 additions and 30 deletions
|
|
@ -10,7 +10,7 @@ import math
|
|||
import re
|
||||
from collections import OrderedDict, defaultdict
|
||||
from collections.abc import Mapping, Iterable
|
||||
from numbers import Real
|
||||
from numbers import Real, Integral
|
||||
from warnings import warn
|
||||
|
||||
from openmc.checkvalue import check_type, check_greater_than
|
||||
|
|
@ -500,7 +500,7 @@ class Chain(object):
|
|||
# Gain
|
||||
for _, target, branching_ratio in nuc.decay_modes:
|
||||
# Allow for total annihilation for debug purposes
|
||||
if target != 'Nothing':
|
||||
if target is not None:
|
||||
branch_val = branching_ratio * decay_constant
|
||||
|
||||
if branch_val != 0.0:
|
||||
|
|
@ -525,17 +525,16 @@ class Chain(object):
|
|||
matrix[i, i] -= path_rate
|
||||
|
||||
# Gain term; allow for total annihilation for debug purposes
|
||||
if target != 'Nothing':
|
||||
if r_type != 'fission':
|
||||
if path_rate != 0.0:
|
||||
k = self.nuclide_dict[target]
|
||||
matrix[k, i] += path_rate * br
|
||||
else:
|
||||
for product, y in fission_yields[nuc.name].items():
|
||||
yield_val = y * path_rate
|
||||
if yield_val != 0.0:
|
||||
k = self.nuclide_dict[product]
|
||||
matrix[k, i] += yield_val
|
||||
if r_type != 'fission':
|
||||
if target is not None and path_rate != 0.0:
|
||||
k = self.nuclide_dict[target]
|
||||
matrix[k, i] += path_rate * br
|
||||
else:
|
||||
for product, y in fission_yields[nuc.name].items():
|
||||
yield_val = y * path_rate
|
||||
if yield_val != 0.0:
|
||||
k = self.nuclide_dict[product]
|
||||
matrix[k, i] += yield_val
|
||||
|
||||
# Clear set of reactions
|
||||
reactions.clear()
|
||||
|
|
@ -821,3 +820,168 @@ class Chain(object):
|
|||
return stat
|
||||
valid = valid and stat
|
||||
return valid
|
||||
|
||||
def reduce(self, initial_isotopes, level=None):
|
||||
"""Reduce the size of the chain by following transmutation paths
|
||||
|
||||
As an example, consider a simple chain with the following
|
||||
isotopes and transmutation paths::
|
||||
|
||||
U235 (n,gamma) U236
|
||||
(n,fission) (Xe135, I135, Cs135)
|
||||
I135 (beta decay) Xe135 (beta decay) Cs135
|
||||
Xe135 (n,gamma) Xe136
|
||||
|
||||
Calling ``chain.reduce(["I135"])`` will produce a depletion
|
||||
chain that contains only isotopes that would originate from
|
||||
I135: I135, Xe135, Cs135, and Xe136. U235 and U236 will not
|
||||
be included, but multiple isotopes can be used to start
|
||||
the search.
|
||||
|
||||
The ``level`` value controls the depth of the search.
|
||||
``chain.reduce(["U235"], level=1)`` would return a chain
|
||||
with all isotopes except Xe136, since it is two transmutations
|
||||
removed from U235 in this case.
|
||||
|
||||
While targets will not be included in the new chain, the
|
||||
total destruction rate and decay rate of included isotopes
|
||||
will be preserved.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
initial_isotopes : iterable of str
|
||||
Start the search based on the contents of these isotopes
|
||||
level : int, optional
|
||||
Depth of transmuation path to follow. Must be greater than
|
||||
or equal to zero. A value of zero returns a chain with
|
||||
``initial_isotopes``. The default value of None implies
|
||||
that all isotopes that appear in the transmutation paths
|
||||
of the initial isotopes and their progeny should be
|
||||
explored
|
||||
|
||||
Returns
|
||||
-------
|
||||
Chain
|
||||
Depletion chain containing isotopes that would appear
|
||||
after following up to ``level`` reactions and decay paths
|
||||
|
||||
"""
|
||||
check_type("initial_isotopes", initial_isotopes, Iterable, str)
|
||||
if level is None:
|
||||
level = math.inf
|
||||
else:
|
||||
check_type("level", level, Integral)
|
||||
check_greater_than("level", level, 0, equality=True)
|
||||
|
||||
all_isotopes = self._follow(set(initial_isotopes), level)
|
||||
|
||||
# Avoid re-sorting for fission yields
|
||||
name_sort = sorted(all_isotopes)
|
||||
|
||||
nuclides = []
|
||||
nuclide_dict = {}
|
||||
reactions = set()
|
||||
|
||||
for idx, iso in enumerate(sorted(all_isotopes, key=openmc.data.zam)):
|
||||
previous = self[iso]
|
||||
new_nuclide = Nuclide(previous.name)
|
||||
new_nuclide.half_life = previous.half_life
|
||||
new_nuclide.decay_energy = new_nuclide.decay_energy
|
||||
|
||||
new_decay = []
|
||||
for mode in previous.decay_modes:
|
||||
if mode.target in all_isotopes:
|
||||
new_decay.append(mode)
|
||||
else:
|
||||
new_decay.append(DecayTuple(
|
||||
mode.type, None, mode.branching_ratio))
|
||||
new_nuclide.decay_modes = new_decay
|
||||
|
||||
new_reactions = []
|
||||
for rxn in previous.reactions:
|
||||
if rxn.target in all_isotopes:
|
||||
new_reactions.append(rxn)
|
||||
reactions.add(rxn.type)
|
||||
elif rxn.type == "fission":
|
||||
new_yields = new_nuclide.yield_data = (
|
||||
previous.yield_data.restrict_products(name_sort))
|
||||
if new_yields is not None:
|
||||
new_reactions.append(rxn)
|
||||
reactions.add("fission")
|
||||
# Maintain total destruction rates but set no target
|
||||
else:
|
||||
new_reactions.append(ReactionTuple(
|
||||
rxn.type, None, rxn.Q, rxn.branching_ratio))
|
||||
reactions.add(rxn.type)
|
||||
|
||||
new_nuclide.reactions = new_reactions
|
||||
|
||||
nuclides.append(new_nuclide)
|
||||
nuclide_dict[iso] = idx
|
||||
|
||||
new_chain = type(self)()
|
||||
new_chain.nuclides = nuclides
|
||||
new_chain.nuclide_dict = nuclide_dict
|
||||
|
||||
# Doesn't appear that the ordering matters for the reactions,
|
||||
# just the contents
|
||||
new_chain.reactions = sorted(reactions)
|
||||
|
||||
return new_chain
|
||||
|
||||
def _follow(self, isotopes, level):
|
||||
"""Return all isotopes present up to depth level"""
|
||||
found = isotopes.copy()
|
||||
remaining = set(self.nuclide_dict)
|
||||
if not found.issubset(remaining):
|
||||
raise IndexError(
|
||||
"The following isotopes were not found in the chain: "
|
||||
"{}".format(", ".join(found - remaining)))
|
||||
|
||||
if level == 0:
|
||||
return found
|
||||
|
||||
remaining -= found
|
||||
|
||||
depth = 0
|
||||
next_iso = set()
|
||||
|
||||
while depth < level and remaining:
|
||||
# Exhaust all isotopes at this level
|
||||
while isotopes:
|
||||
iso = isotopes.pop()
|
||||
found.add(iso)
|
||||
nuclide = self[iso]
|
||||
|
||||
# Follow all transmutation paths for this nuclide
|
||||
for rxn in nuclide.reactions + nuclide.decay_modes:
|
||||
if rxn.type == "fission" or rxn.target is None:
|
||||
continue
|
||||
# Skip if we've already come across this isotope
|
||||
elif (rxn.target in next_iso
|
||||
or rxn.target in found or rxn.target in isotopes):
|
||||
continue
|
||||
next_iso.add(rxn.target)
|
||||
|
||||
if nuclide.yield_data is not None:
|
||||
for product in nuclide.yield_data.products:
|
||||
if (product in next_iso
|
||||
or product in found or product in isotopes):
|
||||
continue
|
||||
next_iso.add(product)
|
||||
|
||||
if not next_iso:
|
||||
# No additional isotopes to process, nor to update the
|
||||
# current set of discovered isotopes
|
||||
return found
|
||||
|
||||
# Prepare for next dig
|
||||
depth += 1
|
||||
isotopes |= next_iso
|
||||
remaining -= next_iso
|
||||
next_iso.clear()
|
||||
|
||||
# Process isotope that would have started next depth
|
||||
found.update(isotopes)
|
||||
|
||||
return found
|
||||
|
|
|
|||
|
|
@ -13,7 +13,7 @@ try:
|
|||
except ImportError:
|
||||
import xml.etree.ElementTree as ET
|
||||
|
||||
from numpy import empty
|
||||
from numpy import empty, searchsorted
|
||||
|
||||
from openmc.checkvalue import check_type
|
||||
|
||||
|
|
@ -30,8 +30,10 @@ Parameters
|
|||
----------
|
||||
type : str
|
||||
Type of the decay mode, e.g., 'beta-'
|
||||
target : str
|
||||
Nuclide resulting from decay
|
||||
target : str or None
|
||||
Nuclide resulting from decay. A value of ``None`` implies the
|
||||
target does not exist in the currently configured depletion
|
||||
chain
|
||||
branching_ratio : float
|
||||
Branching ratio of the decay mode
|
||||
|
||||
|
|
@ -53,8 +55,11 @@ Parameters
|
|||
----------
|
||||
type : str
|
||||
Type of the reaction, e.g., 'fission'
|
||||
target : str
|
||||
nuclide resulting from reaction
|
||||
target : str or None
|
||||
Nuclide resulting from reaction. A value of ``None``
|
||||
implies either no single target, e.g. from fission,
|
||||
or that the target nuclide is not considered
|
||||
in the current depletion chain
|
||||
Q : float
|
||||
Q value of the reaction in [eV]
|
||||
branching_ratio : float
|
||||
|
|
@ -179,6 +184,8 @@ class Nuclide(object):
|
|||
for decay_elem in element.iter('decay'):
|
||||
d_type = decay_elem.get('type')
|
||||
target = decay_elem.get('target')
|
||||
if target is not None and target.lower() == "nothing":
|
||||
target = None
|
||||
branching_ratio = float(decay_elem.get('branching_ratio'))
|
||||
nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio))
|
||||
|
||||
|
|
@ -192,6 +199,8 @@ class Nuclide(object):
|
|||
# just set null values
|
||||
if r_type != 'fission':
|
||||
target = reaction_elem.get('target')
|
||||
if target is not None and target.lower() == "nothing":
|
||||
target = None
|
||||
else:
|
||||
target = None
|
||||
if fission_q is not None:
|
||||
|
|
@ -226,7 +235,7 @@ class Nuclide(object):
|
|||
for mode, daughter, br in self.decay_modes:
|
||||
mode_elem = ET.SubElement(elem, 'decay')
|
||||
mode_elem.set('type', mode)
|
||||
mode_elem.set('target', daughter)
|
||||
mode_elem.set('target', daughter or "Nothing")
|
||||
mode_elem.set('branching_ratio', str(br))
|
||||
|
||||
elem.set('reactions', str(len(self.reactions)))
|
||||
|
|
@ -234,7 +243,7 @@ class Nuclide(object):
|
|||
rx_elem = ET.SubElement(elem, 'reaction')
|
||||
rx_elem.set('type', rx)
|
||||
rx_elem.set('Q', str(Q))
|
||||
if rx != 'fission':
|
||||
if rx != 'fission' or daughter is not None:
|
||||
rx_elem.set('target', daughter)
|
||||
if br != 1.0:
|
||||
rx_elem.set('branching_ratio', str(br))
|
||||
|
|
@ -397,9 +406,7 @@ class FissionYieldDistribution(Mapping):
|
|||
for g_index, energy in enumerate(energies):
|
||||
prod_map = fission_yields[energy]
|
||||
for prod_ix, product in enumerate(ordered_prod):
|
||||
yield_val = prod_map.get(product)
|
||||
yield_matrix[g_index, prod_ix] = (
|
||||
0.0 if yield_val is None else yield_val)
|
||||
yield_matrix[g_index, prod_ix] = prod_map.get(product, 0.0)
|
||||
self.energies = tuple(energies)
|
||||
self.products = tuple(ordered_prod)
|
||||
self.yield_matrix = yield_matrix
|
||||
|
|
@ -435,7 +442,7 @@ class FissionYieldDistribution(Mapping):
|
|||
FissionYieldDistribution
|
||||
"""
|
||||
all_yields = {}
|
||||
for elem_index, yield_elem in enumerate(element.iter("fission_yields")):
|
||||
for yield_elem in element.iter("fission_yields"):
|
||||
energy = float(yield_elem.get("energy"))
|
||||
products = yield_elem.find("products").text.split()
|
||||
yields = map(float, yield_elem.find("data").text.split())
|
||||
|
|
@ -460,6 +467,37 @@ class FissionYieldDistribution(Mapping):
|
|||
data_elem = ET.SubElement(yield_element, "data")
|
||||
data_elem.text = " ".join(map(str, yield_obj.yields))
|
||||
|
||||
def restrict_products(self, possible_products):
|
||||
"""Return a new distribution with select products
|
||||
|
||||
Parameters
|
||||
----------
|
||||
possible_products : iterable of str
|
||||
Candidate pool of fission products. Existing products
|
||||
not contained here will not exist in the new instance
|
||||
|
||||
Returns
|
||||
-------
|
||||
FissionYieldDistribution or None
|
||||
A value of None indicates no values in
|
||||
``possible_products`` exist in :attr:`products`
|
||||
|
||||
"""
|
||||
|
||||
overlap = set(self.products).intersection(possible_products)
|
||||
if not overlap:
|
||||
return None
|
||||
|
||||
products = sorted(overlap)
|
||||
indices = searchsorted(self.products, products)
|
||||
|
||||
# coerce back to dictionary to pass back to __init__
|
||||
new_yields = {}
|
||||
for ene, yields in zip(self.energies, self.yield_matrix.copy()):
|
||||
new_yields[ene] = dict(zip(products, yields[indices]))
|
||||
|
||||
return type(self)(new_yields)
|
||||
|
||||
|
||||
class FissionYield(Mapping):
|
||||
"""Mapping for fission yields of a parent at a specific energy
|
||||
|
|
|
|||
|
|
@ -10,13 +10,10 @@ densities is all done in-memory instead of through the filesystem.
|
|||
import sys
|
||||
import copy
|
||||
from collections import OrderedDict
|
||||
from itertools import chain
|
||||
import os
|
||||
import time
|
||||
import xml.etree.ElementTree as ET
|
||||
from warnings import warn
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
from uncertainties import ufloat
|
||||
|
||||
|
|
@ -113,6 +110,13 @@ class Operator(TransportOperator):
|
|||
``fission_yield_mode``. Will be passed directly on to the
|
||||
helper. Passing a value of None will use the defaults for
|
||||
the associated helper.
|
||||
reduce_chain : bool, optional
|
||||
If True, use :meth:`openmc.deplete.Chain.reduce` to reduce the
|
||||
depletion chain up to ``reduce_chain_level``. Default is False.
|
||||
reduce_chain_level : int, optional
|
||||
Depth of the search when reducing the depletion chain. Only used
|
||||
if ``reduce_chain`` evaluates to true. The default value of
|
||||
``None`` implies no limit on the depth.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
|
|
@ -158,7 +162,8 @@ class Operator(TransportOperator):
|
|||
def __init__(self, geometry, settings, chain_file=None, prev_results=None,
|
||||
diff_burnable_mats=False, energy_mode="fission-q",
|
||||
fission_q=None, dilute_initial=1.0e3,
|
||||
fission_yield_mode="constant", fission_yield_opts=None):
|
||||
fission_yield_mode="constant", fission_yield_opts=None,
|
||||
reduce_chain=False, reduce_chain_level=None):
|
||||
if fission_yield_mode not in self._fission_helpers:
|
||||
raise KeyError(
|
||||
"fission_yield_mode must be one of {}, not {}".format(
|
||||
|
|
@ -179,6 +184,16 @@ class Operator(TransportOperator):
|
|||
self.geometry = geometry
|
||||
self.diff_burnable_mats = diff_burnable_mats
|
||||
|
||||
# Reduce the chain before we create more materials
|
||||
if reduce_chain:
|
||||
all_isotopes = set()
|
||||
for material in geometry.get_all_materials().values():
|
||||
if not material.depletable:
|
||||
continue
|
||||
for name, _dens_percent, _dens_type in material.nuclides:
|
||||
all_isotopes.add(name)
|
||||
self.chain = self.chain.reduce(all_isotopes, reduce_chain_level)
|
||||
|
||||
# Differentiate burnable materials with multiple instances
|
||||
if self.diff_burnable_mats:
|
||||
self._differentiate_burnable_mats()
|
||||
|
|
|
|||
|
|
@ -6,7 +6,6 @@ from pathlib import Path
|
|||
from itertools import product
|
||||
|
||||
import numpy as np
|
||||
from openmc.data import zam, ATOMIC_SYMBOL
|
||||
from openmc.deplete import comm, Chain, reaction_rates, nuclide, cram
|
||||
import pytest
|
||||
|
||||
|
|
@ -405,7 +404,8 @@ def test_fission_yield_attribute(simple_chain):
|
|||
dummy_conc = [[1, 2]] * (len(empty_chain.fission_yields) + 1)
|
||||
with pytest.raises(
|
||||
ValueError, match="fission yield.*not equal.*compositions"):
|
||||
cram.deplete(empty_chain, dummy_conc, None, 0.5)
|
||||
cram.deplete(empty_chain, dummy_conc, None, 0.5)
|
||||
|
||||
|
||||
def test_validate(simple_chain):
|
||||
"""Test the validate method"""
|
||||
|
|
@ -457,3 +457,95 @@ def test_validate_inputs():
|
|||
|
||||
with pytest.raises(ValueError, match="tolerance"):
|
||||
c.validate(tolerance=-1)
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def gnd_simple_chain():
|
||||
chainfile = Path(__file__).parents[1] / "chain_simple.xml"
|
||||
return Chain.from_xml(chainfile)
|
||||
|
||||
|
||||
def test_reduce(gnd_simple_chain):
|
||||
ref_U5 = gnd_simple_chain["U235"]
|
||||
ref_iodine = gnd_simple_chain["I135"]
|
||||
ref_U5_yields = ref_U5.yield_data
|
||||
|
||||
no_depth = gnd_simple_chain.reduce(["U235", "I135"], 0)
|
||||
# We should get a chain just containing U235 and I135
|
||||
assert len(no_depth) == 2
|
||||
assert set(no_depth.reactions) == set(gnd_simple_chain.reactions)
|
||||
|
||||
u5_round0 = no_depth["U235"]
|
||||
assert u5_round0.n_decay_modes == ref_U5.n_decay_modes
|
||||
for newmode, refmode in zip(u5_round0.decay_modes, ref_U5.decay_modes):
|
||||
assert newmode.target is None
|
||||
assert newmode.type == refmode.type
|
||||
assert newmode.branching_ratio == refmode.branching_ratio
|
||||
|
||||
assert u5_round0.n_reaction_paths == ref_U5.n_reaction_paths
|
||||
for newrxn, refrxn in zip(u5_round0.reactions, ref_U5.reactions):
|
||||
assert newrxn.target is None
|
||||
assert newrxn.type == refrxn.type
|
||||
assert newrxn.Q == refrxn.Q
|
||||
assert newrxn.branching_ratio == refrxn.branching_ratio
|
||||
|
||||
assert u5_round0.yield_data is not None
|
||||
assert u5_round0.yield_data.products == ("I135",)
|
||||
assert u5_round0.yield_data.yield_matrix == (
|
||||
ref_U5_yields.yield_matrix[:, ref_U5_yields.products.index("I135")]
|
||||
)
|
||||
|
||||
bareI5 = no_depth["I135"]
|
||||
assert bareI5.n_decay_modes == ref_iodine.n_decay_modes
|
||||
for newmode, refmode in zip(bareI5.decay_modes, ref_iodine.decay_modes):
|
||||
assert newmode.target is None
|
||||
assert newmode.type == refmode.type
|
||||
assert newmode.branching_ratio == refmode.branching_ratio
|
||||
|
||||
assert bareI5.n_reaction_paths == ref_iodine.n_reaction_paths
|
||||
for newrxn, refrxn in zip(bareI5.reactions, ref_iodine.reactions):
|
||||
assert newrxn.target is None
|
||||
assert newrxn.type == refrxn.type
|
||||
assert newrxn.Q == refrxn.Q
|
||||
assert newrxn.branching_ratio == refrxn.branching_ratio
|
||||
|
||||
follow_u5 = gnd_simple_chain.reduce(["U235"], 1)
|
||||
u5_round1 = follow_u5["U235"]
|
||||
assert u5_round1.decay_modes == ref_U5.decay_modes
|
||||
assert u5_round1.reactions == ref_U5.reactions
|
||||
assert u5_round1.yield_data is not None
|
||||
assert (
|
||||
u5_round1.yield_data.yield_matrix == ref_U5_yields.yield_matrix
|
||||
).all()
|
||||
|
||||
# Per the chain_simple.xml
|
||||
# I135 -> Xe135 -> Cs135
|
||||
# I135 -> Xe136
|
||||
# No limit on depth
|
||||
iodine_chain = gnd_simple_chain.reduce(["I135"])
|
||||
truncated_iodine = gnd_simple_chain.reduce(["I135"], 1)
|
||||
assert len(iodine_chain) == 4
|
||||
assert len(truncated_iodine) == 3
|
||||
assert set(iodine_chain.nuclide_dict) == {
|
||||
"I135", "Xe135", "Xe136", "Cs135"}
|
||||
assert set(truncated_iodine.nuclide_dict) == {"I135", "Xe135", "Xe136"}
|
||||
assert iodine_chain.reactions == ["(n,gamma)"]
|
||||
assert iodine_chain["I135"].decay_modes == ref_iodine.decay_modes
|
||||
assert iodine_chain["I135"].reactions == ref_iodine.reactions
|
||||
for mode in truncated_iodine["Xe135"].decay_modes:
|
||||
assert mode.target is None
|
||||
|
||||
# Test that no FissionYieldDistribution is made if there are no
|
||||
# fission products
|
||||
u5_noyields = gnd_simple_chain.reduce(["U235"], 0)["U235"]
|
||||
assert u5_noyields.yield_data is None
|
||||
|
||||
# Check early termination if the eventual full chain
|
||||
# is specified by using the iodine isotopes
|
||||
new_iodine = gnd_simple_chain.reduce(set(iodine_chain.nuclide_dict))
|
||||
assert set(iodine_chain.nuclide_dict) == set(new_iodine.nuclide_dict)
|
||||
|
||||
# Failure if some requested isotopes not in chain
|
||||
|
||||
with pytest.raises(IndexError, match=".*not found.*Xx999"):
|
||||
gnd_simple_chain.reduce(["U235", "Xx999"])
|
||||
|
|
|
|||
|
|
@ -176,6 +176,21 @@ def test_fission_yield_distribution():
|
|||
with pytest.raises(TypeError):
|
||||
orig_yields *= similar
|
||||
|
||||
# Test restriction of fission products
|
||||
strict_restrict = yield_dist.restrict_products(["Xe135", "Sm149"])
|
||||
with_extras = yield_dist.restrict_products(
|
||||
["Xe135", "Sm149", "H1", "U235"])
|
||||
|
||||
assert strict_restrict.products == ("Sm149", "Xe135")
|
||||
assert strict_restrict.energies == yield_dist.energies
|
||||
assert with_extras.products == ("Sm149", "Xe135")
|
||||
assert with_extras.energies == yield_dist.energies
|
||||
for ene, new_yields in strict_restrict.items():
|
||||
for product in strict_restrict.products:
|
||||
assert new_yields[product] == yield_dist[ene][product]
|
||||
assert with_extras[ene][product] == yield_dist[ene][product]
|
||||
|
||||
assert yield_dist.restrict_products(["U235"]) is None
|
||||
|
||||
def test_validate():
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue