Merge branch 'develop' into chain-fission-yields

This commit is contained in:
Andrew Johnson 2019-08-21 17:32:14 -05:00
commit 93d68a24df
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
17 changed files with 404 additions and 50 deletions

View file

@ -1,28 +0,0 @@
.. _developers:
================
Development Team
================
The following people have contributed to development of the OpenMC Monte Carlo
code:
* `Paul Romano <https://github.com/paulromano>`_
* `Bryan Herman <https://github.com/bhermanmit>`_
* `Nick Horelik <https://github.com/nhorelik>`_
* `Adam Nelson <https://github.com/nelsonag>`_
* `Jon Walsh <https://github.com/walshjon>`_
* `Sterling Harper <https://github.com/smharper>`_
* `Will Boyd <https://github.com/wbinventor>`_
* `Samuel Shaner <https://github.com/samuelshaner>`_
* `Jingang Liang <https://github.com/liangjg>`_
* `Colin Josey <https://github.com/cjosey>`_
* `Amanda Lund <https://github.com/amandalund>`_
* `Guillaume Giudicelli <https://github.com/giudgiud>`_
* `Isaac Meyer <https://github.com/icmeyer>`_
* `Patrick Shriwise <https://github.com/pshriwise>`_
* `Shikhar Kumar <https://github.com/shikhar413>`_
* `Andrew Davis <https://github.com/makeclean>`_
* `Benoit Forget <http://web.mit.edu/nse/people/faculty/forget.html>`_
* `Kord Smith <http://web.mit.edu/nse/people/faculty/smith.html>`_
* `Andrew Siegel <http://www.mcs.anl.gov/person/andrew-siegel>`_

View file

@ -44,4 +44,3 @@ list <https://groups.google.com/forum/?fromgroups=#!forum/openmc-users>`_.
io_formats/index
publications
license
developers

View file

@ -4,7 +4,7 @@
License Agreement
=================
Copyright © 2011-2018 Massachusetts Institute of Technology and OpenMC contributors
Copyright © 2011-2019 Massachusetts Institute of Technology and OpenMC contributors
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in

View file

@ -57,6 +57,9 @@ Benchmarking
Coupling and Multi-physics
--------------------------
- Miriam A. Kreher, Benoit Forget, and Kord Smith, "Single-Batch Monte Carlo
Multiphysics Coupling," *Proc. M&C*, Portland, Oregon, Aug. 25-29 (2019).
- Ze-Long Zhao, Yongwei Yang, and Shuang Hong, "`Application of FLUKA and OpenMC
in coupled physics calculation of target and subcritical reactor for ADS
<https://doi.org/10.1007/s41365-018-0539-1>`_," *Nucl. Sci. Tech.*, **30**: 10
@ -113,6 +116,10 @@ Coupling and Multi-physics
Geometry and Visualization
--------------------------
- Sterling Harper, Paul Romano, Benoit Forget, and Kord Smith, "Efficient
dynamic threadsafe neighbor lists for Monte Carlo ray tracing," *Proc. M&C*,
Portland, Oregon, Aug. 25-29 (2019).
- Jin-Yang Li, Long Gu, Hu-Shan Xu, Nadezha Korepanova, Rui Yu, Yan-Lei Zhu, and
Chang-Ping Qin, "`CAD modeling study on FLUKA and OpenMC for accelerator
driven system simulation <https://doi.org/10.1016/j.anucene.2017.12.050>`_",
@ -303,6 +310,11 @@ Multigroup Cross Section Generation
Doppler Broadening
------------------
- Jonathan A. Walsh, Benoit Forget, Kord S. Smith, and Forrest B. Brown,
"`On-the-fly Doppler broadening of unresolved resonance region cross sections
<https://doi.org/10.1016/j.pnucene.2017.05.032>`_," *Prog. Nucl. Energy*,
**101**, 444-460 (2017).
- Colin Josey, Pablo Ducru, Benoit Forget, and Kord Smith, "`Windowed multipole
for cross section Doppler broadening
<https://doi.org/10.1016/j.jcp.2015.08.013>`_," *J. Comput. Phys.*, **307**,
@ -313,6 +325,12 @@ Doppler Broadening
via Probability Band Interpolation," *Proc. PHYSOR*, Sun Valley, Idaho, May
1-5, 2016.
- Jonathan A. Walsh, Benoit Forget, Kord S. Smith, Brian C. Kiedrowski, and
Forrest B. Brown, "`Direct, on-the-fly calculation of unresolved resonance
region cross sections in Monte Carlo simulations
<http://hdl.handle.net/1721.1/108644>`_," *Proc. Joint Int. Conf. M&C+SNA+MC*,
Nashville, Tennessee, Apr. 19--23 (2015).
- Colin Josey, Benoit Forget, and Kord Smith, "`Windowed multipole sensitivity
to target accuracy of the optimization procedure
<https://doi.org/10.1080/00223131.2015.1035353>`_,"
@ -366,11 +384,6 @@ Nuclear Data
<https://doi.org/10.1016/j.cpc.2015.05.025>`_", *Comput. Phys. Commun.*,
**196**, 134-142 (2015).
- Jonathan A. Walsh, Benoit Forget, Kord S. Smith, Brian C. Kiedrowski, and
Forrest B. Brown, "Direct, on-the-fly calculation of unresolved resonance
region cross sections in Monte Carlo simulations," *Proc. Joint
Int. Conf. M&C+SNA+MC*, Nashville, Tennessee, Apr. 19--23 (2015).
- Amanda L. Lund, Andrew R. Siegel, Benoit Forget, Colin Josey, and
Paul K. Romano, "Using fractional cascading to accelerate cross section
lookups in Monte Carlo particle transport calculations," *Proc. Joint
@ -481,6 +494,10 @@ Parallelism
Depletion
---------
- Jose L. Salcedo-Perez, Benoit Forget, Kord Smith, and Paul Romano, "Hybrid
tallies to improve performance in depletion Monte Carlo simulations," *Proc.
M&C*, Aug. 25-29 (2019).
- Zhao-Qing Liu, Ze-Long Zhao, Yong-Wei Yang, Yu-Cui Gao, Hai-Yan Meng, and
Qing-Yu Gao, "`Development and validation of depletion code system IMPC-Burnup
for ADS <https://doi.org/10.1007/s41365-019-0560-z>`_," *Nucl. Sci. Tech.*,

View file

@ -106,6 +106,37 @@ Compton profile data using an existing data library from `Geant4
<http://geant4.cern.ch/>`_. Note that OpenMC includes this data file by default
so it should not be necessary in practice to generate it yourself.
.. _scripts_depletion_chain:
-------------------------------
``openmc-make-depletion-chain``
-------------------------------
This script generates a depletion chain file called ``chain_endfb71.xml``
using ENDF/B-VII.1 nuclear data. If the :envvar:`OPENMC_ENDF_DATA` variable
is not set, and ``"neutron"``, ``"decay"``, ``"nfy"`` directories
do not exist, then ENDF/B-VII.1 data will be downloaded.
.. _scripts_depletion_chain_casl:
------------------------------------
``openmc-make-depletion-chain-casl``
------------------------------------
This script generates a depletion chain called ``chain_casl.xml``
using ENDF/B-VII.1 nuclear data for a simplified chain.
The nuclides were chosen by CASL-ORIGEN, which can be found in
Appendix A of Kang Seog Kim, `"Specification for the VERA Depletion
Benchmark Suite" <https://doi.org/10.2172/1256820>`_,
CASL-U-2015-1014-000, Rev. 0, ORNL/TM-2016/53, 2016.
``Te129`` has been added into this chain due to its link to
``I129`` production.
If the :envvar:`OPENMC_ENDF_DATA` variable is not set,
and ``"neutron"``, ``"decay"``, ``"nfy"`` directories
to not exist, then ENDF/B-VII.1 data will be downloaded.
.. _scripts_stopping:
-------------------------------

View file

@ -126,6 +126,12 @@ namespace data {
extern std::array<double, 2> energy_min;
extern std::array<double, 2> energy_max;
//! Minimum temperature in [K] that nuclide data is available at
extern double temperature_min;
//! Maximum temperature in [K] that nuclide data is available at
extern double temperature_max;
extern std::vector<std::unique_ptr<Nuclide>> nuclides;
extern std::unordered_map<std::string, int> nuclide_map;

View file

@ -54,7 +54,7 @@ Indicates the default path to an HDF5 file that contains multi-group cross
section libraries if the user has not specified the <cross_sections> tag in
.I materials.xml\fP.
.SH LICENSE
Copyright \(co 2011-2018 Massachusetts Institute of Technology and OpenMC
Copyright \(co 2011-2019 Massachusetts Institute of Technology and OpenMC
contributors.
.PP
Permission is hereby granted, free of charge, to any person obtaining a copy of

View file

@ -10,9 +10,10 @@ import math
import re
from collections import OrderedDict, defaultdict
from collections.abc import Mapping, Iterable
from numbers import Real
from warnings import warn
from openmc.checkvalue import check_type, check_less_than
from openmc.checkvalue import check_type, check_less_than, check_greater_than
from openmc.data import gnd_name, zam
from .nuclide import FissionYieldDistribution
@ -666,3 +667,53 @@ class Chain(object):
yields = [yields]
check_type("fission_yields", yields, Iterable, Mapping)
self._fission_yields = yields
def validate(self, strict=True, quiet=False, tolerance=1e-4):
"""Search for possible inconsistencies
The following checks are performed for all nuclides present:
1) For all non-fission reactions, does the sum of branching
ratios equal about one?
2) For fission reactions, does the sum of fission yield
fractions equal about two?
Parameters
----------
strict : bool, optional
Raise exceptions at the first inconsistency if true.
Otherwise mark a warning
quiet : bool, optional
Flag to suppress warnings and return immediately at
the first inconsistency. Used only if
``strict`` does not evaluate to ``True``.
tolerance : float, optional
Absolute tolerance for comparisons. Used to compare computed
value ``x`` to intended value ``y`` as::
valid = (y - tolerance <= x <= y + tolerance)
Returns
-------
valid : bool
True if no inconsistencies were found
Raises
------
ValueError
If ``strict`` evaluates to ``True`` and an inconistency was
found
See Also
--------
openmc.deplete.Nuclide.validate
"""
check_type("tolerance", tolerance, Real)
check_greater_than("tolerance", tolerance, 0.0, True)
valid = True
# Sort through nuclides by name
for name in sorted(self.nuclide_dict):
stat = self[name].validate(strict, quiet, tolerance)
if quiet and not stat:
return stat
valid = valid and stat
return valid

View file

@ -3,9 +3,9 @@
Contains the per-nuclide components of a depletion chain.
"""
from collections import namedtuple
from collections.abc import Mapping
from collections import namedtuple, defaultdict
from warnings import warn
try:
import lxml.etree as ET
except ImportError:
@ -241,6 +241,103 @@ class Nuclide(object):
return elem
def validate(self, strict=True, quiet=False, tolerance=1e-4):
"""Search for possible inconsistencies
The following checks are performed:
1) for all non-fission reactions and decay modes,
does the sum of branching ratios equal about one?
2) for fission reactions, does the sum of fission yield
fractions equal about two?
Parameters
----------
strict : bool, optional
Raise exceptions at the first inconsistency if true.
Otherwise mark a warning
quiet : bool, optional
Flag to suppress warnings and return immediately at
the first inconsistency. Used only if
``strict`` does not evaluate to ``True``.
tolerance : float, optional
Absolute tolerance for comparisons. Used to compare computed
value ``x`` to intended value ``y`` as::
valid = (y - tolerance <= x <= y + tolerance)
Returns
-------
valid : bool
True if no inconsistencies were found
Raises
------
ValueError
If ``strict`` evaluates to ``True`` and an inconistency was
found
See Also
--------
openmc.deplete.Chain.validate
"""
msg_func = ("Nuclide {name} has {prop} that sum to {actual} "
"instead of {expected} +/- {tol:7.4e}").format
valid = True
# check decay modes
if self.decay_modes:
sum_br = sum(m.branching_ratio for m in self.decay_modes)
stat = 1.0 - tolerance <= sum_br <= 1.0 + tolerance
if not stat:
msg = msg_func(
name=self.name, actual=sum_br, expected=1.0, tol=tolerance,
prop="decay mode branch ratios")
if strict:
raise ValueError(msg)
elif quiet:
return False
warn(msg)
valid = False
if self.reactions:
type_map = defaultdict(set)
for reaction in self.reactions:
type_map[reaction.type].add(reaction)
for rxn_type, reactions in type_map.items():
sum_rxn = sum(rx.branching_ratio for rx in reactions)
stat = 1.0 - tolerance <= sum_rxn <= 1.0 + tolerance
if stat:
continue
msg = msg_func(
name=self.name, actual=sum_br, expected=1.0, tol=tolerance,
prop="{} reaction branch ratios".format(rxn_type))
if strict:
raise ValueError(msg)
elif quiet:
return False
warn(msg)
valid = False
if self.yield_data:
for energy, yield_list in self.yield_data.items():
sum_yield = sum(y[1] for y in yield_list)
stat = 2.0 - tolerance <= sum_yield <= 2.0 + tolerance
if stat:
continue
msg = msg_func(
name=self.name, actual=sum_yield,
expected=2.0, tol=tolerance,
prop="fission yields (E = {:7.4e} eV)".format(energy))
if strict:
raise ValueError(msg)
elif quiet:
return False
warn(msg)
valid = False
return valid
class FissionYieldDistribution(Mapping):
"""Energy-dependent fission product yields for a single nuclide

View file

@ -6,6 +6,8 @@
# Note 32 of the 255 nuclides appeare twice as they are both activation
# nuclides (category 1) and fission product nuclides (category 3).
# Te129 has been added due to it's link to I129 production.
CASL_CHAIN = {
# Nuclide: (Stable, CAT, IFPY, Special yield treatment)
# Stable: True if nuclide has no decay reactions
@ -187,6 +189,7 @@ CASL_CHAIN = {
'Sb127': (False, 3, 2, None),
'Te127': (False, 3, -1, None),
'Te127_m1': (False, 3, -1, None),
'Te129': (False, 3, 1, None),
'Te129_m1': (False, 3, 2, None),
'Te132': (False, 3, 2, None),
'I127': (True, 3, 1, None),

View file

@ -1,7 +1,7 @@
#!/usr/bin/env python3
import glob
import os
from pathlib import Path
from zipfile import ZipFile
from openmc._utils import download
@ -15,15 +15,28 @@ URLS = [
]
def main():
for url in URLS:
basename = download(url)
with ZipFile(basename, 'r') as zf:
print('Extracting {}...'.format(basename))
zf.extractall()
endf_dir = os.environ.get("OPENMC_ENDF_DATA")
if endf_dir is not None:
endf_dir = Path(endf_dir)
elif all(os.path.isdir(lib) for lib in ("neutrons", "decay", "nfy")):
endf_dir = Path(".")
else:
for url in URLS:
basename = download(url)
with ZipFile(basename, 'r') as zf:
print('Extracting {}...'.format(basename))
zf.extractall()
endf_dir = Path(".")
decay_files = glob.glob(os.path.join('decay', '*.endf'))
nfy_files = glob.glob(os.path.join('nfy', '*.endf'))
neutron_files = glob.glob(os.path.join('neutrons', '*.endf'))
decay_files = tuple((endf_dir / "decay").glob("*endf"))
neutron_files = tuple((endf_dir / "neutrons").glob("*endf"))
nfy_files = tuple((endf_dir / "nfy").glob("*endf"))
# check files exist
for flist, ftype in [(decay_files, "decay"), (neutron_files, "neutron"),
(nfy_files, "neutron fission product yield")]:
if not flist:
raise IOError("No {} endf files found in {}".format(ftype, endf_dir))
chain = openmc.deplete.Chain.from_endf(decay_files, nfy_files, neutron_files)
chain.export_to_xml('chain_endfb71.xml')

View file

@ -16,6 +16,7 @@
#include "openmc/hdf5_interface.h"
#include "openmc/lattice.h"
#include "openmc/material.h"
#include "openmc/nuclide.h"
#include "openmc/settings.h"
#include "openmc/surface.h"
#include "openmc/xml_interface.h"
@ -244,6 +245,16 @@ Cell::temperature(int32_t instance) const
void
Cell::set_temperature(double T, int32_t instance)
{
if (settings::temperature_method == TEMPERATURE_INTERPOLATION) {
if (T < data::temperature_min) {
throw std::runtime_error{"Temperature is below minimum temperature at "
"which data is available."};
} else if (T > data::temperature_max) {
throw std::runtime_error{"Temperature is above maximum temperature at "
"which data is available."};
}
}
if (instance >= 0) {
sqrtkT_.at(instance) = std::sqrt(K_BOLTZMANN * T);
} else {
@ -1057,7 +1068,6 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
return 0;
}
//TODO: make sure data is loaded for this temperature
extern "C" int
openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance)
{

View file

@ -357,6 +357,12 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
}
}
// Show minimum/maximum temperature
write_message("Minimum neutron data temperature: " +
std::to_string(data::temperature_min) + " K");
write_message("Maximum neutron data temperature: " +
std::to_string(data::temperature_max) + " K");
// If the user wants multipole, make sure we found a multipole library.
if (settings::temperature_multipole) {
bool mp_found = false;

View file

@ -117,6 +117,8 @@ int openmc_finalize()
data::energy_max = {INFTY, INFTY};
data::energy_min = {0.0, 0.0};
data::temperature_min = 0.0;
data::temperature_max = INFTY;
model::root_universe = -1;
openmc::openmc_set_seed(DEFAULT_SEED);

View file

@ -18,7 +18,7 @@
#include "xtensor/xbuilder.hpp"
#include "xtensor/xview.hpp"
#include <algorithm> // for sort
#include <algorithm> // for sort, min_element
#include <string> // for to_string, stoi
namespace openmc {
@ -30,6 +30,8 @@ namespace openmc {
namespace data {
std::array<double, 2> energy_min {0.0, 0.0};
std::array<double, 2> energy_max {INFTY, INFTY};
double temperature_min {0.0};
double temperature_max {INFTY};
std::vector<std::unique_ptr<Nuclide>> nuclides;
std::unordered_map<std::string, int> nuclide_map;
} // namespace data
@ -154,6 +156,12 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature, int i_nucl
// Sort temperatures to read
std::sort(temps_to_read.begin(), temps_to_read.end());
double T_min_read = *std::min_element(temps_to_read.cbegin(), temps_to_read.cend());
double T_max_read = *std::max_element(temps_to_read.cbegin(), temps_to_read.cend());
data::temperature_min = std::max(data::temperature_min, T_min_read);
data::temperature_max = std::min(data::temperature_max, T_max_read);
hid_t energy_group = open_group(group, "energy");
for (const auto& T : temps_to_read) {
std::string dset {std::to_string(T) + "K"};

View file

@ -360,3 +360,54 @@ def test_fission_yield_attribute(simple_chain):
with pytest.raises(
ValueError, match="fission yield.*not equal.*compositions"):
cram.deplete(empty_chain, dummy_conc, None, 0.5)
def test_validate(simple_chain):
"""Test the validate method"""
# current chain is invalid
# fission yields do not sum to 2.0
with pytest.raises(ValueError, match="Nuclide C.*fission yields"):
simple_chain.validate(strict=True, tolerance=0.0)
with pytest.warns(UserWarning) as record:
assert not simple_chain.validate(strict=False, quiet=False, tolerance=0.0)
assert not simple_chain.validate(strict=False, quiet=True, tolerance=0.0)
assert len(record) == 1
assert "Nuclide C" in record[0].message.args[0]
# Fix fission yields but keep to restore later
old_yields = simple_chain["C"].yield_data
simple_chain["C"].yield_data = {0.0253: [("A", 1.4), ("B", 0.6)]}
assert simple_chain.validate(strict=True, tolerance=0.0)
with pytest.warns(None) as record:
assert simple_chain.validate(strict=False, quiet=False, tolerance=0.0)
assert len(record) == 0
# Mess up "earlier" nuclide's reactions
decay_mode = simple_chain["A"].decay_modes.pop()
with pytest.raises(ValueError, match="Nuclide A.*decay mode"):
simple_chain.validate(strict=True, tolerance=0.0)
# restore old fission yields
simple_chain["C"].yield_data = old_yields
with pytest.warns(UserWarning) as record:
assert not simple_chain.validate(strict=False, quiet=False, tolerance=0.0)
assert len(record) == 2
assert "Nuclide A" in record[0].message.args[0]
assert "Nuclide C" in record[1].message.args[0]
# restore decay modes
simple_chain["A"].decay_modes.append(decay_mode)
def test_validate_inputs():
c = Chain()
with pytest.raises(TypeError, match="tolerance"):
c.validate(tolerance=None)
with pytest.raises(ValueError, match="tolerance"):
c.validate(tolerance=-1)

View file

@ -153,3 +153,91 @@ def test_fission_yield_distribution():
assert numpy.array_equal(orig_yield_obj.yields * 2, mod_yields.yields)
mod_yields += orig_yield_obj
assert numpy.array_equal(orig_yield_obj.yields * 3, mod_yields.yields)
def test_validate():
nuc = nuclide.Nuclide()
nuc.name = "Test"
# decay modes: type, target, branching_ratio
nuc.decay_modes = [
nuclide.DecayTuple("type 0", "0", 0.5),
nuclide.DecayTuple("type 1", "1", 0.5),
]
# reactions: type, target, Q, branching_ratio
nuc.reactions = [
nuclide.ReactionTuple("0", "0", 1000, 0.3),
nuclide.ReactionTuple("0", "1", 1000, 0.3),
nuclide.ReactionTuple("1", "2", 1000, 1.0),
nuclide.ReactionTuple("0", "3", 1000, 0.4),
]
# fission yields
nuc.yield_data = {
0.0253: [("0", 1.5), ("1", 0.5)],
1e6: [("0", 1.5), ("1", 0.5)],
}
# nuclide is good and should have no warnings raise
with pytest.warns(None) as record:
assert nuc.validate(strict=True, quiet=False, tolerance=0.0)
assert len(record) == 0
# invalidate decay modes
decay = nuc.decay_modes.pop()
with pytest.raises(ValueError, match="decay mode"):
nuc.validate(strict=True, quiet=False, tolerance=0.0)
with pytest.warns(UserWarning) as record:
assert not nuc.validate(strict=False, quiet=False, tolerance=0.0)
assert not nuc.validate(strict=False, quiet=True, tolerance=0.0)
assert len(record) == 1
assert "decay mode" in record[0].message.args[0]
# restore decay modes, invalidate reactions
nuc.decay_modes.append(decay)
reaction = nuc.reactions.pop()
with pytest.raises(ValueError, match="0 reaction"):
nuc.validate(strict=True, quiet=False, tolerance=0.0)
with pytest.warns(UserWarning) as record:
assert not nuc.validate(strict=False, quiet=False, tolerance=0.0)
assert not nuc.validate(strict=False, quiet=True, tolerance=0.0)
assert len(record) == 1
assert "0 reaction" in record[0].message.args[0]
# restore reactions, invalidate fission yields
nuc.reactions.append(reaction)
nuc.yield_data[1e6].pop()
with pytest.raises(ValueError, match=r"fission yields.*1\.0*e"):
nuc.validate(strict=True, quiet=False, tolerance=0.0)
with pytest.warns(UserWarning) as record:
assert not nuc.validate(strict=False, quiet=False, tolerance=0.0)
assert not nuc.validate(strict=False, quiet=True, tolerance=0.0)
assert len(record) == 1
assert "1.0" in record[0].message.args[0]
# invalidate everything, check that error is raised at decay modes
decay = nuc.decay_modes.pop()
reaction = nuc.reactions.pop()
with pytest.raises(ValueError, match="decay mode"):
nuc.validate(strict=True, quiet=False, tolerance=0.0)
# check for warnings
# should be one warning for decay modes, reactions, fission yields
with pytest.warns(UserWarning) as record:
assert not nuc.validate(strict=False, quiet=False, tolerance=0.0)
assert not nuc.validate(strict=False, quiet=True, tolerance=0.0)
assert len(record) == 3
assert "decay mode" in record[0].message.args[0]
assert "0 reaction" in record[1].message.args[0]
assert "1.0" in record[2].message.args[0]