mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-29 06:35:48 -04:00
Merge pull request #976 from paulromano/deplete
Implement depletion via Python
This commit is contained in:
Colin Josey
GitHub
commit
d29bc1e1d7
55 changed files with 4714 additions and 58 deletions
7
.gitignore
vendored
7
.gitignore
vendored
|
|
@ -42,11 +42,8 @@ results_error.dat
|
|||
inputs_error.dat
|
||||
results_test.dat
|
||||
|
||||
# Test build files
|
||||
tests/build/
|
||||
tests/coverage/
|
||||
tests/memcheck/
|
||||
tests/ctestscript.run
|
||||
# Test
|
||||
.pytest_cache/
|
||||
|
||||
# HDF5 files
|
||||
*.h5
|
||||
|
|
|
|||
|
|
@ -496,7 +496,8 @@ add_custom_command(TARGET libopenmc POST_BUILD
|
|||
|
||||
install(TARGETS ${program} libopenmc
|
||||
RUNTIME DESTINATION bin
|
||||
LIBRARY DESTINATION lib)
|
||||
LIBRARY DESTINATION lib
|
||||
ARCHIVE DESTINATION lib)
|
||||
install(DIRECTORY src/relaxng DESTINATION share/openmc)
|
||||
install(FILES man/man1/openmc.1 DESTINATION share/man/man1)
|
||||
install(FILES LICENSE DESTINATION "share/doc/${program}" RENAME copyright)
|
||||
|
|
|
|||
|
|
@ -21,15 +21,18 @@ on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
|
|||
from unittest.mock import MagicMock
|
||||
|
||||
|
||||
MOCK_MODULES = ['numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial',
|
||||
'numpy.ctypeslib', 'scipy', 'scipy.sparse', 'scipy.interpolate',
|
||||
'scipy.integrate', 'scipy.optimize', 'scipy.special',
|
||||
'scipy.stats', 'scipy.spatial', 'h5py', 'pandas', 'uncertainties',
|
||||
'matplotlib', 'matplotlib.pyplot','openmoc',
|
||||
'openmc.data.reconstruct']
|
||||
MOCK_MODULES = [
|
||||
'numpy', 'numpy.polynomial', 'numpy.polynomial.polynomial',
|
||||
'numpy.ctypeslib', 'scipy', 'scipy.sparse', 'scipy.sparse.linalg',
|
||||
'scipy.interpolate', 'scipy.integrate', 'scipy.optimize', 'scipy.special',
|
||||
'scipy.stats', 'scipy.spatial', 'h5py', 'pandas', 'uncertainties',
|
||||
'matplotlib', 'matplotlib.pyplot', 'openmoc',
|
||||
'openmc.data.reconstruct'
|
||||
]
|
||||
sys.modules.update((mod_name, MagicMock()) for mod_name in MOCK_MODULES)
|
||||
|
||||
import numpy as np
|
||||
np.ndarray = MagicMock
|
||||
np.polynomial.Polynomial = MagicMock
|
||||
|
||||
|
||||
|
|
|
|||
102
docs/source/io_formats/depletion_chain.rst
Normal file
102
docs/source/io_formats/depletion_chain.rst
Normal file
|
|
@ -0,0 +1,102 @@
|
|||
.. _io_depletion_chain:
|
||||
|
||||
============================
|
||||
Depletion Chain -- chain.xml
|
||||
============================
|
||||
|
||||
A depletion chain file has a ``<depletion_chain>`` root element with one or more
|
||||
``<nuclide>`` child elements. The decay, reaction, and fission product data for
|
||||
each nuclide appears as child elements of ``<nuclide>``.
|
||||
|
||||
---------------------
|
||||
``<nuclide>`` Element
|
||||
---------------------
|
||||
|
||||
The ``<nuclide>`` element contains information on the decay modes, reactions,
|
||||
and fission product yields for a given nuclide in the depletion chain. This
|
||||
element may have the following attributes:
|
||||
|
||||
:name:
|
||||
Name of the nuclide
|
||||
|
||||
:half_life:
|
||||
Half-life of the nuclide in [s]
|
||||
|
||||
:decay_modes:
|
||||
Number of decay modes present
|
||||
|
||||
:decay_energy:
|
||||
Decay energy released in [eV]
|
||||
|
||||
:reactions:
|
||||
Number of reactions present
|
||||
|
||||
For each decay mode, a :ref:`io_chain_decay` appears as a child of
|
||||
``<nuclide>``. For each reaction present, a :ref:`io_chain_reaction` appears as
|
||||
a child of ``<nuclide>``. If the nuclide is fissionable, a :ref:`io_chain_nfy`
|
||||
appears as well.
|
||||
|
||||
.. _io_chain_decay:
|
||||
|
||||
-------------------
|
||||
``<decay>`` Element
|
||||
-------------------
|
||||
|
||||
The ``<decay>`` element represents a single decay mode and has the following
|
||||
attributes:
|
||||
|
||||
:type:
|
||||
The type of the decay, e.g. 'ec/beta+'
|
||||
|
||||
:target:
|
||||
The daughter nuclide produced from the decay
|
||||
|
||||
:branching_ratio:
|
||||
The branching ratio for this decay mode
|
||||
|
||||
.. _io_chain_reaction:
|
||||
|
||||
----------------------
|
||||
``<reaction>`` Element
|
||||
----------------------
|
||||
|
||||
The ``<reaction>`` element represents a single transmutation reaction. This
|
||||
element has the following attributes:
|
||||
|
||||
:type:
|
||||
The type of the reaction, e.g., '(n,gamma)'
|
||||
|
||||
:Q:
|
||||
The Q value of the reaction in [eV]
|
||||
|
||||
:target:
|
||||
The nuclide produced in the reaction (absent if the type is 'fission')
|
||||
|
||||
:branching_ratio:
|
||||
The branching ratio for the reaction
|
||||
|
||||
.. _io_chain_nfy:
|
||||
|
||||
------------------------------------
|
||||
``<neutron_fission_yields>`` Element
|
||||
------------------------------------
|
||||
|
||||
The ``<neutron_fission_yields>`` element provides yields of fission products for
|
||||
fissionable nuclides. It has the follow sub-elements:
|
||||
|
||||
:energies:
|
||||
Energies in [eV] at which yields for products are tabulated
|
||||
|
||||
:fission_yields:
|
||||
|
||||
Fission product yields for a single energy point. This element itself has a
|
||||
number of attributes/sub-elements:
|
||||
|
||||
:energy:
|
||||
Energy in [eV] at which yields are tabulated
|
||||
|
||||
:products:
|
||||
Names of fission products
|
||||
|
||||
:data:
|
||||
Independent yields for each fission product
|
||||
42
docs/source/io_formats/depletion_results.rst
Normal file
42
docs/source/io_formats/depletion_results.rst
Normal file
|
|
@ -0,0 +1,42 @@
|
|||
.. _io_depletion_results:
|
||||
|
||||
=============================
|
||||
Depletion Results File Format
|
||||
=============================
|
||||
|
||||
The current version of the depletion results file format is 1.0.
|
||||
|
||||
**/**
|
||||
|
||||
:Attributes: - **filetype** (*char[]*) -- String indicating the type of file.
|
||||
- **version** (*int[2]*) -- Major and minor version of the
|
||||
statepoint file format.
|
||||
|
||||
:Datasets: - **eigenvalues** (*double[][]*) -- k-eigenvalues at each
|
||||
time/stage. This array has shape (number of timesteps, number of
|
||||
stages).
|
||||
- **number** (*double[][][][]*) -- Total number of atoms. This array
|
||||
has shape (number of timesteps, number of stages, number of
|
||||
materials, number of nuclides).
|
||||
- **reaction rates** (*double[][][][][]*) -- Reaction rates used to
|
||||
build depletion matrices. This array has shape (number of
|
||||
timesteps, number of stages, number of materials, number of
|
||||
nuclides, number of reactions).
|
||||
- **time** (*double[][2]*) -- Time in [s] at beginning/end of each
|
||||
step.
|
||||
|
||||
**/materials/<id>/**
|
||||
|
||||
:Attributes: - **index** (*int*) -- Index used in results for this material
|
||||
- **volume** (*double*) -- Volume of this material in [cm^3]
|
||||
|
||||
**/nuclides/<name>/**
|
||||
|
||||
:Attributes: - **atom number index** (*int*) -- Index in array of total atoms
|
||||
for this nuclide
|
||||
- **reaction rate index** (*int*) -- Index in array of reaction
|
||||
rates for this nuclide
|
||||
|
||||
**/reactions/<name>/**
|
||||
|
||||
:Attributes: - **index** (*int*) -- Index user in results for this reaction
|
||||
|
|
@ -12,7 +12,7 @@ Input Files
|
|||
|
||||
.. toctree::
|
||||
:numbered:
|
||||
:maxdepth: 2
|
||||
:maxdepth: 1
|
||||
|
||||
geometry
|
||||
materials
|
||||
|
|
@ -27,9 +27,10 @@ Data Files
|
|||
|
||||
.. toctree::
|
||||
:numbered:
|
||||
:maxdepth: 2
|
||||
:maxdepth: 1
|
||||
|
||||
cross_sections
|
||||
depletion_chain
|
||||
nuclear_data
|
||||
mgxs_library
|
||||
data_wmp
|
||||
|
|
@ -41,11 +42,12 @@ Output Files
|
|||
|
||||
.. toctree::
|
||||
:numbered:
|
||||
:maxdepth: 2
|
||||
:maxdepth: 1
|
||||
|
||||
statepoint
|
||||
source
|
||||
summary
|
||||
depletion_results
|
||||
particle_restart
|
||||
track
|
||||
voxel
|
||||
|
|
|
|||
|
|
@ -32,10 +32,12 @@ Core Functions
|
|||
:template: myfunction.rst
|
||||
|
||||
openmc.data.atomic_mass
|
||||
openmc.data.gnd_name
|
||||
openmc.data.linearize
|
||||
openmc.data.thin
|
||||
openmc.data.water_density
|
||||
openmc.data.write_compact_458_library
|
||||
openmc.data.zam
|
||||
|
||||
Angle-Energy Distributions
|
||||
--------------------------
|
||||
|
|
|
|||
85
docs/source/pythonapi/deplete.rst
Normal file
85
docs/source/pythonapi/deplete.rst
Normal file
|
|
@ -0,0 +1,85 @@
|
|||
.. _pythonapi_deplete:
|
||||
|
||||
----------------------------------
|
||||
:mod:`openmc.deplete` -- Depletion
|
||||
----------------------------------
|
||||
|
||||
.. module:: openmc.deplete
|
||||
|
||||
Two functions are provided that implement different time-integration algorithms
|
||||
for depletion calculations.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated
|
||||
:nosignatures:
|
||||
:template: myfunction.rst
|
||||
|
||||
integrator.predictor
|
||||
integrator.cecm
|
||||
|
||||
Each of these functions expects a "transport operator" to be passed. An operator
|
||||
specific to OpenMC is available using the following class:
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated
|
||||
:nosignatures:
|
||||
:template: myclass.rst
|
||||
|
||||
Operator
|
||||
|
||||
When running in parallel using `mpi4py <http://mpi4py.scipy.org>`_, the MPI
|
||||
intercommunicator used can be changed by modifying the following module
|
||||
variable. If it is not explicitly modified, it defaults to
|
||||
``mpi4py.MPI.COMM_WORLD``.
|
||||
|
||||
.. data:: comm
|
||||
|
||||
MPI intercommunicator used to call OpenMC library
|
||||
|
||||
:type: mpi4py.MPI.Comm
|
||||
|
||||
Internal Classes and Functions
|
||||
------------------------------
|
||||
|
||||
During a depletion calculation, the depletion chain, reaction rates, and number
|
||||
densities are managed through a series of internal classes that are not normally
|
||||
visible to a user. However, should you find yourself wondering about these
|
||||
classes (e.g., if you want to know what decay modes or reactions are present in
|
||||
a depletion chain), they are documented here. The following classes store data
|
||||
for a depletion chain:
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated
|
||||
:nosignatures:
|
||||
:template: myclass.rst
|
||||
|
||||
Chain
|
||||
DecayTuple
|
||||
Nuclide
|
||||
ReactionTuple
|
||||
|
||||
The following classes are used during a depletion simulation and store auxiliary
|
||||
data, such as number densities and reaction rates for each material.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated
|
||||
:nosignatures:
|
||||
:template: myclass.rst
|
||||
|
||||
AtomNumber
|
||||
OperatorResult
|
||||
ReactionRates
|
||||
Results
|
||||
ResultsList
|
||||
TransportOperator
|
||||
|
||||
Each of the integrator functions also relies on a number of "helper" functions
|
||||
as follows:
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated
|
||||
:nosignatures:
|
||||
:template: myfunction.rst
|
||||
|
||||
integrator.CRAM16
|
||||
integrator.CRAM48
|
||||
|
|
@ -15,14 +15,15 @@ there are many substantial benefits to using the Python API, including:
|
|||
- The ability to define dimensions using variables.
|
||||
- Availability of standard-library modules for working with files.
|
||||
- An entire ecosystem of third-party packages for scientific computing.
|
||||
- Ability to create materials based on natural elements or uranium enrichment
|
||||
- Automated multi-group cross section generation (:mod:`openmc.mgxs`)
|
||||
- A fully-featured nuclear data interface (:mod:`openmc.data`)
|
||||
- Depletion capability (:mod:`openmc.deplete`)
|
||||
- Convenience functions (e.g., a function returning a hexagonal region)
|
||||
- Ability to plot individual universes as geometry is being created
|
||||
- A :math:`k_\text{eff}` search function (:func:`openmc.search_for_keff`)
|
||||
- Random sphere packing for generating TRISO particle locations
|
||||
(:func:`openmc.model.pack_trisos`)
|
||||
- A fully-featured nuclear data interface (:mod:`openmc.data`)
|
||||
- Ability to create materials based on natural elements or uranium enrichment
|
||||
|
||||
For those new to Python, there are many good tutorials available online. We
|
||||
recommend going through the modules from `Codecademy
|
||||
|
|
@ -45,6 +46,7 @@ Modules
|
|||
base
|
||||
model
|
||||
examples
|
||||
deplete
|
||||
mgxs
|
||||
stats
|
||||
data
|
||||
|
|
|
|||
|
|
@ -452,6 +452,11 @@ distributions.
|
|||
.. admonition:: Optional
|
||||
:class: note
|
||||
|
||||
`mpi4py <http://mpi4py.scipy.org/>`_
|
||||
mpi4py provides Python bindings to MPI for running distributed-memory
|
||||
parallel runs. This package is needed if you plan on running depletion
|
||||
simulations in parallel using MPI.
|
||||
|
||||
`Cython <http://cython.org/>`_
|
||||
Cython is used for resonance reconstruction for ENDF data converted to
|
||||
:class:`openmc.data.IncidentNeutron`.
|
||||
|
|
|
|||
49
openmc/_utils.py
Normal file
49
openmc/_utils.py
Normal file
|
|
@ -0,0 +1,49 @@
|
|||
import os.path
|
||||
from pathlib import Path
|
||||
from urllib.parse import urlparse
|
||||
from urllib.request import urlopen
|
||||
|
||||
_BLOCK_SIZE = 16384
|
||||
|
||||
|
||||
def download(url):
|
||||
"""Download file from a URL
|
||||
|
||||
Parameters
|
||||
----------
|
||||
url : str
|
||||
URL from which to download
|
||||
|
||||
Returns
|
||||
-------
|
||||
basename : str
|
||||
Name of file written locally
|
||||
|
||||
"""
|
||||
req = urlopen(url)
|
||||
|
||||
# Get file size from header
|
||||
file_size = req.length
|
||||
|
||||
# Check if file already downloaded
|
||||
basename = Path(urlparse(url).path).name
|
||||
if os.path.exists(basename):
|
||||
if os.path.getsize(basename) == file_size:
|
||||
print('Skipping {}, already downloaded'.format(basename))
|
||||
return basename
|
||||
|
||||
# Copy file to disk in chunks
|
||||
print('Downloading {}... '.format(basename), end='')
|
||||
downloaded = 0
|
||||
with open(basename, 'wb') as fh:
|
||||
while True:
|
||||
chunk = req.read(_BLOCK_SIZE)
|
||||
if not chunk:
|
||||
break
|
||||
fh.write(chunk)
|
||||
downloaded += len(chunk)
|
||||
status = '{:10} [{:3.2f}%]'.format(
|
||||
downloaded, downloaded * 100. / file_size)
|
||||
print(status + '\b'*len(status), end='')
|
||||
print('')
|
||||
return basename
|
||||
|
|
@ -15,7 +15,6 @@ objects in the :mod:`openmc.capi` subpackage, for example:
|
|||
from ctypes import CDLL
|
||||
import os
|
||||
import sys
|
||||
from warnings import warn
|
||||
|
||||
import pkg_resources
|
||||
|
||||
|
|
@ -36,10 +35,7 @@ else:
|
|||
# available. Instead, we create a mock object so that when the modules
|
||||
# within the openmc.capi package try to configure arguments and return
|
||||
# values for symbols, no errors occur
|
||||
try:
|
||||
from unittest.mock import Mock
|
||||
except ImportError:
|
||||
from mock import Mock
|
||||
from unittest.mock import Mock
|
||||
_dll = Mock()
|
||||
|
||||
from .error import *
|
||||
|
|
@ -50,6 +46,3 @@ from .cell import *
|
|||
from .filter import *
|
||||
from .tally import *
|
||||
from .settings import settings
|
||||
|
||||
warn("The Python bindings to OpenMC's C API are still unstable "
|
||||
"and may change substantially in future releases.", FutureWarning)
|
||||
|
|
|
|||
|
|
@ -246,9 +246,9 @@ def check_filetype_version(obj, expected_type, expected_version):
|
|||
----------
|
||||
obj : h5py.File
|
||||
HDF5 file to check
|
||||
expected_type
|
||||
expected_type : str
|
||||
Expected file type, e.g. 'statepoint'
|
||||
expected_version
|
||||
expected_version : int
|
||||
Expected major version number.
|
||||
|
||||
"""
|
||||
|
|
|
|||
|
|
@ -313,14 +313,63 @@ def water_density(temperature, pressure=0.1013):
|
|||
return coeff / pi / gamma1_pi
|
||||
|
||||
|
||||
def gnd_name(Z, A, m=0):
|
||||
"""Return nuclide name using GND convention
|
||||
|
||||
Parameters
|
||||
----------
|
||||
Z : int
|
||||
Atomic number
|
||||
A : int
|
||||
Mass number
|
||||
m : int, optional
|
||||
Metastable state
|
||||
|
||||
Returns
|
||||
-------
|
||||
str
|
||||
Nuclide name in GND convention, e.g., 'Am242_m1'
|
||||
|
||||
"""
|
||||
if m > 0:
|
||||
return '{}{}_m{}'.format(ATOMIC_SYMBOL[Z], A, m)
|
||||
else:
|
||||
return '{}{}'.format(ATOMIC_SYMBOL[Z], A)
|
||||
|
||||
|
||||
def zam(name):
|
||||
"""Return tuple of (atomic number, mass number, metastable state)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
name : str
|
||||
Name of nuclide using GND convention, e.g., 'Am242_m1'
|
||||
|
||||
Returns
|
||||
-------
|
||||
3-tuple of int
|
||||
Atomic number, mass number, and metastable state
|
||||
|
||||
"""
|
||||
try:
|
||||
symbol, A, state = re.match(r'([A-Zn][a-z]*)(\d+)((?:_[em]\d+)?)',
|
||||
name).groups()
|
||||
except AttributeError:
|
||||
raise ValueError("'{}' does not appear to be a nuclide name in GND "
|
||||
"format.".format(name))
|
||||
metastable = int(state[2:]) if state else 0
|
||||
return (ATOMIC_NUMBER[symbol], int(A), metastable)
|
||||
|
||||
|
||||
# Values here are from the Committee on Data for Science and Technology
|
||||
# (CODATA) 2014 recommendation (doi:10.1103/RevModPhys.88.035009).
|
||||
|
||||
# The value of the Boltzman constant in units of eV / K
|
||||
K_BOLTZMANN = 8.6173303e-5
|
||||
|
||||
# Used for converting units in ACE data
|
||||
# Unit conversions
|
||||
EV_PER_MEV = 1.0e6
|
||||
JOULE_PER_EV = 1.6021766208e-19
|
||||
|
||||
# Avogadro's constant
|
||||
AVOGADRO = 6.022140857e23
|
||||
|
|
|
|||
|
|
@ -457,6 +457,7 @@ class Decay(EqualityMixin):
|
|||
items, values = get_list_record(file_obj)
|
||||
self.nuclide['spin'] = items[0]
|
||||
self.nuclide['parity'] = items[1]
|
||||
self.half_life = ufloat(float('inf'), float('inf'))
|
||||
|
||||
@property
|
||||
def decay_constant(self):
|
||||
|
|
|
|||
|
|
@ -10,13 +10,14 @@ import io
|
|||
import re
|
||||
import os
|
||||
from math import pi
|
||||
from pathlib import PurePath
|
||||
from collections import OrderedDict
|
||||
from collections.abc import Iterable
|
||||
|
||||
import numpy as np
|
||||
from numpy.polynomial.polynomial import Polynomial
|
||||
|
||||
from .data import ATOMIC_SYMBOL
|
||||
from .data import ATOMIC_SYMBOL, gnd_name
|
||||
from .function import Tabulated1D, INTERPOLATION_SCHEME
|
||||
from openmc.stats.univariate import Uniform, Tabular, Legendre
|
||||
|
||||
|
|
@ -248,6 +249,7 @@ def get_tab2_record(file_obj):
|
|||
|
||||
return params, Tabulated2D(breakpoints, interpolation)
|
||||
|
||||
|
||||
def get_evaluations(filename):
|
||||
"""Return a list of all evaluations within an ENDF file.
|
||||
|
||||
|
|
@ -299,8 +301,8 @@ class Evaluation(object):
|
|||
|
||||
"""
|
||||
def __init__(self, filename_or_obj):
|
||||
if isinstance(filename_or_obj, str):
|
||||
fh = open(filename_or_obj, 'r')
|
||||
if isinstance(filename_or_obj, (str, PurePath)):
|
||||
fh = open(str(filename_or_obj), 'r')
|
||||
else:
|
||||
fh = filename_or_obj
|
||||
self.section = {}
|
||||
|
|
@ -423,13 +425,9 @@ class Evaluation(object):
|
|||
|
||||
@property
|
||||
def gnd_name(self):
|
||||
symbol = ATOMIC_SYMBOL[self.target['atomic_number']]
|
||||
A = self.target['mass_number']
|
||||
m = self.target['isomeric_state']
|
||||
if m > 0:
|
||||
return '{}{}_m{}'.format(symbol, A, m)
|
||||
else:
|
||||
return '{}{}'.format(symbol, A)
|
||||
return gnd_name(self.target['atomic_number'],
|
||||
self.target['mass_number'],
|
||||
self.target['isomeric_state'])
|
||||
|
||||
|
||||
class Tabulated2D(object):
|
||||
|
|
|
|||
24
openmc/deplete/__init__.py
Normal file
24
openmc/deplete/__init__.py
Normal file
|
|
@ -0,0 +1,24 @@
|
|||
"""
|
||||
openmc.deplete
|
||||
==============
|
||||
|
||||
A depletion front-end tool.
|
||||
"""
|
||||
|
||||
from .dummy_comm import DummyCommunicator
|
||||
try:
|
||||
from mpi4py import MPI
|
||||
comm = MPI.COMM_WORLD
|
||||
have_mpi = True
|
||||
except ImportError:
|
||||
comm = DummyCommunicator()
|
||||
have_mpi = False
|
||||
|
||||
from .nuclide import *
|
||||
from .chain import *
|
||||
from .operator import *
|
||||
from .reaction_rates import *
|
||||
from .abc import *
|
||||
from .results import *
|
||||
from .results_list import *
|
||||
from .integrator import *
|
||||
142
openmc/deplete/abc.py
Normal file
142
openmc/deplete/abc.py
Normal file
|
|
@ -0,0 +1,142 @@
|
|||
"""function module.
|
||||
|
||||
This module contains the Operator class, which is then passed to an integrator
|
||||
to run a full depletion simulation.
|
||||
"""
|
||||
|
||||
from collections import namedtuple
|
||||
import os
|
||||
from pathlib import Path
|
||||
from abc import ABCMeta, abstractmethod
|
||||
|
||||
from .chain import Chain
|
||||
|
||||
OperatorResult = namedtuple('OperatorResult', ['k', 'rates'])
|
||||
OperatorResult.__doc__ = """\
|
||||
Result of applying transport operator
|
||||
|
||||
Parameters
|
||||
----------
|
||||
k : float
|
||||
Resulting eigenvalue
|
||||
rates : openmc.deplete.ReactionRates
|
||||
Resulting reaction rates
|
||||
|
||||
"""
|
||||
try:
|
||||
OperatorResult.k.__doc__ = None
|
||||
OperatorResult.rates.__doc__ = None
|
||||
except AttributeError:
|
||||
# Can't set __doc__ on properties on Python 3.4
|
||||
pass
|
||||
|
||||
|
||||
class TransportOperator(metaclass=ABCMeta):
|
||||
"""Abstract class defining a transport operator
|
||||
|
||||
Each depletion integrator is written to work with a generic transport
|
||||
operator that takes a vector of material compositions and returns an
|
||||
eigenvalue and reaction rates. This abstract class sets the requirements for
|
||||
such a transport operator. Users should instantiate
|
||||
:class:`openmc.deplete.Operator` rather than this class.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
chain_file : str, optional
|
||||
Path to the depletion chain XML file. Defaults to the
|
||||
:envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
dilute_initial : float
|
||||
Initial atom density to add for nuclides that are zero in initial
|
||||
condition to ensure they exist in the decay chain. Only done for
|
||||
nuclides with reaction rates. Defaults to 1.0e3.
|
||||
|
||||
"""
|
||||
def __init__(self, chain_file=None):
|
||||
self.dilute_initial = 1.0e3
|
||||
self.output_dir = '.'
|
||||
|
||||
# Read depletion chain
|
||||
if chain_file is None:
|
||||
chain_file = os.environ.get("OPENMC_DEPLETE_CHAIN", None)
|
||||
if chain_file is None:
|
||||
raise IOError("No chain specified, either manually or in "
|
||||
"environment variable OPENMC_DEPLETE_CHAIN.")
|
||||
self.chain = Chain.from_xml(chain_file)
|
||||
|
||||
@abstractmethod
|
||||
def __call__(self, vec, print_out=True):
|
||||
"""Runs a simulation.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
vec : list of numpy.ndarray
|
||||
Total atoms to be used in function.
|
||||
print_out : bool, optional
|
||||
Whether or not to print out time.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.deplete.OperatorResult
|
||||
Eigenvalue and reaction rates resulting from transport operator
|
||||
|
||||
"""
|
||||
pass
|
||||
|
||||
def __enter__(self):
|
||||
# Save current directory and move to specific output directory
|
||||
self._orig_dir = os.getcwd()
|
||||
if not self.output_dir.exists():
|
||||
self.output_dir.mkdir() # exist_ok parameter is 3.5+
|
||||
|
||||
# In Python 3.6+, chdir accepts a Path directly
|
||||
os.chdir(str(self.output_dir))
|
||||
|
||||
return self.initial_condition()
|
||||
|
||||
def __exit__(self, exc_type, exc_value, traceback):
|
||||
self.finalize()
|
||||
os.chdir(self._orig_dir)
|
||||
|
||||
@property
|
||||
def output_dir(self):
|
||||
return self._output_dir
|
||||
|
||||
@output_dir.setter
|
||||
def output_dir(self, output_dir):
|
||||
self._output_dir = Path(output_dir)
|
||||
|
||||
@abstractmethod
|
||||
def initial_condition(self):
|
||||
"""Performs final setup and returns initial condition.
|
||||
|
||||
Returns
|
||||
-------
|
||||
list of numpy.ndarray
|
||||
Total density for initial conditions.
|
||||
"""
|
||||
|
||||
pass
|
||||
|
||||
@abstractmethod
|
||||
def get_results_info(self):
|
||||
"""Returns volume list, cell lists, and nuc lists.
|
||||
|
||||
Returns
|
||||
-------
|
||||
volume : list of float
|
||||
Volumes corresponding to materials in burn_list
|
||||
nuc_list : list of str
|
||||
A list of all nuclide names. Used for sorting the simulation.
|
||||
burn_list : list of int
|
||||
A list of all cell IDs to be burned. Used for sorting the simulation.
|
||||
full_burn_list : list of int
|
||||
All burnable materials in the geometry.
|
||||
"""
|
||||
|
||||
pass
|
||||
|
||||
def finalize(self):
|
||||
pass
|
||||
214
openmc/deplete/atom_number.py
Normal file
214
openmc/deplete/atom_number.py
Normal file
|
|
@ -0,0 +1,214 @@
|
|||
"""AtomNumber module.
|
||||
|
||||
An ndarray to store atom densities with string, integer, or slice indexing.
|
||||
"""
|
||||
from collections import OrderedDict
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
class AtomNumber(object):
|
||||
"""Stores local material compositions (atoms of each nuclide).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
local_mats : list of str
|
||||
Material IDs
|
||||
nuclides : list of str
|
||||
Nuclides to be tracked
|
||||
volume : dict
|
||||
Volume of each material in [cm^3]
|
||||
n_nuc_burn : int
|
||||
Number of nuclides to be burned.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
index_mat : dict
|
||||
A dictionary mapping material ID as string to index.
|
||||
index_nuc : dict
|
||||
A dictionary mapping nuclide name to index.
|
||||
volume : numpy.ndarray
|
||||
Volume of each material in [cm^3]. If a volume is not found, it defaults
|
||||
to 1 so that reading density still works correctly.
|
||||
number : numpy.ndarray
|
||||
Array storing total atoms for each material/nuclide
|
||||
materials : list of str
|
||||
Material IDs as strings
|
||||
nuclides : list of str
|
||||
All nuclide names
|
||||
burnable_nuclides : list of str
|
||||
Burnable nuclides names. Used for sorting the simulation.
|
||||
n_nuc_burn : int
|
||||
Number of burnable nuclides.
|
||||
n_nuc : int
|
||||
Number of nuclides.
|
||||
|
||||
"""
|
||||
def __init__(self, local_mats, nuclides, volume, n_nuc_burn):
|
||||
self.index_mat = OrderedDict((mat, i) for i, mat in enumerate(local_mats))
|
||||
self.index_nuc = OrderedDict((nuc, i) for i, nuc in enumerate(nuclides))
|
||||
|
||||
self.volume = np.ones(len(local_mats))
|
||||
for mat, val in volume.items():
|
||||
if mat in self.index_mat:
|
||||
ind = self.index_mat[mat]
|
||||
self.volume[ind] = val
|
||||
|
||||
self.n_nuc_burn = n_nuc_burn
|
||||
|
||||
self.number = np.zeros((len(local_mats), len(nuclides)))
|
||||
|
||||
def __getitem__(self, pos):
|
||||
"""Retrieves total atom number from AtomNumber.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
pos : tuple
|
||||
A two-length tuple containing a material index and a nuc index.
|
||||
These indexes can be strings (which get converted to integers via
|
||||
the dictionaries), integers used directly, or slices.
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.ndarray
|
||||
The value indexed from self.number.
|
||||
"""
|
||||
|
||||
mat, nuc = pos
|
||||
if isinstance(mat, str):
|
||||
mat = self.index_mat[mat]
|
||||
if isinstance(nuc, str):
|
||||
nuc = self.index_nuc[nuc]
|
||||
|
||||
return self.number[mat, nuc]
|
||||
|
||||
def __setitem__(self, pos, val):
|
||||
"""Sets total atom number into AtomNumber.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
pos : tuple
|
||||
A two-length tuple containing a material index and a nuc index.
|
||||
These indexes can be strings (which get converted to integers via
|
||||
the dictionaries), integers used directly, or slices.
|
||||
val : float
|
||||
The value [atom] to set the array to.
|
||||
|
||||
"""
|
||||
mat, nuc = pos
|
||||
if isinstance(mat, str):
|
||||
mat = self.index_mat[mat]
|
||||
if isinstance(nuc, str):
|
||||
nuc = self.index_nuc[nuc]
|
||||
|
||||
self.number[mat, nuc] = val
|
||||
|
||||
@property
|
||||
def materials(self):
|
||||
return self.index_mat.keys()
|
||||
|
||||
@property
|
||||
def nuclides(self):
|
||||
return self.index_nuc.keys()
|
||||
|
||||
@property
|
||||
def n_nuc(self):
|
||||
return len(self.index_nuc)
|
||||
|
||||
@property
|
||||
def burnable_nuclides(self):
|
||||
return [nuc for nuc, ind in self.index_nuc.items()
|
||||
if ind < self.n_nuc_burn]
|
||||
|
||||
def get_atom_density(self, mat, nuc):
|
||||
"""Accesses atom density instead of total number.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mat : str, int or slice
|
||||
Material index.
|
||||
nuc : str, int or slice
|
||||
Nuclide index.
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.ndarray
|
||||
Density in [atom/cm^3]
|
||||
|
||||
"""
|
||||
if isinstance(mat, str):
|
||||
mat = self.index_mat[mat]
|
||||
if isinstance(nuc, str):
|
||||
nuc = self.index_nuc[nuc]
|
||||
|
||||
return self[mat, nuc] / self.volume[mat]
|
||||
|
||||
def set_atom_density(self, mat, nuc, val):
|
||||
"""Sets atom density instead of total number.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mat : str, int or slice
|
||||
Material index.
|
||||
nuc : str, int or slice
|
||||
Nuclide index.
|
||||
val : numpy.ndarray
|
||||
Array of densities to set in [atom/cm^3]
|
||||
|
||||
"""
|
||||
if isinstance(mat, str):
|
||||
mat = self.index_mat[mat]
|
||||
if isinstance(nuc, str):
|
||||
nuc = self.index_nuc[nuc]
|
||||
|
||||
self[mat, nuc] = val * self.volume[mat]
|
||||
|
||||
def get_mat_slice(self, mat):
|
||||
"""Gets atom quantity indexed by mats for all burned nuclides
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mat : str, int or slice
|
||||
Material index.
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.ndarray
|
||||
The slice requested in [atom].
|
||||
|
||||
"""
|
||||
if isinstance(mat, str):
|
||||
mat = self.index_mat[mat]
|
||||
|
||||
return self[mat, :self.n_nuc_burn]
|
||||
|
||||
def set_mat_slice(self, mat, val):
|
||||
"""Sets atom quantity indexed by mats for all burned nuclides
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mat : str, int or slice
|
||||
Material index.
|
||||
val : numpy.ndarray
|
||||
The slice to set in [atom]
|
||||
|
||||
"""
|
||||
if isinstance(mat, str):
|
||||
mat = self.index_mat[mat]
|
||||
|
||||
self[mat, :self.n_nuc_burn] = val
|
||||
|
||||
def set_density(self, total_density):
|
||||
"""Sets density.
|
||||
|
||||
Sets the density in the exact same order as total_density_list outputs,
|
||||
allowing for internal consistency
|
||||
|
||||
Parameters
|
||||
----------
|
||||
total_density : list of numpy.ndarray
|
||||
Total atoms.
|
||||
|
||||
"""
|
||||
for i, density_slice in enumerate(total_density):
|
||||
self.set_mat_slice(i, density_slice)
|
||||
440
openmc/deplete/chain.py
Normal file
440
openmc/deplete/chain.py
Normal file
|
|
@ -0,0 +1,440 @@
|
|||
"""chain module.
|
||||
|
||||
This module contains information about a depletion chain. A depletion chain is
|
||||
loaded from an .xml file and all the nuclides are linked together.
|
||||
"""
|
||||
|
||||
from collections import OrderedDict, defaultdict
|
||||
from io import StringIO
|
||||
from itertools import chain
|
||||
import math
|
||||
import re
|
||||
import os
|
||||
|
||||
# Try to use lxml if it is available. It preserves the order of attributes and
|
||||
# provides a pretty-printer by default. If not available, use OpenMC function to
|
||||
# pretty print.
|
||||
try:
|
||||
import lxml.etree as ET
|
||||
_have_lxml = True
|
||||
except ImportError:
|
||||
import xml.etree.ElementTree as ET
|
||||
_have_lxml = False
|
||||
import scipy.sparse as sp
|
||||
|
||||
import openmc.data
|
||||
from openmc.clean_xml import clean_xml_indentation
|
||||
from .nuclide import Nuclide, DecayTuple, ReactionTuple
|
||||
|
||||
|
||||
# tuple of (reaction name, possible MT values, (dA, dZ)) where dA is the change
|
||||
# in the mass number and dZ is the change in the atomic number
|
||||
_REACTIONS = [
|
||||
('(n,2n)', set(chain([16], range(875, 892))), (-1, 0)),
|
||||
('(n,3n)', {17}, (-2, 0)),
|
||||
('(n,4n)', {37}, (-3, 0)),
|
||||
('(n,gamma)', {102}, (1, 0)),
|
||||
('(n,p)', set(chain([103], range(600, 650))), (0, -1)),
|
||||
('(n,a)', set(chain([107], range(800, 850))), (-3, -2))
|
||||
]
|
||||
|
||||
|
||||
def replace_missing(product, decay_data):
|
||||
"""Replace missing product with suitable decay daughter.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
product : str
|
||||
Name of product in GND format, e.g. 'Y86_m1'.
|
||||
decay_data : dict
|
||||
Dictionary of decay data
|
||||
|
||||
Returns
|
||||
-------
|
||||
product : str
|
||||
Replacement for missing product in GND format.
|
||||
|
||||
"""
|
||||
# Determine atomic number, mass number, and metastable state
|
||||
Z, A, state = openmc.data.zam(product)
|
||||
symbol = openmc.data.ATOMIC_SYMBOL[Z]
|
||||
|
||||
# Replace neutron with proton
|
||||
if Z == 0 and A == 1:
|
||||
return 'H1'
|
||||
|
||||
# First check if ground state is available
|
||||
if state:
|
||||
product = '{}{}'.format(symbol, A)
|
||||
|
||||
# Find isotope with longest half-life
|
||||
half_life = 0.0
|
||||
for nuclide, data in decay_data.items():
|
||||
m = re.match(r'{}(\d+)(?:_m\d+)?'.format(symbol), nuclide)
|
||||
if m:
|
||||
# If we find a stable nuclide, stop search
|
||||
if data.nuclide['stable']:
|
||||
mass_longest_lived = int(m.group(1))
|
||||
break
|
||||
if data.half_life.nominal_value > half_life:
|
||||
mass_longest_lived = int(m.group(1))
|
||||
half_life = data.half_life.nominal_value
|
||||
|
||||
# If mass number of longest-lived isotope is less than that of missing
|
||||
# product, assume it undergoes beta-. Otherwise assume beta+.
|
||||
beta_minus = (mass_longest_lived < A)
|
||||
|
||||
# Iterate until we find an existing nuclide
|
||||
while product not in decay_data:
|
||||
if Z > 98:
|
||||
Z -= 2
|
||||
A -= 4
|
||||
else:
|
||||
if beta_minus:
|
||||
Z += 1
|
||||
else:
|
||||
Z -= 1
|
||||
product = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
|
||||
|
||||
return product
|
||||
|
||||
|
||||
class Chain(object):
|
||||
"""Full representation of a depletion chain.
|
||||
|
||||
A depletion chain can be created by using the :meth:`from_endf` method which
|
||||
requires a list of ENDF incident neutron, decay, and neutron fission product
|
||||
yield sublibrary files. The depletion chain used during a depletion
|
||||
simulation is indicated by either an argument to
|
||||
:class:`openmc.deplete.Operator` or through the
|
||||
:envvar:`OPENMC_DEPLETE_CHAIN` environment variable.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
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.
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self):
|
||||
self.nuclides = []
|
||||
self.reactions = []
|
||||
self.nuclide_dict = OrderedDict()
|
||||
|
||||
def __contains__(self, nuclide):
|
||||
return nuclide in self.nuclide_dict
|
||||
|
||||
def __getitem__(self, name):
|
||||
"""Get a Nuclide by name."""
|
||||
return self.nuclides[self.nuclide_dict[name]]
|
||||
|
||||
def __len__(self):
|
||||
"""Number of nuclides in chain."""
|
||||
return len(self.nuclides)
|
||||
|
||||
@classmethod
|
||||
def from_endf(cls, decay_files, fpy_files, neutron_files):
|
||||
"""Create a depletion chain from ENDF files.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
decay_files : list of str
|
||||
List of ENDF decay sub-library files
|
||||
fpy_files : list of str
|
||||
List of ENDF neutron-induced fission product yield sub-library files
|
||||
neutron_files : list of str
|
||||
List of ENDF neutron reaction sub-library files
|
||||
|
||||
"""
|
||||
chain = cls()
|
||||
|
||||
# Create dictionary mapping target to filename
|
||||
print('Processing neutron sub-library files...')
|
||||
reactions = {}
|
||||
for f in neutron_files:
|
||||
evaluation = openmc.data.endf.Evaluation(f)
|
||||
name = evaluation.gnd_name
|
||||
reactions[name] = {}
|
||||
for mf, mt, nc, mod in evaluation.reaction_list:
|
||||
if mf == 3:
|
||||
file_obj = StringIO(evaluation.section[3, mt])
|
||||
openmc.data.endf.get_head_record(file_obj)
|
||||
q_value = openmc.data.endf.get_cont_record(file_obj)[1]
|
||||
reactions[name][mt] = q_value
|
||||
|
||||
# Determine what decay and FPY nuclides are available
|
||||
print('Processing decay sub-library files...')
|
||||
decay_data = {}
|
||||
for f in decay_files:
|
||||
data = openmc.data.Decay(f)
|
||||
# Skip decay data for neutron itself
|
||||
if data.nuclide['atomic_number'] == 0:
|
||||
continue
|
||||
decay_data[data.nuclide['name']] = data
|
||||
|
||||
print('Processing fission product yield sub-library files...')
|
||||
fpy_data = {}
|
||||
for f in fpy_files:
|
||||
data = openmc.data.FissionProductYields(f)
|
||||
fpy_data[data.nuclide['name']] = data
|
||||
|
||||
print('Creating depletion_chain...')
|
||||
missing_daughter = []
|
||||
missing_rx_product = []
|
||||
missing_fpy = []
|
||||
missing_fp = []
|
||||
|
||||
for idx, parent in enumerate(sorted(decay_data, key=openmc.data.zam)):
|
||||
data = decay_data[parent]
|
||||
|
||||
nuclide = Nuclide()
|
||||
nuclide.name = parent
|
||||
|
||||
chain.nuclides.append(nuclide)
|
||||
chain.nuclide_dict[parent] = idx
|
||||
|
||||
if not data.nuclide['stable'] and data.half_life.nominal_value != 0.0:
|
||||
nuclide.half_life = data.half_life.nominal_value
|
||||
nuclide.decay_energy = sum(E.nominal_value for E in
|
||||
data.average_energies.values())
|
||||
sum_br = 0.0
|
||||
for i, mode in enumerate(data.modes):
|
||||
type_ = ','.join(mode.modes)
|
||||
if mode.daughter in decay_data:
|
||||
target = mode.daughter
|
||||
else:
|
||||
print('missing {} {} {}'.format(parent, ','.join(mode.modes), mode.daughter))
|
||||
target = replace_missing(mode.daughter, decay_data)
|
||||
|
||||
# Write branching ratio, taking care to ensure sum is unity
|
||||
br = mode.branching_ratio.nominal_value
|
||||
sum_br += br
|
||||
if i == len(data.modes) - 1 and sum_br != 1.0:
|
||||
br = 1.0 - sum(m.branching_ratio.nominal_value
|
||||
for m in data.modes[:-1])
|
||||
|
||||
# Append decay mode
|
||||
nuclide.decay_modes.append(DecayTuple(type_, target, br))
|
||||
|
||||
if parent in reactions:
|
||||
reactions_available = set(reactions[parent].keys())
|
||||
for name, mts, changes in _REACTIONS:
|
||||
if mts & reactions_available:
|
||||
delta_A, delta_Z = changes
|
||||
A = data.nuclide['mass_number'] + delta_A
|
||||
Z = data.nuclide['atomic_number'] + delta_Z
|
||||
daughter = '{}{}'.format(openmc.data.ATOMIC_SYMBOL[Z], A)
|
||||
|
||||
if name not in chain.reactions:
|
||||
chain.reactions.append(name)
|
||||
|
||||
if daughter not in decay_data:
|
||||
missing_rx_product.append((parent, name, daughter))
|
||||
|
||||
# Store Q value
|
||||
for mt in sorted(mts):
|
||||
if mt in reactions[parent]:
|
||||
q_value = reactions[parent][mt]
|
||||
break
|
||||
else:
|
||||
q_value = 0.0
|
||||
|
||||
nuclide.reactions.append(ReactionTuple(
|
||||
name, daughter, q_value, 1.0))
|
||||
|
||||
if any(mt in reactions_available for mt in [18, 19, 20, 21, 38]):
|
||||
if parent in fpy_data:
|
||||
q_value = reactions[parent][18]
|
||||
nuclide.reactions.append(
|
||||
ReactionTuple('fission', 0, q_value, 1.0))
|
||||
|
||||
if 'fission' not in chain.reactions:
|
||||
chain.reactions.append('fission')
|
||||
else:
|
||||
missing_fpy.append(parent)
|
||||
|
||||
if parent in fpy_data:
|
||||
fpy = fpy_data[parent]
|
||||
|
||||
if fpy.energies is not None:
|
||||
nuclide.yield_energies = fpy.energies
|
||||
else:
|
||||
nuclide.yield_energies = [0.0]
|
||||
|
||||
for E, table in zip(nuclide.yield_energies, fpy.independent):
|
||||
yield_replace = 0.0
|
||||
yields = defaultdict(float)
|
||||
for product, y in table.items():
|
||||
# Handle fission products that have no decay data available
|
||||
if product not in decay_data:
|
||||
daughter = replace_missing(product, decay_data)
|
||||
product = daughter
|
||||
yield_replace += y.nominal_value
|
||||
|
||||
yields[product] += y.nominal_value
|
||||
|
||||
if yield_replace > 0.0:
|
||||
missing_fp.append((parent, E, yield_replace))
|
||||
|
||||
nuclide.yield_data[E] = []
|
||||
for k in sorted(yields, key=openmc.data.zam):
|
||||
nuclide.yield_data[E].append((k, yields[k]))
|
||||
|
||||
# Display warnings
|
||||
if missing_daughter:
|
||||
print('The following decay modes have daughters with no decay data:')
|
||||
for mode in missing_daughter:
|
||||
print(' {}'.format(mode))
|
||||
print('')
|
||||
|
||||
if missing_rx_product:
|
||||
print('The following reaction products have no decay data:')
|
||||
for vals in missing_rx_product:
|
||||
print('{} {} -> {}'.format(*vals))
|
||||
print('')
|
||||
|
||||
if missing_fpy:
|
||||
print('The following fissionable nuclides have no fission product yields:')
|
||||
for parent in missing_fpy:
|
||||
print(' ' + parent)
|
||||
print('')
|
||||
|
||||
if missing_fp:
|
||||
print('The following nuclides have fission products with no decay data:')
|
||||
for vals in missing_fp:
|
||||
print(' {}, E={} eV (total yield={})'.format(*vals))
|
||||
|
||||
return chain
|
||||
|
||||
@classmethod
|
||||
def from_xml(cls, filename):
|
||||
"""Reads a depletion chain XML file.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str
|
||||
The path to the depletion chain XML file.
|
||||
|
||||
"""
|
||||
chain = cls()
|
||||
|
||||
# Load XML tree
|
||||
root = ET.parse(str(filename))
|
||||
|
||||
for i, nuclide_elem in enumerate(root.findall('nuclide')):
|
||||
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.reactions:
|
||||
chain.reactions.append(rx.type)
|
||||
|
||||
chain.nuclides.append(nuc)
|
||||
|
||||
return chain
|
||||
|
||||
def export_to_xml(self, filename):
|
||||
"""Writes a depletion chain XML file.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str
|
||||
The path to the depletion chain XML file.
|
||||
|
||||
"""
|
||||
|
||||
root_elem = ET.Element('depletion_chain')
|
||||
for nuclide in self.nuclides:
|
||||
root_elem.append(nuclide.to_xml_element())
|
||||
|
||||
tree = ET.ElementTree(root_elem)
|
||||
if _have_lxml:
|
||||
tree.write(str(filename), encoding='utf-8', pretty_print=True)
|
||||
else:
|
||||
clean_xml_indentation(root_elem)
|
||||
tree.write(str(filename), encoding='utf-8')
|
||||
|
||||
def form_matrix(self, rates):
|
||||
"""Forms depletion matrix.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
rates : numpy.ndarray
|
||||
2D array indexed by (nuclide, reaction)
|
||||
|
||||
Returns
|
||||
-------
|
||||
scipy.sparse.csr_matrix
|
||||
Sparse matrix representing depletion.
|
||||
|
||||
"""
|
||||
matrix = defaultdict(float)
|
||||
reactions = set()
|
||||
|
||||
for i, nuc in enumerate(self.nuclides):
|
||||
|
||||
if nuc.n_decay_modes != 0:
|
||||
# Decay paths
|
||||
# Loss
|
||||
decay_constant = math.log(2) / nuc.half_life
|
||||
|
||||
if decay_constant != 0.0:
|
||||
matrix[i, i] -= decay_constant
|
||||
|
||||
# Gain
|
||||
for _, target, branching_ratio in nuc.decay_modes:
|
||||
# Allow for total annihilation for debug purposes
|
||||
if target != 'Nothing':
|
||||
branch_val = branching_ratio * decay_constant
|
||||
|
||||
if branch_val != 0.0:
|
||||
k = self.nuclide_dict[target]
|
||||
matrix[k, i] += branch_val
|
||||
|
||||
if nuc.name in rates.index_nuc:
|
||||
# Extract all reactions for this nuclide in this cell
|
||||
nuc_ind = rates.index_nuc[nuc.name]
|
||||
nuc_rates = rates[nuc_ind, :]
|
||||
|
||||
for r_type, target, _, br in nuc.reactions:
|
||||
# Extract reaction index, and then final reaction rate
|
||||
r_id = rates.index_rx[r_type]
|
||||
path_rate = nuc_rates[r_id]
|
||||
|
||||
# Loss term -- make sure we only count loss once for
|
||||
# reactions with branching ratios
|
||||
if r_type not in reactions:
|
||||
reactions.add(r_type)
|
||||
if path_rate != 0.0:
|
||||
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:
|
||||
# Assume that we should always use thermal fission
|
||||
# yields. At some point it would be nice to account
|
||||
# for the energy-dependence..
|
||||
energy, data = sorted(nuc.yield_data.items())[0]
|
||||
for product, y in data:
|
||||
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()
|
||||
|
||||
# Use DOK matrix as intermediate representation, then convert to CSR and return
|
||||
n = len(self)
|
||||
matrix_dok = sp.dok_matrix((n, n))
|
||||
dict.update(matrix_dok, matrix)
|
||||
return matrix_dok.tocsr()
|
||||
27
openmc/deplete/dummy_comm.py
Normal file
27
openmc/deplete/dummy_comm.py
Normal file
|
|
@ -0,0 +1,27 @@
|
|||
class DummyCommunicator(object):
|
||||
rank = 0
|
||||
size = 1
|
||||
|
||||
def allgather(self, sendobj):
|
||||
return [sendobj]
|
||||
|
||||
def allreduce(self, sendobj, op=None):
|
||||
return sendobj
|
||||
|
||||
def barrier(self):
|
||||
pass
|
||||
|
||||
def bcast(self, obj, root=0):
|
||||
return obj
|
||||
|
||||
def gather(self, sendobj, root=0):
|
||||
return [sendobj]
|
||||
|
||||
def py2f(self):
|
||||
return 0
|
||||
|
||||
def reduce(self, sendobj, op=None, root=0):
|
||||
return sendobj
|
||||
|
||||
def scatter(self, sendobj, root=0):
|
||||
return sendobj[0]
|
||||
10
openmc/deplete/integrator/__init__.py
Normal file
10
openmc/deplete/integrator/__init__.py
Normal file
|
|
@ -0,0 +1,10 @@
|
|||
"""
|
||||
Integrator
|
||||
===========
|
||||
|
||||
The integrator subcomponents.
|
||||
"""
|
||||
|
||||
from .cecm import *
|
||||
from .cram import *
|
||||
from .predictor import *
|
||||
78
openmc/deplete/integrator/cecm.py
Normal file
78
openmc/deplete/integrator/cecm.py
Normal file
|
|
@ -0,0 +1,78 @@
|
|||
"""The CE/CM integrator."""
|
||||
|
||||
import copy
|
||||
from collections.abc import Iterable
|
||||
|
||||
from .cram import deplete
|
||||
from ..results import Results
|
||||
|
||||
|
||||
def cecm(operator, timesteps, power, print_out=True):
|
||||
r"""Deplete using the CE/CM algorithm.
|
||||
|
||||
Implements the second order `CE/CM predictor-corrector algorithm
|
||||
<https://doi.org/10.13182/NSE14-92>`_. This algorithm is mathematically
|
||||
defined as:
|
||||
|
||||
.. math::
|
||||
y' &= A(y, t) y(t)
|
||||
|
||||
A_p &= A(y_n, t_n)
|
||||
|
||||
y_m &= \text{expm}(A_p h/2) y_n
|
||||
|
||||
A_c &= A(y_m, t_n + h/2)
|
||||
|
||||
y_{n+1} &= \text{expm}(A_c h) y_n
|
||||
|
||||
Parameters
|
||||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
The operator object to simulate on.
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not cumulative.
|
||||
power : float or iterable of float
|
||||
Power of the reactor in [W]. A single value indicates that the power is
|
||||
constant over all timesteps. An iterable indicates potentially different
|
||||
power levels for each timestep. For a 2D problem, the power can be given
|
||||
in [W/cm] as long as the "volume" assigned to a depletion material is
|
||||
actually an area in [cm^2].
|
||||
print_out : bool, optional
|
||||
Whether or not to print out time.
|
||||
|
||||
"""
|
||||
if not isinstance(power, Iterable):
|
||||
power = [power]*len(timesteps)
|
||||
|
||||
# Generate initial conditions
|
||||
with operator as vec:
|
||||
chain = operator.chain
|
||||
t = 0.0
|
||||
for i, (dt, p) in enumerate(zip(timesteps, power)):
|
||||
# Get beginning-of-timestep reaction rates
|
||||
x = [copy.deepcopy(vec)]
|
||||
op_results = [operator(x[0], p)]
|
||||
|
||||
# Deplete for first half of timestep
|
||||
x_middle = deplete(chain, x[0], op_results[0], dt/2, print_out)
|
||||
|
||||
# Get middle-of-timestep reaction rates
|
||||
x.append(x_middle)
|
||||
op_results.append(operator(x_middle, p))
|
||||
|
||||
# Deplete for full timestep using beginning-of-step materials
|
||||
x_end = deplete(chain, x[0], op_results[1], dt, print_out)
|
||||
|
||||
# Create results, write to disk
|
||||
Results.save(operator, x, op_results, [t, t + dt], i)
|
||||
|
||||
# Advance time, update vector
|
||||
t += dt
|
||||
vec = copy.deepcopy(x_end)
|
||||
|
||||
# Perform one last simulation
|
||||
x = [copy.deepcopy(vec)]
|
||||
op_results = [operator(x[0], power[-1])]
|
||||
|
||||
# Create results, write to disk
|
||||
Results.save(operator, x, op_results, [t, t], len(timesteps))
|
||||
227
openmc/deplete/integrator/cram.py
Normal file
227
openmc/deplete/integrator/cram.py
Normal file
|
|
@ -0,0 +1,227 @@
|
|||
"""Chebyshev Rational Approximation Method module
|
||||
|
||||
Implements two different forms of CRAM for use in openmc.deplete.
|
||||
"""
|
||||
|
||||
from itertools import repeat
|
||||
from multiprocessing import Pool
|
||||
import time
|
||||
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
import scipy.sparse.linalg as sla
|
||||
|
||||
from .. import comm
|
||||
|
||||
|
||||
def deplete(chain, x, op_result, dt, print_out):
|
||||
"""Deplete materials using given reaction rates for a specified time
|
||||
|
||||
Parameters
|
||||
----------
|
||||
chain : openmc.deplete.Chain
|
||||
Depletion chain
|
||||
x : list of numpy.ndarray
|
||||
Atom number vectors for each material
|
||||
op_result : openmc.deplete.OperatorResult
|
||||
Result of applying transport operator (contains reaction rates)
|
||||
dt : float
|
||||
Time in [s] to deplete for
|
||||
print_out : bool
|
||||
Whether to show elapsed time
|
||||
|
||||
Returns
|
||||
-------
|
||||
x_result : list of numpy.ndarray
|
||||
Updated atom number vectors for each material
|
||||
|
||||
"""
|
||||
t_start = time.time()
|
||||
|
||||
# Set up iterators
|
||||
n_mats = len(x)
|
||||
chains = repeat(chain, n_mats)
|
||||
vecs = (x[i] for i in range(n_mats))
|
||||
rates = (op_result.rates[i, :, :] for i in range(n_mats))
|
||||
dts = repeat(dt, n_mats)
|
||||
|
||||
# Use multiprocessing pool to distribute work
|
||||
with Pool() as pool:
|
||||
iters = zip(chains, vecs, rates, dts)
|
||||
x_result = list(pool.starmap(_cram_wrapper, iters))
|
||||
|
||||
t_end = time.time()
|
||||
if comm.rank == 0:
|
||||
if print_out:
|
||||
print("Time to matexp: ", t_end - t_start)
|
||||
|
||||
return x_result
|
||||
|
||||
|
||||
def _cram_wrapper(chain, n0, rates, dt):
|
||||
"""Wraps depletion matrix creation / CRAM solve for multiprocess execution
|
||||
|
||||
Parameters
|
||||
----------
|
||||
chain : DepletionChain
|
||||
Depletion chain used to construct the burnup matrix
|
||||
n0 : numpy.array
|
||||
Vector to operate a matrix exponent on.
|
||||
rates : numpy.ndarray
|
||||
2D array indexed by nuclide then by cell.
|
||||
dt : float
|
||||
Time to integrate to.
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.array
|
||||
Results of the matrix exponent.
|
||||
"""
|
||||
A = chain.form_matrix(rates)
|
||||
return CRAM48(A, n0, dt)
|
||||
|
||||
|
||||
def CRAM16(A, n0, dt):
|
||||
"""Chebyshev Rational Approximation Method, order 16
|
||||
|
||||
Algorithm is the 16th order Chebyshev Rational Approximation Method,
|
||||
implemented in the more stable `incomplete partial fraction (IPF)
|
||||
<https://doi.org/10.13182/NSE15-26>`_ form.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
A : scipy.linalg.csr_matrix
|
||||
Matrix to take exponent of.
|
||||
n0 : numpy.array
|
||||
Vector to operate a matrix exponent on.
|
||||
dt : float
|
||||
Time to integrate to.
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.array
|
||||
Results of the matrix exponent.
|
||||
|
||||
"""
|
||||
|
||||
alpha = np.array([+2.124853710495224e-16,
|
||||
+5.464930576870210e+3 - 3.797983575308356e+4j,
|
||||
+9.045112476907548e+1 - 1.115537522430261e+3j,
|
||||
+2.344818070467641e+2 - 4.228020157070496e+2j,
|
||||
+9.453304067358312e+1 - 2.951294291446048e+2j,
|
||||
+7.283792954673409e+2 - 1.205646080220011e+5j,
|
||||
+3.648229059594851e+1 - 1.155509621409682e+2j,
|
||||
+2.547321630156819e+1 - 2.639500283021502e+1j,
|
||||
+2.394538338734709e+1 - 5.650522971778156e+0j],
|
||||
dtype=np.complex128)
|
||||
theta = np.array([+0.0,
|
||||
+3.509103608414918 + 8.436198985884374j,
|
||||
+5.948152268951177 + 3.587457362018322j,
|
||||
-5.264971343442647 + 16.22022147316793j,
|
||||
+1.419375897185666 + 10.92536348449672j,
|
||||
+6.416177699099435 + 1.194122393370139j,
|
||||
+4.993174737717997 + 5.996881713603942j,
|
||||
-1.413928462488886 + 13.49772569889275j,
|
||||
-10.84391707869699 + 19.27744616718165j],
|
||||
dtype=np.complex128)
|
||||
|
||||
n = A.shape[0]
|
||||
|
||||
alpha0 = 2.124853710495224e-16
|
||||
|
||||
k = 8
|
||||
|
||||
y = np.array(n0, dtype=np.float64)
|
||||
for l in range(1, k+1):
|
||||
y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y
|
||||
|
||||
y *= alpha0
|
||||
return y
|
||||
|
||||
|
||||
def CRAM48(A, n0, dt):
|
||||
"""Chebyshev Rational Approximation Method, order 48
|
||||
|
||||
Algorithm is the 48th order Chebyshev Rational Approximation Method,
|
||||
implemented in the more stable `incomplete partial fraction (IPF)
|
||||
<https://doi.org/10.13182/NSE15-26>`_ form.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
A : scipy.linalg.csr_matrix
|
||||
Matrix to take exponent of.
|
||||
n0 : numpy.array
|
||||
Vector to operate a matrix exponent on.
|
||||
dt : float
|
||||
Time to integrate to.
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.array
|
||||
Results of the matrix exponent.
|
||||
|
||||
"""
|
||||
|
||||
theta_r = np.array([-4.465731934165702e+1, -5.284616241568964e+0,
|
||||
-8.867715667624458e+0, +3.493013124279215e+0,
|
||||
+1.564102508858634e+1, +1.742097597385893e+1,
|
||||
-2.834466755180654e+1, +1.661569367939544e+1,
|
||||
+8.011836167974721e+0, -2.056267541998229e+0,
|
||||
+1.449208170441839e+1, +1.853807176907916e+1,
|
||||
+9.932562704505182e+0, -2.244223871767187e+1,
|
||||
+8.590014121680897e-1, -1.286192925744479e+1,
|
||||
+1.164596909542055e+1, +1.806076684783089e+1,
|
||||
+5.870672154659249e+0, -3.542938819659747e+1,
|
||||
+1.901323489060250e+1, +1.885508331552577e+1,
|
||||
-1.734689708174982e+1, +1.316284237125190e+1])
|
||||
theta_i = np.array([+6.233225190695437e+1, +4.057499381311059e+1,
|
||||
+4.325515754166724e+1, +3.281615453173585e+1,
|
||||
+1.558061616372237e+1, +1.076629305714420e+1,
|
||||
+5.492841024648724e+1, +1.316994930024688e+1,
|
||||
+2.780232111309410e+1, +3.794824788914354e+1,
|
||||
+1.799988210051809e+1, +5.974332563100539e+0,
|
||||
+2.532823409972962e+1, +5.179633600312162e+1,
|
||||
+3.536456194294350e+1, +4.600304902833652e+1,
|
||||
+2.287153304140217e+1, +8.368200580099821e+0,
|
||||
+3.029700159040121e+1, +5.834381701800013e+1,
|
||||
+1.194282058271408e+0, +3.583428564427879e+0,
|
||||
+4.883941101108207e+1, +2.042951874827759e+1])
|
||||
theta = np.array(theta_r + theta_i * 1j, dtype=np.complex128)
|
||||
|
||||
alpha_r = np.array([+6.387380733878774e+2, +1.909896179065730e+2,
|
||||
+4.236195226571914e+2, +4.645770595258726e+2,
|
||||
+7.765163276752433e+2, +1.907115136768522e+3,
|
||||
+2.909892685603256e+3, +1.944772206620450e+2,
|
||||
+1.382799786972332e+5, +5.628442079602433e+3,
|
||||
+2.151681283794220e+2, +1.324720240514420e+3,
|
||||
+1.617548476343347e+4, +1.112729040439685e+2,
|
||||
+1.074624783191125e+2, +8.835727765158191e+1,
|
||||
+9.354078136054179e+1, +9.418142823531573e+1,
|
||||
+1.040012390717851e+2, +6.861882624343235e+1,
|
||||
+8.766654491283722e+1, +1.056007619389650e+2,
|
||||
+7.738987569039419e+1, +1.041366366475571e+2])
|
||||
alpha_i = np.array([-6.743912502859256e+2, -3.973203432721332e+2,
|
||||
-2.041233768918671e+3, -1.652917287299683e+3,
|
||||
-1.783617639907328e+4, -5.887068595142284e+4,
|
||||
-9.953255345514560e+3, -1.427131226068449e+3,
|
||||
-3.256885197214938e+6, -2.924284515884309e+4,
|
||||
-1.121774011188224e+3, -6.370088443140973e+4,
|
||||
-1.008798413156542e+6, -8.837109731680418e+1,
|
||||
-1.457246116408180e+2, -6.388286188419360e+1,
|
||||
-2.195424319460237e+2, -6.719055740098035e+2,
|
||||
-1.693747595553868e+2, -1.177598523430493e+1,
|
||||
-4.596464999363902e+3, -1.738294585524067e+3,
|
||||
-4.311715386228984e+1, -2.777743732451969e+2])
|
||||
alpha = np.array(alpha_r + alpha_i * 1j, dtype=np.complex128)
|
||||
n = A.shape[0]
|
||||
|
||||
alpha0 = 2.258038182743983e-47
|
||||
|
||||
k = 24
|
||||
|
||||
y = np.array(n0, dtype=np.float64)
|
||||
for l in range(k):
|
||||
y = 2.0*np.real(alpha[l]*sla.spsolve(A*dt - theta[l]*sp.eye(n), y)) + y
|
||||
|
||||
y *= alpha0
|
||||
return y
|
||||
66
openmc/deplete/integrator/predictor.py
Normal file
66
openmc/deplete/integrator/predictor.py
Normal file
|
|
@ -0,0 +1,66 @@
|
|||
"""First-order predictor algorithm."""
|
||||
|
||||
import copy
|
||||
from collections.abc import Iterable
|
||||
|
||||
from .cram import deplete
|
||||
from ..results import Results
|
||||
|
||||
|
||||
def predictor(operator, timesteps, power, print_out=True):
|
||||
r"""Deplete using a first-order predictor algorithm.
|
||||
|
||||
Implements the first-order predictor algorithm. This algorithm is
|
||||
mathematically defined as:
|
||||
|
||||
.. math::
|
||||
y' &= A(y, t) y(t)
|
||||
|
||||
A_p &= A(y_n, t_n)
|
||||
|
||||
y_{n+1} &= \text{expm}(A_p h) y_n
|
||||
|
||||
Parameters
|
||||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
The operator object to simulate on.
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not cumulative.
|
||||
power : float or iterable of float
|
||||
Power of the reactor in [W]. A single value indicates that the power is
|
||||
constant over all timesteps. An iterable indicates potentially different
|
||||
power levels for each timestep. For a 2D problem, the power can be given
|
||||
in [W/cm] as long as the "volume" assigned to a depletion material is
|
||||
actually an area in [cm^2].
|
||||
print_out : bool, optional
|
||||
Whether or not to print out time.
|
||||
|
||||
"""
|
||||
if not isinstance(power, Iterable):
|
||||
power = [power]*len(timesteps)
|
||||
|
||||
# Generate initial conditions
|
||||
with operator as vec:
|
||||
chain = operator.chain
|
||||
t = 0.0
|
||||
for i, (dt, p) in enumerate(zip(timesteps, power)):
|
||||
# Get beginning-of-timestep reaction rates
|
||||
x = [copy.deepcopy(vec)]
|
||||
op_results = [operator(x[0], p)]
|
||||
|
||||
# Create results, write to disk
|
||||
Results.save(operator, x, op_results, [t, t + dt], i)
|
||||
|
||||
# Deplete for full timestep
|
||||
x_end = deplete(chain, x[0], op_results[0], dt, print_out)
|
||||
|
||||
# Advance time, update vector
|
||||
t += dt
|
||||
vec = copy.deepcopy(x_end)
|
||||
|
||||
# Perform one last simulation
|
||||
x = [copy.deepcopy(vec)]
|
||||
op_results = [operator(x[0], power[-1])]
|
||||
|
||||
# Create results, write to disk
|
||||
Results.save(operator, x, op_results, [t, t], len(timesteps))
|
||||
219
openmc/deplete/nuclide.py
Normal file
219
openmc/deplete/nuclide.py
Normal file
|
|
@ -0,0 +1,219 @@
|
|||
"""Nuclide module.
|
||||
|
||||
Contains the per-nuclide components of a depletion chain.
|
||||
"""
|
||||
|
||||
from collections import namedtuple
|
||||
try:
|
||||
import lxml.etree as ET
|
||||
except ImportError:
|
||||
import xml.etree.ElementTree as ET
|
||||
|
||||
|
||||
DecayTuple = namedtuple('DecayTuple', 'type target branching_ratio')
|
||||
DecayTuple.__doc__ = """\
|
||||
Decay mode information
|
||||
|
||||
Parameters
|
||||
----------
|
||||
type : str
|
||||
Type of the decay mode, e.g., 'beta-'
|
||||
target : str
|
||||
Nuclide resulting from decay
|
||||
branching_ratio : float
|
||||
Branching ratio of the decay mode
|
||||
|
||||
"""
|
||||
try:
|
||||
DecayTuple.type.__doc__ = None
|
||||
DecayTuple.target.__doc__ = None
|
||||
DecayTuple.branching_ratio.__doc__ = None
|
||||
except AttributeError:
|
||||
# Can't set __doc__ on properties on Python 3.4
|
||||
pass
|
||||
|
||||
|
||||
ReactionTuple = namedtuple('ReactionTuple', 'type target Q branching_ratio')
|
||||
ReactionTuple.__doc__ = """\
|
||||
Transmutation reaction information
|
||||
|
||||
Parameters
|
||||
----------
|
||||
type : str
|
||||
Type of the reaction, e.g., 'fission'
|
||||
target : str
|
||||
nuclide resulting from reaction
|
||||
Q : float
|
||||
Q value of the reaction in [eV]
|
||||
branching_ratio : float
|
||||
Branching ratio of the reaction
|
||||
|
||||
"""
|
||||
try:
|
||||
ReactionTuple.type.__doc__ = None
|
||||
ReactionTuple.target.__doc__ = None
|
||||
ReactionTuple.Q.__doc__ = None
|
||||
ReactionTuple.branching_ratio.__doc__ = None
|
||||
except AttributeError:
|
||||
pass
|
||||
|
||||
|
||||
class Nuclide(object):
|
||||
"""Decay modes, reactions, and fission yields for a single nuclide.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
name : str
|
||||
Name of nuclide.
|
||||
half_life : float
|
||||
Half life of nuclide in [s].
|
||||
decay_energy : float
|
||||
Energy deposited from decay in [eV].
|
||||
n_decay_modes : int
|
||||
Number of decay pathways.
|
||||
decay_modes : list of openmc.deplete.DecayTuple
|
||||
Decay mode information. Each element of the list is a named tuple with
|
||||
attributes 'type', 'target', and 'branching_ratio'.
|
||||
n_reaction_paths : int
|
||||
Number of possible reaction pathways.
|
||||
reactions : list of openmc.deplete.ReactionTuple
|
||||
Reaction information. Each element of the list is a named tuple with
|
||||
attribute 'type', 'target', 'Q', and 'branching_ratio'.
|
||||
yield_data : dict of float to list
|
||||
Maps tabulated energy to list of (product, yield) for all
|
||||
neutron-induced fission products.
|
||||
yield_energies : list of float
|
||||
Energies at which fission product yiels exist
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self):
|
||||
# Information about the nuclide
|
||||
self.name = None
|
||||
self.half_life = None
|
||||
self.decay_energy = 0.0
|
||||
|
||||
# Decay paths
|
||||
self.decay_modes = []
|
||||
|
||||
# Reaction paths
|
||||
self.reactions = []
|
||||
|
||||
# Neutron fission yields, if present
|
||||
self.yield_data = {}
|
||||
self.yield_energies = []
|
||||
|
||||
@property
|
||||
def n_decay_modes(self):
|
||||
return len(self.decay_modes)
|
||||
|
||||
@property
|
||||
def n_reaction_paths(self):
|
||||
return len(self.reactions)
|
||||
|
||||
@classmethod
|
||||
def from_xml(cls, element):
|
||||
"""Read nuclide from an XML element.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
element : xml.etree.ElementTree.Element
|
||||
XML element to write nuclide data to
|
||||
|
||||
Returns
|
||||
-------
|
||||
nuc : openmc.deplete.Nuclide
|
||||
Instance of a nuclide
|
||||
|
||||
"""
|
||||
nuc = cls()
|
||||
nuc.name = element.get('name')
|
||||
|
||||
# Check for half-life
|
||||
if 'half_life' in element.attrib:
|
||||
nuc.half_life = float(element.get('half_life'))
|
||||
nuc.decay_energy = float(element.get('decay_energy', '0'))
|
||||
|
||||
# Check for decay paths
|
||||
for decay_elem in element.iter('decay'):
|
||||
d_type = decay_elem.get('type')
|
||||
target = decay_elem.get('target')
|
||||
branching_ratio = float(decay_elem.get('branching_ratio'))
|
||||
nuc.decay_modes.append(DecayTuple(d_type, target, branching_ratio))
|
||||
|
||||
# Check for reaction paths
|
||||
for reaction_elem in element.iter('reaction'):
|
||||
r_type = reaction_elem.get('type')
|
||||
Q = float(reaction_elem.get('Q', '0'))
|
||||
branching_ratio = float(reaction_elem.get('branching_ratio', '1'))
|
||||
|
||||
# If the type is not fission, get target and Q value, otherwise
|
||||
# just set null values
|
||||
if r_type != 'fission':
|
||||
target = reaction_elem.get('target')
|
||||
else:
|
||||
target = None
|
||||
|
||||
# Append reaction
|
||||
nuc.reactions.append(ReactionTuple(
|
||||
r_type, target, Q, branching_ratio))
|
||||
|
||||
fpy_elem = element.find('neutron_fission_yields')
|
||||
if fpy_elem is not None:
|
||||
for yields_elem in fpy_elem.iter('fission_yields'):
|
||||
E = float(yields_elem.get('energy'))
|
||||
products = yields_elem.find('products').text.split()
|
||||
yields = [float(y) for y in
|
||||
yields_elem.find('data').text.split()]
|
||||
nuc.yield_data[E] = list(zip(products, yields))
|
||||
nuc.yield_energies = list(sorted(nuc.yield_data.keys()))
|
||||
|
||||
return nuc
|
||||
|
||||
def to_xml_element(self):
|
||||
"""Write nuclide to XML element.
|
||||
|
||||
Returns
|
||||
-------
|
||||
elem : xml.etree.ElementTree.Element
|
||||
XML element to write nuclide data to
|
||||
|
||||
"""
|
||||
elem = ET.Element('nuclide')
|
||||
elem.set('name', self.name)
|
||||
|
||||
if self.half_life is not None:
|
||||
elem.set('half_life', str(self.half_life))
|
||||
elem.set('decay_modes', str(len(self.decay_modes)))
|
||||
elem.set('decay_energy', str(self.decay_energy))
|
||||
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('branching_ratio', str(br))
|
||||
|
||||
elem.set('reactions', str(len(self.reactions)))
|
||||
for rx, daughter, Q, br in self.reactions:
|
||||
rx_elem = ET.SubElement(elem, 'reaction')
|
||||
rx_elem.set('type', rx)
|
||||
rx_elem.set('Q', str(Q))
|
||||
if rx != 'fission':
|
||||
rx_elem.set('target', daughter)
|
||||
if br != 1.0:
|
||||
rx_elem.set('branching_ratio', str(br))
|
||||
|
||||
if self.yield_data:
|
||||
fpy_elem = ET.SubElement(elem, 'neutron_fission_yields')
|
||||
energy_elem = ET.SubElement(fpy_elem, 'energies')
|
||||
energy_elem.text = ' '.join(str(E) for E in self.yield_energies)
|
||||
|
||||
for E in self.yield_energies:
|
||||
yields_elem = ET.SubElement(fpy_elem, 'fission_yields')
|
||||
yields_elem.set('energy', str(E))
|
||||
|
||||
products_elem = ET.SubElement(yields_elem, 'products')
|
||||
products_elem.text = ' '.join(x[0] for x in self.yield_data[E])
|
||||
data_elem = ET.SubElement(yields_elem, 'data')
|
||||
data_elem.text = ' '.join(str(x[1]) for x in self.yield_data[E])
|
||||
|
||||
return elem
|
||||
562
openmc/deplete/operator.py
Normal file
562
openmc/deplete/operator.py
Normal file
|
|
@ -0,0 +1,562 @@
|
|||
"""OpenMC transport operator
|
||||
|
||||
This module implements a transport operator for OpenMC so that it can be used by
|
||||
depletion integrators. The implementation makes use of the Python bindings to
|
||||
OpenMC's C API so that reading tally results and updating material number
|
||||
densities is all done in-memory instead of through the filesystem.
|
||||
|
||||
"""
|
||||
|
||||
import copy
|
||||
from collections import OrderedDict
|
||||
from itertools import chain
|
||||
import os
|
||||
import time
|
||||
import xml.etree.ElementTree as ET
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
|
||||
import openmc
|
||||
import openmc.capi
|
||||
from openmc.data import JOULE_PER_EV
|
||||
from . import comm
|
||||
from .abc import TransportOperator, OperatorResult
|
||||
from .atom_number import AtomNumber
|
||||
from .reaction_rates import ReactionRates
|
||||
|
||||
|
||||
def _distribute(items):
|
||||
"""Distribute items across MPI communicator
|
||||
|
||||
Parameters
|
||||
----------
|
||||
items : list
|
||||
List of items of distribute
|
||||
|
||||
Returns
|
||||
-------
|
||||
list
|
||||
Items assigned to process that called
|
||||
|
||||
"""
|
||||
min_size, extra = divmod(len(items), comm.size)
|
||||
j = 0
|
||||
for i in range(comm.size):
|
||||
chunk_size = min_size + int(i < extra)
|
||||
if comm.rank == i:
|
||||
return items[j:j + chunk_size]
|
||||
j += chunk_size
|
||||
|
||||
|
||||
class Operator(TransportOperator):
|
||||
"""OpenMC transport operator for depletion.
|
||||
|
||||
Instances of this class can be used to perform depletion using OpenMC as the
|
||||
transport operator. Normally, a user needn't call methods of this class
|
||||
directly. Instead, an instance of this class is passed to an integrator
|
||||
function, such as :func:`openmc.deplete.integrator.cecm`.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
geometry : openmc.Geometry
|
||||
OpenMC geometry object
|
||||
settings : openmc.Settings
|
||||
OpenMC Settings object
|
||||
chain_file : str, optional
|
||||
Path to the depletion chain XML file. Defaults to the
|
||||
:envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
geometry : openmc.Geometry
|
||||
OpenMC geometry object
|
||||
settings : openmc.Settings
|
||||
OpenMC settings object
|
||||
dilute_initial : float
|
||||
Initial atom density to add for nuclides that are zero in initial
|
||||
condition to ensure they exist in the decay chain. Only done for
|
||||
nuclides with reaction rates. Defaults to 1.0e3.
|
||||
output_dir : pathlib.Path
|
||||
Path to output directory to save results.
|
||||
round_number : bool
|
||||
Whether or not to round output to OpenMC to 8 digits.
|
||||
Useful in testing, as OpenMC is incredibly sensitive to exact values.
|
||||
number : openmc.deplete.AtomNumber
|
||||
Total number of atoms in simulation.
|
||||
nuclides_with_data : set of str
|
||||
A set listing all unique nuclides available from cross_sections.xml.
|
||||
chain : openmc.deplete.Chain
|
||||
The depletion chain information necessary to form matrices and tallies.
|
||||
reaction_rates : openmc.deplete.ReactionRates
|
||||
Reaction rates from the last operator step.
|
||||
burnable_mats : list of str
|
||||
All burnable material IDs
|
||||
local_mats : list of str
|
||||
All burnable material IDs being managed by a single process
|
||||
|
||||
"""
|
||||
def __init__(self, geometry, settings, chain_file=None):
|
||||
super().__init__(chain_file)
|
||||
self.round_number = False
|
||||
self.settings = settings
|
||||
self.geometry = geometry
|
||||
|
||||
# Clear out OpenMC, create task lists, distribute
|
||||
openmc.reset_auto_ids()
|
||||
self.burnable_mats, volume, nuclides = self._get_burnable_mats()
|
||||
self.local_mats = _distribute(self.burnable_mats)
|
||||
|
||||
# Determine which nuclides have incident neutron data
|
||||
self.nuclides_with_data = self._get_nuclides_with_data()
|
||||
self._burnable_nucs = [nuc for nuc in self.nuclides_with_data
|
||||
if nuc in self.chain]
|
||||
|
||||
# Extract number densities from the geometry
|
||||
self._extract_number(self.local_mats, volume, nuclides)
|
||||
|
||||
# Create reaction rates array
|
||||
self.reaction_rates = ReactionRates(
|
||||
self.local_mats, self._burnable_nucs, self.chain.reactions)
|
||||
|
||||
def __call__(self, vec, power, print_out=True):
|
||||
"""Runs a simulation.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
vec : list of numpy.ndarray
|
||||
Total atoms to be used in function.
|
||||
power : float
|
||||
Power of the reactor in [W]
|
||||
print_out : bool, optional
|
||||
Whether or not to print out time.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.deplete.OperatorResult
|
||||
Eigenvalue and reaction rates resulting from transport operator
|
||||
|
||||
"""
|
||||
# Prevent OpenMC from complaining about re-creating tallies
|
||||
openmc.reset_auto_ids()
|
||||
|
||||
# Update status
|
||||
self.number.set_density(vec)
|
||||
|
||||
time_start = time.time()
|
||||
|
||||
# Update material compositions and tally nuclides
|
||||
self._update_materials()
|
||||
self._tally.nuclides = self._get_tally_nuclides()
|
||||
|
||||
# Run OpenMC
|
||||
openmc.capi.reset()
|
||||
openmc.capi.run()
|
||||
|
||||
time_openmc = time.time()
|
||||
|
||||
# Extract results
|
||||
op_result = self._unpack_tallies_and_normalize(power)
|
||||
|
||||
if comm.rank == 0:
|
||||
time_unpack = time.time()
|
||||
|
||||
if print_out:
|
||||
print("Time to openmc: ", time_openmc - time_start)
|
||||
print("Time to unpack: ", time_unpack - time_openmc)
|
||||
|
||||
return copy.deepcopy(op_result)
|
||||
|
||||
def _get_burnable_mats(self):
|
||||
"""Determine depletable materials, volumes, and nuclids
|
||||
|
||||
Returns
|
||||
-------
|
||||
burnable_mats : list of str
|
||||
List of burnable material IDs
|
||||
volume : OrderedDict of str to float
|
||||
Volume of each material in [cm^3]
|
||||
nuclides : list of str
|
||||
Nuclides in order of how they'll appear in the simulation.
|
||||
|
||||
"""
|
||||
burnable_mats = set()
|
||||
model_nuclides = set()
|
||||
volume = OrderedDict()
|
||||
|
||||
# Iterate once through the geometry to get dictionaries
|
||||
for mat in self.geometry.get_all_materials().values():
|
||||
for nuclide in mat.get_nuclides():
|
||||
model_nuclides.add(nuclide)
|
||||
if mat.depletable:
|
||||
burnable_mats.add(str(mat.id))
|
||||
if mat.volume is None:
|
||||
raise RuntimeError("Volume not specified for depletable "
|
||||
"material with ID={}.".format(mat.id))
|
||||
volume[str(mat.id)] = mat.volume
|
||||
|
||||
# Make sure there are burnable materials
|
||||
if not burnable_mats:
|
||||
raise RuntimeError(
|
||||
"No depletable materials were found in the model.")
|
||||
|
||||
# Sort the sets
|
||||
burnable_mats = sorted(burnable_mats, key=int)
|
||||
model_nuclides = sorted(model_nuclides)
|
||||
|
||||
# Construct a global nuclide dictionary, burned first
|
||||
nuclides = list(self.chain.nuclide_dict)
|
||||
for nuc in model_nuclides:
|
||||
if nuc not in nuclides:
|
||||
nuclides.append(nuc)
|
||||
|
||||
return burnable_mats, volume, nuclides
|
||||
|
||||
def _extract_number(self, local_mats, volume, nuclides):
|
||||
"""Construct AtomNumber using geometry
|
||||
|
||||
Parameters
|
||||
----------
|
||||
local_mats : list of str
|
||||
Material IDs to be managed by this process
|
||||
volume : OrderedDict of str to float
|
||||
Volumes for the above materials in [cm^3]
|
||||
nuclides : list of str
|
||||
Nuclides to be used in the simulation.
|
||||
|
||||
"""
|
||||
self.number = AtomNumber(local_mats, nuclides, volume, len(self.chain))
|
||||
|
||||
if self.dilute_initial != 0.0:
|
||||
for nuc in self._burnable_nucs:
|
||||
self.number.set_atom_density(np.s_[:], nuc, self.dilute_initial)
|
||||
|
||||
# Now extract the number densities and store
|
||||
for mat in self.geometry.get_all_materials().values():
|
||||
if str(mat.id) in local_mats:
|
||||
self._set_number_from_mat(mat)
|
||||
|
||||
def _set_number_from_mat(self, mat):
|
||||
"""Extracts material and number densities from openmc.Material
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mat : openmc.Material
|
||||
The material to read from
|
||||
|
||||
"""
|
||||
mat_id = str(mat.id)
|
||||
|
||||
for nuclide, density in mat.get_nuclide_atom_densities().values():
|
||||
number = density * 1.0e24
|
||||
self.number.set_atom_density(mat_id, nuclide, number)
|
||||
|
||||
def initial_condition(self):
|
||||
"""Performs final setup and returns initial condition.
|
||||
|
||||
Returns
|
||||
-------
|
||||
list of numpy.ndarray
|
||||
Total density for initial conditions.
|
||||
"""
|
||||
|
||||
# Create XML files
|
||||
if comm.rank == 0:
|
||||
self.geometry.export_to_xml()
|
||||
self.settings.export_to_xml()
|
||||
self._generate_materials_xml()
|
||||
|
||||
# Initialize OpenMC library
|
||||
comm.barrier()
|
||||
openmc.capi.init(comm)
|
||||
|
||||
# Generate tallies in memory
|
||||
self._generate_tallies()
|
||||
|
||||
# Return number density vector
|
||||
return list(self.number.get_mat_slice(np.s_[:]))
|
||||
|
||||
def finalize(self):
|
||||
"""Finalize a depletion simulation and release resources."""
|
||||
openmc.capi.finalize()
|
||||
|
||||
def _update_materials(self):
|
||||
"""Updates material compositions in OpenMC on all processes."""
|
||||
|
||||
for rank in range(comm.size):
|
||||
number_i = comm.bcast(self.number, root=rank)
|
||||
|
||||
for mat in number_i.materials:
|
||||
nuclides = []
|
||||
densities = []
|
||||
for nuc in number_i.nuclides:
|
||||
if nuc in self.nuclides_with_data:
|
||||
val = 1.0e-24 * number_i.get_atom_density(mat, nuc)
|
||||
|
||||
# If nuclide is zero, do not add to the problem.
|
||||
if val > 0.0:
|
||||
if self.round_number:
|
||||
val_magnitude = np.floor(np.log10(val))
|
||||
val_scaled = val / 10**val_magnitude
|
||||
val_round = round(val_scaled, 8)
|
||||
|
||||
val = val_round * 10**val_magnitude
|
||||
|
||||
nuclides.append(nuc)
|
||||
densities.append(val)
|
||||
else:
|
||||
# Only output warnings if values are significantly
|
||||
# negative. CRAM does not guarantee positive values.
|
||||
if val < -1.0e-21:
|
||||
print("WARNING: nuclide ", nuc, " in material ", mat,
|
||||
" is negative (density = ", val, " at/barn-cm)")
|
||||
number_i[mat, nuc] = 0.0
|
||||
|
||||
mat_internal = openmc.capi.materials[int(mat)]
|
||||
mat_internal.set_densities(nuclides, densities)
|
||||
|
||||
def _generate_materials_xml(self):
|
||||
"""Creates materials.xml from self.number.
|
||||
|
||||
Due to uncertainty with how MPI interacts with OpenMC API, this
|
||||
constructs the XML manually. The long term goal is to do this
|
||||
through direct memory writing.
|
||||
|
||||
"""
|
||||
materials = openmc.Materials(self.geometry.get_all_materials()
|
||||
.values())
|
||||
|
||||
# Sort nuclides according to order in AtomNumber object
|
||||
nuclides = list(self.number.nuclides)
|
||||
for mat in materials:
|
||||
mat._nuclides.sort(key=lambda x: nuclides.index(x[0]))
|
||||
|
||||
materials.export_to_xml()
|
||||
|
||||
def _get_tally_nuclides(self):
|
||||
"""Determine nuclides that should be tallied for reaction rates.
|
||||
|
||||
This method returns a list of all nuclides that have neutron data and
|
||||
are listed in the depletion chain. Technically, we should tally nuclides
|
||||
that may not appear in the depletion chain because we still need to get
|
||||
the fission reaction rate for these nuclides in order to normalize
|
||||
power, but that is left as a future exercise.
|
||||
|
||||
Returns
|
||||
-------
|
||||
list of str
|
||||
Tally nuclides
|
||||
|
||||
"""
|
||||
nuc_set = set()
|
||||
|
||||
# Create the set of all nuclides in the decay chain in materials marked
|
||||
# for burning in which the number density is greater than zero.
|
||||
for nuc in self.number.nuclides:
|
||||
if nuc in self.nuclides_with_data:
|
||||
if np.sum(self.number[:, nuc]) > 0.0:
|
||||
nuc_set.add(nuc)
|
||||
|
||||
# Communicate which nuclides have nonzeros to rank 0
|
||||
if comm.rank == 0:
|
||||
for i in range(1, comm.size):
|
||||
nuc_newset = comm.recv(source=i, tag=i)
|
||||
nuc_set |= nuc_newset
|
||||
else:
|
||||
comm.send(nuc_set, dest=0, tag=comm.rank)
|
||||
|
||||
if comm.rank == 0:
|
||||
# Sort nuclides in the same order as self.number
|
||||
nuc_list = [nuc for nuc in self.number.nuclides
|
||||
if nuc in nuc_set]
|
||||
else:
|
||||
nuc_list = None
|
||||
|
||||
# Store list of tally nuclides on each process
|
||||
nuc_list = comm.bcast(nuc_list)
|
||||
return [nuc for nuc in nuc_list if nuc in self.chain]
|
||||
|
||||
def _generate_tallies(self):
|
||||
"""Generates depletion tallies.
|
||||
|
||||
Using information from the depletion chain as well as the nuclides
|
||||
currently in the problem, this function automatically generates a
|
||||
tally.xml for the simulation.
|
||||
|
||||
"""
|
||||
# Create tallies for depleting regions
|
||||
materials = [openmc.capi.materials[int(i)]
|
||||
for i in self.burnable_mats]
|
||||
mat_filter = openmc.capi.MaterialFilter(materials)
|
||||
|
||||
# Set up a tally that has a material filter covering each depletable
|
||||
# material and scores corresponding to all reactions that cause
|
||||
# transmutation. The nuclides for the tally are set later when eval() is
|
||||
# called.
|
||||
self._tally = openmc.capi.Tally()
|
||||
self._tally.scores = self.chain.reactions
|
||||
self._tally.filters = [mat_filter]
|
||||
|
||||
def _unpack_tallies_and_normalize(self, power):
|
||||
"""Unpack tallies from OpenMC and return an operator result
|
||||
|
||||
This method uses OpenMC's C API bindings to determine the k-effective
|
||||
value and reaction rates from the simulation. The reaction rates are
|
||||
normalized by the user-specified power, summing the product of the
|
||||
fission reaction rate times the fission Q value for each material.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
power : float
|
||||
Power of the reactor in [W]
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.deplete.OperatorResult
|
||||
Eigenvalue and reaction rates resulting from transport operator
|
||||
|
||||
"""
|
||||
rates = self.reaction_rates
|
||||
rates[:, :, :] = 0.0
|
||||
|
||||
k_combined = openmc.capi.keff()[0]
|
||||
|
||||
# Extract tally bins
|
||||
materials = self.burnable_mats
|
||||
nuclides = self._tally.nuclides
|
||||
|
||||
# Form fast map
|
||||
nuc_ind = [rates.index_nuc[nuc] for nuc in nuclides]
|
||||
react_ind = [rates.index_rx[react] for react in self.chain.reactions]
|
||||
|
||||
# Compute fission power
|
||||
# TODO : improve this calculation
|
||||
|
||||
# Keep track of energy produced from all reactions in eV per source
|
||||
# particle
|
||||
energy = 0.0
|
||||
|
||||
# Create arrays to store fission Q values, reaction rates, and nuclide
|
||||
# numbers
|
||||
fission_Q = np.zeros(rates.n_nuc)
|
||||
rates_expanded = np.zeros((rates.n_nuc, rates.n_react))
|
||||
number = np.zeros(rates.n_nuc)
|
||||
|
||||
fission_ind = rates.index_rx["fission"]
|
||||
|
||||
for nuclide in self.chain.nuclides:
|
||||
if nuclide.name in rates.index_nuc:
|
||||
for rx in nuclide.reactions:
|
||||
if rx.type == 'fission':
|
||||
ind = rates.index_nuc[nuclide.name]
|
||||
fission_Q[ind] = rx.Q
|
||||
break
|
||||
|
||||
# Extract results
|
||||
for i, mat in enumerate(self.local_mats):
|
||||
# Get tally index
|
||||
slab = materials.index(mat)
|
||||
|
||||
# Get material results hyperslab
|
||||
results = self._tally.results[slab, :, 1]
|
||||
|
||||
# Zero out reaction rates and nuclide numbers
|
||||
rates_expanded[:] = 0.0
|
||||
number[:] = 0.0
|
||||
|
||||
# Expand into our memory layout
|
||||
j = 0
|
||||
for nuc, i_nuc_results in zip(nuclides, nuc_ind):
|
||||
number[i_nuc_results] = self.number[mat, nuc]
|
||||
for react in react_ind:
|
||||
rates_expanded[i_nuc_results, react] = results[j]
|
||||
j += 1
|
||||
|
||||
# Accumulate energy from fission
|
||||
energy += np.dot(rates_expanded[:, fission_ind], fission_Q)
|
||||
|
||||
# Divide by total number and store
|
||||
for i_nuc_results in nuc_ind:
|
||||
if number[i_nuc_results] != 0.0:
|
||||
for react in react_ind:
|
||||
rates_expanded[i_nuc_results, react] /= number[i_nuc_results]
|
||||
|
||||
rates[i, :, :] = rates_expanded
|
||||
|
||||
# Reduce energy produced from all processes
|
||||
energy = comm.allreduce(energy)
|
||||
|
||||
# Determine power in eV/s
|
||||
power /= JOULE_PER_EV
|
||||
|
||||
# Scale reaction rates to obtain units of reactions/sec
|
||||
rates *= power / energy
|
||||
|
||||
return OperatorResult(k_combined, rates)
|
||||
|
||||
def _get_nuclides_with_data(self):
|
||||
"""Loads a cross_sections.xml file to find participating nuclides.
|
||||
|
||||
This allows for nuclides that are important in the decay chain but not
|
||||
important neutronically, or have no cross section data.
|
||||
"""
|
||||
|
||||
# Reads cross_sections.xml to create a dictionary containing
|
||||
# participating (burning and not just decaying) nuclides.
|
||||
|
||||
try:
|
||||
filename = os.environ["OPENMC_CROSS_SECTIONS"]
|
||||
except KeyError:
|
||||
filename = None
|
||||
|
||||
nuclides = set()
|
||||
|
||||
try:
|
||||
tree = ET.parse(filename)
|
||||
except Exception:
|
||||
if filename is None:
|
||||
msg = "No cross_sections.xml specified in materials."
|
||||
else:
|
||||
msg = 'Cross section file "{}" is invalid.'.format(filename)
|
||||
raise IOError(msg)
|
||||
|
||||
root = tree.getroot()
|
||||
for nuclide_node in root.findall('library'):
|
||||
mats = nuclide_node.get('materials')
|
||||
if not mats:
|
||||
continue
|
||||
for name in mats.split():
|
||||
# Make a burn list of the union of nuclides in cross_sections.xml
|
||||
# and nuclides in depletion chain.
|
||||
if name not in nuclides:
|
||||
nuclides.add(name)
|
||||
|
||||
return nuclides
|
||||
|
||||
def get_results_info(self):
|
||||
"""Returns volume list, material lists, and nuc lists.
|
||||
|
||||
Returns
|
||||
-------
|
||||
volume : dict of str float
|
||||
Volumes corresponding to materials in full_burn_dict
|
||||
nuc_list : list of str
|
||||
A list of all nuclide names. Used for sorting the simulation.
|
||||
burn_list : list of int
|
||||
A list of all material IDs to be burned. Used for sorting the simulation.
|
||||
full_burn_list : list
|
||||
List of all burnable material IDs
|
||||
|
||||
"""
|
||||
nuc_list = self.number.burnable_nuclides
|
||||
burn_list = self.local_mats
|
||||
|
||||
volume = {}
|
||||
for i, mat in enumerate(burn_list):
|
||||
volume[mat] = self.number.volume[i]
|
||||
|
||||
# Combine volume dictionaries across processes
|
||||
volume_list = comm.allgather(volume)
|
||||
volume = {k: v for d in volume_list for k, v in d.items()}
|
||||
|
||||
return volume, nuc_list, burn_list, self.burnable_mats
|
||||
139
openmc/deplete/reaction_rates.py
Normal file
139
openmc/deplete/reaction_rates.py
Normal file
|
|
@ -0,0 +1,139 @@
|
|||
"""ReactionRates module.
|
||||
|
||||
An ndarray to store reaction rates with string, integer, or slice indexing.
|
||||
"""
|
||||
from collections import OrderedDict
|
||||
|
||||
import numpy as np
|
||||
|
||||
|
||||
class ReactionRates(np.ndarray):
|
||||
"""Reaction rates resulting from a transport operator call
|
||||
|
||||
This class is a subclass of :class:`numpy.ndarray` with a few custom
|
||||
attributes that make it easy to determine what index corresponds to a given
|
||||
material, nuclide, and reaction rate.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
local_mats : list of str
|
||||
Material IDs
|
||||
nuclides : list of str
|
||||
Depletable nuclides
|
||||
reactions : list of str
|
||||
Transmutation reactions being tracked
|
||||
|
||||
Attributes
|
||||
----------
|
||||
index_mat : OrderedDict of str to int
|
||||
A dictionary mapping material ID as string to index.
|
||||
index_nuc : OrderedDict of str to int
|
||||
A dictionary mapping nuclide name as string to index.
|
||||
index_rx : OrderedDict of str to int
|
||||
A dictionary mapping reaction name as string to index.
|
||||
n_mat : int
|
||||
Number of materials.
|
||||
n_nuc : int
|
||||
Number of nucs.
|
||||
n_react : int
|
||||
Number of reactions.
|
||||
|
||||
"""
|
||||
|
||||
# NumPy arrays can be created 1) explicitly 2) using view casting, and 3) by
|
||||
# slicing an existing array. Because of these possibilities, it's necessary
|
||||
# to put initialization logic in __new__ rather than __init__. Additionally,
|
||||
# subclasses need to handle the multiple ways of creating arrays by using
|
||||
# the __array_finalize__ method (discussed here:
|
||||
# https://docs.scipy.org/doc/numpy/user/basics.subclassing.html)
|
||||
|
||||
def __new__(cls, local_mats, nuclides, reactions):
|
||||
# Create appropriately-sized zeroed-out ndarray
|
||||
shape = (len(local_mats), len(nuclides), len(reactions))
|
||||
obj = super().__new__(cls, shape)
|
||||
obj[:] = 0.0
|
||||
|
||||
# Add mapping attributes
|
||||
obj.index_mat = {mat: i for i, mat in enumerate(local_mats)}
|
||||
obj.index_nuc = {nuc: i for i, nuc in enumerate(nuclides)}
|
||||
obj.index_rx = {rx: i for i, rx in enumerate(reactions)}
|
||||
|
||||
return obj
|
||||
|
||||
def __array_finalize__(self, obj):
|
||||
if obj is None:
|
||||
return
|
||||
self.index_mat = getattr(obj, 'index_mat', None)
|
||||
self.index_nuc = getattr(obj, 'index_nuc', None)
|
||||
self.index_rx = getattr(obj, 'index_rx', None)
|
||||
|
||||
# Reaction rates are distributed to other processes via multiprocessing,
|
||||
# which entails pickling the objects. In order to preserve the custom
|
||||
# attributes, we have to modify how the ndarray is pickled as described
|
||||
# here: https://stackoverflow.com/a/26599346/1572453
|
||||
|
||||
def __reduce__(self):
|
||||
state = super().__reduce__()
|
||||
new_state = state[2] + (self.index_mat, self.index_nuc, self.index_rx)
|
||||
return (state[0], state[1], new_state)
|
||||
|
||||
def __setstate__(self, state):
|
||||
self.index_mat = state[-3]
|
||||
self.index_nuc = state[-2]
|
||||
self.index_rx = state[-1]
|
||||
super().__setstate__(state[0:-3])
|
||||
|
||||
@property
|
||||
def n_mat(self):
|
||||
return len(self.index_mat)
|
||||
|
||||
@property
|
||||
def n_nuc(self):
|
||||
return len(self.index_nuc)
|
||||
|
||||
@property
|
||||
def n_react(self):
|
||||
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
|
||||
397
openmc/deplete/results.py
Normal file
397
openmc/deplete/results.py
Normal file
|
|
@ -0,0 +1,397 @@
|
|||
"""The results module.
|
||||
|
||||
Contains results generation and saving capabilities.
|
||||
"""
|
||||
|
||||
from collections import OrderedDict
|
||||
import copy
|
||||
|
||||
import numpy as np
|
||||
import h5py
|
||||
|
||||
from . import comm, have_mpi
|
||||
from .reaction_rates import ReactionRates
|
||||
|
||||
_VERSION_RESULTS = (1, 0)
|
||||
|
||||
|
||||
class Results(object):
|
||||
"""Output of a depletion run
|
||||
|
||||
Attributes
|
||||
----------
|
||||
k : list of float
|
||||
Eigenvalue for each substep.
|
||||
time : list of float
|
||||
Time at beginning, end of step, in seconds.
|
||||
n_mat : int
|
||||
Number of mats.
|
||||
n_nuc : int
|
||||
Number of nuclides.
|
||||
rates : list of ReactionRates
|
||||
The reaction rates for each substep.
|
||||
volume : OrderedDict of int to float
|
||||
Dictionary mapping mat id to volume.
|
||||
mat_to_ind : OrderedDict of str to int
|
||||
A dictionary mapping mat ID as string to index.
|
||||
nuc_to_ind : OrderedDict of str to int
|
||||
A dictionary mapping nuclide name as string to index.
|
||||
mat_to_hdf5_ind : OrderedDict of str to int
|
||||
A dictionary mapping mat ID as string to global index.
|
||||
n_hdf5_mats : int
|
||||
Number of materials in entire geometry.
|
||||
n_stages : int
|
||||
Number of stages in simulation.
|
||||
data : numpy.ndarray
|
||||
Atom quantity, stored by stage, mat, then by nuclide.
|
||||
|
||||
"""
|
||||
def __init__(self):
|
||||
self.k = None
|
||||
self.time = None
|
||||
self.rates = None
|
||||
self.volume = None
|
||||
|
||||
self.mat_to_ind = None
|
||||
self.nuc_to_ind = None
|
||||
self.mat_to_hdf5_ind = None
|
||||
|
||||
self.data = None
|
||||
|
||||
def __getitem__(self, pos):
|
||||
"""Retrieves an item from results.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
pos : tuple
|
||||
A three-length tuple containing a stage index, mat index and a nuc
|
||||
index. All can be integers or slices. The second two can be
|
||||
strings corresponding to their respective dictionary.
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
The atoms for stage, mat, nuc
|
||||
|
||||
"""
|
||||
stage, mat, nuc = pos
|
||||
if isinstance(mat, str):
|
||||
mat = self.mat_to_ind[mat]
|
||||
if isinstance(nuc, str):
|
||||
nuc = self.nuc_to_ind[nuc]
|
||||
|
||||
return self.data[stage, mat, nuc]
|
||||
|
||||
def __setitem__(self, pos, val):
|
||||
"""Sets an item from results.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
pos : tuple
|
||||
A three-length tuple containing a stage index, mat index and a nuc
|
||||
index. All can be integers or slices. The second two can be
|
||||
strings corresponding to their respective dictionary.
|
||||
|
||||
val : float
|
||||
The value to set data to.
|
||||
|
||||
"""
|
||||
stage, mat, nuc = pos
|
||||
if isinstance(mat, str):
|
||||
mat = self.mat_to_ind[mat]
|
||||
if isinstance(nuc, str):
|
||||
nuc = self.nuc_to_ind[nuc]
|
||||
|
||||
self.data[stage, mat, nuc] = val
|
||||
|
||||
@property
|
||||
def n_mat(self):
|
||||
return len(self.mat_to_ind)
|
||||
|
||||
@property
|
||||
def n_nuc(self):
|
||||
return len(self.nuc_to_ind)
|
||||
|
||||
@property
|
||||
def n_hdf5_mats(self):
|
||||
return len(self.mat_to_hdf5_ind)
|
||||
|
||||
@property
|
||||
def n_stages(self):
|
||||
return self.data.shape[0]
|
||||
|
||||
def allocate(self, volume, nuc_list, burn_list, full_burn_list, stages):
|
||||
"""Allocates memory of Results.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
volume : dict of str float
|
||||
Volumes corresponding to materials in full_burn_dict
|
||||
nuc_list : list of str
|
||||
A list of all nuclide names. Used for sorting the simulation.
|
||||
burn_list : list of int
|
||||
A list of all mat IDs to be burned. Used for sorting the simulation.
|
||||
full_burn_list : list of str
|
||||
List of all burnable material IDs
|
||||
stages : int
|
||||
Number of stages in simulation.
|
||||
|
||||
"""
|
||||
self.volume = copy.deepcopy(volume)
|
||||
self.nuc_to_ind = {nuc: i for i, nuc in enumerate(nuc_list)}
|
||||
self.mat_to_ind = {mat: i for i, mat in enumerate(burn_list)}
|
||||
self.mat_to_hdf5_ind = {mat: i for i, mat in enumerate(full_burn_list)}
|
||||
|
||||
# Create storage array
|
||||
self.data = np.zeros((stages, self.n_mat, self.n_nuc))
|
||||
|
||||
def export_to_hdf5(self, filename, step):
|
||||
"""Export results to an HDF5 file
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str
|
||||
The filename to write to
|
||||
step : int
|
||||
What step is this?
|
||||
|
||||
"""
|
||||
if have_mpi and h5py.get_config().mpi:
|
||||
kwargs = {'driver': 'mpio', 'comm': comm}
|
||||
else:
|
||||
kwargs = {}
|
||||
|
||||
kwargs['mode'] = "w" if step == 0 else "a"
|
||||
|
||||
with h5py.File(filename, **kwargs) as handle:
|
||||
self._to_hdf5(handle, step)
|
||||
|
||||
def _write_hdf5_metadata(self, handle):
|
||||
"""Writes result metadata in HDF5 file
|
||||
|
||||
Parameters
|
||||
----------
|
||||
handle : h5py.File or h5py.Group
|
||||
An hdf5 file or group type to store this in.
|
||||
|
||||
"""
|
||||
# Create and save the 5 dictionaries:
|
||||
# quantities
|
||||
# self.mat_to_ind -> self.volume (TODO: support for changing volumes)
|
||||
# self.nuc_to_ind
|
||||
# reactions
|
||||
# self.rates[0].nuc_to_ind (can be different from above, above is superset)
|
||||
# self.rates[0].react_to_ind
|
||||
# these are shared by every step of the simulation, and should be deduplicated.
|
||||
|
||||
# Store concentration mat and nuclide dictionaries (along with volumes)
|
||||
|
||||
handle.attrs['version'] = np.array(_VERSION_RESULTS)
|
||||
handle.attrs['filetype'] = np.string_('depletion results')
|
||||
|
||||
mat_list = sorted(self.mat_to_hdf5_ind, key=int)
|
||||
nuc_list = sorted(self.nuc_to_ind)
|
||||
rxn_list = sorted(self.rates[0].index_rx)
|
||||
|
||||
n_mats = self.n_hdf5_mats
|
||||
n_nuc_number = len(nuc_list)
|
||||
n_nuc_rxn = len(self.rates[0].index_nuc)
|
||||
n_rxn = len(rxn_list)
|
||||
n_stages = self.n_stages
|
||||
|
||||
mat_group = handle.create_group("materials")
|
||||
|
||||
for mat in mat_list:
|
||||
mat_single_group = mat_group.create_group(mat)
|
||||
mat_single_group.attrs["index"] = self.mat_to_hdf5_ind[mat]
|
||||
mat_single_group.attrs["volume"] = self.volume[mat]
|
||||
|
||||
nuc_group = handle.create_group("nuclides")
|
||||
|
||||
for nuc in nuc_list:
|
||||
nuc_single_group = nuc_group.create_group(nuc)
|
||||
nuc_single_group.attrs["atom number index"] = self.nuc_to_ind[nuc]
|
||||
if nuc in self.rates[0].index_nuc:
|
||||
nuc_single_group.attrs["reaction rate index"] = self.rates[0].index_nuc[nuc]
|
||||
|
||||
rxn_group = handle.create_group("reactions")
|
||||
|
||||
for rxn in rxn_list:
|
||||
rxn_single_group = rxn_group.create_group(rxn)
|
||||
rxn_single_group.attrs["index"] = self.rates[0].index_rx[rxn]
|
||||
|
||||
# Construct array storage
|
||||
|
||||
handle.create_dataset("number", (1, n_stages, n_mats, n_nuc_number),
|
||||
maxshape=(None, n_stages, n_mats, n_nuc_number),
|
||||
chunks=(1, 1, n_mats, n_nuc_number),
|
||||
dtype='float64')
|
||||
|
||||
handle.create_dataset("reaction rates", (1, n_stages, n_mats, n_nuc_rxn, n_rxn),
|
||||
maxshape=(None, n_stages, n_mats, n_nuc_rxn, n_rxn),
|
||||
chunks=(1, 1, n_mats, n_nuc_rxn, n_rxn),
|
||||
dtype='float64')
|
||||
|
||||
handle.create_dataset("eigenvalues", (1, n_stages),
|
||||
maxshape=(None, n_stages), dtype='float64')
|
||||
|
||||
handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
|
||||
|
||||
def _to_hdf5(self, handle, index):
|
||||
"""Converts results object into an hdf5 object.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
handle : h5py.File or h5py.Group
|
||||
An HDF5 file or group type to store this in.
|
||||
index : int
|
||||
What step is this?
|
||||
|
||||
"""
|
||||
if "/number" not in handle:
|
||||
comm.barrier()
|
||||
self._write_hdf5_metadata(handle)
|
||||
|
||||
comm.barrier()
|
||||
|
||||
# Grab handles
|
||||
number_dset = handle["/number"]
|
||||
rxn_dset = handle["/reaction rates"]
|
||||
eigenvalues_dset = handle["/eigenvalues"]
|
||||
time_dset = handle["/time"]
|
||||
|
||||
# Get number of results stored
|
||||
number_shape = list(number_dset.shape)
|
||||
number_results = number_shape[0]
|
||||
|
||||
new_shape = index + 1
|
||||
|
||||
if number_results < new_shape:
|
||||
# Extend first dimension by 1
|
||||
number_shape[0] = new_shape
|
||||
number_dset.resize(number_shape)
|
||||
|
||||
rxn_shape = list(rxn_dset.shape)
|
||||
rxn_shape[0] = new_shape
|
||||
rxn_dset.resize(rxn_shape)
|
||||
|
||||
eigenvalues_shape = list(eigenvalues_dset.shape)
|
||||
eigenvalues_shape[0] = new_shape
|
||||
eigenvalues_dset.resize(eigenvalues_shape)
|
||||
|
||||
time_shape = list(time_dset.shape)
|
||||
time_shape[0] = new_shape
|
||||
time_dset.resize(time_shape)
|
||||
|
||||
# If nothing to write, just return
|
||||
if len(self.mat_to_ind) == 0:
|
||||
return
|
||||
|
||||
# Add data
|
||||
# Note, for the last step, self.n_stages = 1, even if n_stages != 1.
|
||||
n_stages = self.n_stages
|
||||
inds = [self.mat_to_hdf5_ind[mat] for mat in self.mat_to_ind]
|
||||
low = min(inds)
|
||||
high = max(inds)
|
||||
for i in range(n_stages):
|
||||
number_dset[index, i, low:high+1, :] = self.data[i, :, :]
|
||||
rxn_dset[index, i, low:high+1, :, :] = self.rates[i][:, :, :]
|
||||
if comm.rank == 0:
|
||||
eigenvalues_dset[index, i] = self.k[i]
|
||||
if comm.rank == 0:
|
||||
time_dset[index, :] = self.time
|
||||
|
||||
@classmethod
|
||||
def from_hdf5(cls, handle, step):
|
||||
"""Loads results object from HDF5.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
handle : h5py.File or h5py.Group
|
||||
An HDF5 file or group type to load from.
|
||||
step : int
|
||||
What step is this?
|
||||
|
||||
"""
|
||||
results = cls()
|
||||
|
||||
# Grab handles
|
||||
number_dset = handle["/number"]
|
||||
eigenvalues_dset = handle["/eigenvalues"]
|
||||
time_dset = handle["/time"]
|
||||
|
||||
results.data = number_dset[step, :, :, :]
|
||||
results.k = eigenvalues_dset[step, :]
|
||||
results.time = time_dset[step, :]
|
||||
|
||||
# Reconstruct dictionaries
|
||||
results.volume = OrderedDict()
|
||||
results.mat_to_ind = OrderedDict()
|
||||
results.nuc_to_ind = OrderedDict()
|
||||
rxn_nuc_to_ind = OrderedDict()
|
||||
rxn_to_ind = OrderedDict()
|
||||
|
||||
for mat, mat_handle in handle["/materials"].items():
|
||||
vol = mat_handle.attrs["volume"]
|
||||
ind = mat_handle.attrs["index"]
|
||||
|
||||
results.volume[mat] = vol
|
||||
results.mat_to_ind[mat] = ind
|
||||
|
||||
for nuc, nuc_handle in handle["/nuclides"].items():
|
||||
ind_atom = nuc_handle.attrs["atom number index"]
|
||||
results.nuc_to_ind[nuc] = ind_atom
|
||||
|
||||
if "reaction rate index" in nuc_handle.attrs:
|
||||
rxn_nuc_to_ind[nuc] = nuc_handle.attrs["reaction rate index"]
|
||||
|
||||
for rxn, rxn_handle in handle["/reactions"].items():
|
||||
rxn_to_ind[rxn] = rxn_handle.attrs["index"]
|
||||
|
||||
results.rates = []
|
||||
# Reconstruct reactions
|
||||
for i in range(results.n_stages):
|
||||
rate = ReactionRates(results.mat_to_ind, rxn_nuc_to_ind, rxn_to_ind)
|
||||
|
||||
rate[:] = handle["/reaction rates"][step, i, :, :, :]
|
||||
results.rates.append(rate)
|
||||
|
||||
return results
|
||||
|
||||
@staticmethod
|
||||
def save(op, x, op_results, t, step_ind):
|
||||
"""Creates and writes depletion results to disk
|
||||
|
||||
Parameters
|
||||
----------
|
||||
op : openmc.deplete.TransportOperator
|
||||
The operator used to generate these results.
|
||||
x : list of list of numpy.array
|
||||
The prior x vectors. Indexed [i][cell] using the above equation.
|
||||
op_results : list of openmc.deplete.OperatorResult
|
||||
Results of applying transport operator
|
||||
t : list of float
|
||||
Time indices.
|
||||
step_ind : int
|
||||
Step index.
|
||||
|
||||
"""
|
||||
# Get indexing terms
|
||||
vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
|
||||
|
||||
# Create results
|
||||
stages = len(x)
|
||||
results = Results()
|
||||
results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages)
|
||||
|
||||
n_mat = len(burn_list)
|
||||
|
||||
for i in range(stages):
|
||||
for mat_i in range(n_mat):
|
||||
results[i, mat_i, :] = x[i][mat_i][:]
|
||||
|
||||
results.k = [r.k for r in op_results]
|
||||
results.rates = [r.rates for r in op_results]
|
||||
results.time = t
|
||||
|
||||
results.export_to_hdf5("depletion_results.h5", step_ind)
|
||||
105
openmc/deplete/results_list.py
Normal file
105
openmc/deplete/results_list.py
Normal file
|
|
@ -0,0 +1,105 @@
|
|||
import h5py
|
||||
import numpy as np
|
||||
|
||||
from .results import Results, _VERSION_RESULTS
|
||||
from openmc.checkvalue import check_filetype_version
|
||||
|
||||
|
||||
class ResultsList(list):
|
||||
"""A list of openmc.deplete.Results objects
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str
|
||||
The filename to read from.
|
||||
|
||||
"""
|
||||
def __init__(self, filename):
|
||||
super().__init__()
|
||||
with h5py.File(str(filename), "r") as fh:
|
||||
check_filetype_version(fh, 'depletion results', _VERSION_RESULTS[0])
|
||||
|
||||
# Get number of results stored
|
||||
n = fh["number"].value.shape[0]
|
||||
|
||||
for i in range(n):
|
||||
self.append(Results.from_hdf5(fh, i))
|
||||
|
||||
def get_atoms(self, mat, nuc):
|
||||
"""Get nuclide concentration over time from a single material
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mat : str
|
||||
Material name to evaluate
|
||||
nuc : str
|
||||
Nuclide name to evaluate
|
||||
|
||||
Returns
|
||||
-------
|
||||
time : numpy.ndarray
|
||||
Array of times in [s]
|
||||
concentration : numpy.ndarray
|
||||
Total number of atoms for specified nuclide
|
||||
|
||||
"""
|
||||
time = np.empty_like(self)
|
||||
concentration = np.empty_like(self)
|
||||
|
||||
# Evaluate value in each region
|
||||
for i, result in enumerate(self):
|
||||
time[i] = result.time[0]
|
||||
concentration[i] = result[0, mat, nuc]
|
||||
|
||||
return time, concentration
|
||||
|
||||
def get_reaction_rate(self, mat, nuc, rx):
|
||||
"""Get reaction rate in a single material/nuclide over time
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mat : str
|
||||
Material name to evaluate
|
||||
nuc : str
|
||||
Nuclide name to evaluate
|
||||
rx : str
|
||||
Reaction rate to evaluate
|
||||
|
||||
Returns
|
||||
-------
|
||||
time : numpy.ndarray
|
||||
Array of times in [s]
|
||||
rate : numpy.ndarray
|
||||
Array of reaction rates
|
||||
|
||||
"""
|
||||
time = np.empty_like(self)
|
||||
rate = np.empty_like(self)
|
||||
|
||||
# Evaluate value in each region
|
||||
for i, result in enumerate(self):
|
||||
time[i] = result.time[0]
|
||||
rate[i] = result.rates[0].get(mat, nuc, rx) * result[0, mat, nuc]
|
||||
|
||||
return time, rate
|
||||
|
||||
def get_eigenvalue(self):
|
||||
"""Evaluates the eigenvalue from a results list.
|
||||
|
||||
Returns
|
||||
-------
|
||||
time : numpy.ndarray
|
||||
Array of times in [s]
|
||||
eigenvalue : numpy.ndarray
|
||||
k-eigenvalue at each time
|
||||
|
||||
"""
|
||||
time = np.empty_like(self)
|
||||
eigenvalue = np.empty_like(self)
|
||||
|
||||
# Get time/eigenvalue at each point
|
||||
for i, result in enumerate(self):
|
||||
time[i] = result.time[0]
|
||||
eigenvalue[i] = result.k[0]
|
||||
|
||||
return time, eigenvalue
|
||||
|
|
@ -24,7 +24,7 @@ class Material(IDManagerMixin):
|
|||
To create a material, one should create an instance of this class, add
|
||||
nuclides or elements with :meth:`Material.add_nuclide` or
|
||||
`Material.add_element`, respectively, and set the total material density
|
||||
with `Material.export_to_xml()`. The material can then be assigned to a cell
|
||||
with `Material.set_density()`. The material can then be assigned to a cell
|
||||
using the :attr:`Cell.fill` attribute.
|
||||
|
||||
Parameters
|
||||
|
|
@ -52,8 +52,7 @@ class Material(IDManagerMixin):
|
|||
'atom/b-cm', 'atom/cm3', 'sum', or 'macro'. The 'macro' unit only
|
||||
applies in the case of a multi-group calculation.
|
||||
depletable : bool
|
||||
Indicate whether the material is depletable. This attribute can be used
|
||||
by downstream depletion applications.
|
||||
Indicate whether the material is depletable.
|
||||
nuclides : list of tuple
|
||||
List in which each item is a 3-tuple consisting of a nuclide string, the
|
||||
percent density, and the percent type ('ao' or 'wo').
|
||||
|
|
@ -74,6 +73,9 @@ class Material(IDManagerMixin):
|
|||
:meth:`Geometry.determine_paths` method.
|
||||
num_instances : int
|
||||
The number of instances of this material throughout the geometry.
|
||||
fissionable_mass : float
|
||||
Mass of fissionable nuclides in the material in [g]. Requires that the
|
||||
:attr:`volume` attribute is set.
|
||||
|
||||
"""
|
||||
|
||||
|
|
@ -245,6 +247,18 @@ class Material(IDManagerMixin):
|
|||
str)
|
||||
self._isotropic = list(isotropic)
|
||||
|
||||
@property
|
||||
def fissionable_mass(self):
|
||||
if self.volume is None:
|
||||
raise ValueError("Volume must be set in order to determine mass.")
|
||||
density = 0.0
|
||||
for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values():
|
||||
Z = openmc.data.zam(nuc)[0]
|
||||
if Z >= 90:
|
||||
density += 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
|
||||
/ openmc.data.AVOGADRO
|
||||
return density*self.volume
|
||||
|
||||
@classmethod
|
||||
def from_hdf5(cls, group):
|
||||
"""Create material from HDF5 group
|
||||
|
|
@ -687,6 +701,51 @@ class Material(IDManagerMixin):
|
|||
|
||||
return nuclides
|
||||
|
||||
def get_mass_density(self, nuclide=None):
|
||||
"""Return mass density of one or all nuclides
|
||||
|
||||
Parameters
|
||||
----------
|
||||
nuclides : str, optional
|
||||
Nuclide for which density is desired. If not specified, the density
|
||||
for the entire material is given.
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Density of the nuclide/material in [g/cm^3]
|
||||
|
||||
"""
|
||||
mass_density = 0.0
|
||||
for nuc, atoms_per_cc in self.get_nuclide_atom_densities().values():
|
||||
density_i = 1e24 * atoms_per_cc * openmc.data.atomic_mass(nuc) \
|
||||
/ openmc.data.AVOGADRO
|
||||
if nuclide is None or nuclide == nuc:
|
||||
mass_density += density_i
|
||||
return mass_density
|
||||
|
||||
def get_mass(self, nuclide=None):
|
||||
"""Return mass of one or all nuclides.
|
||||
|
||||
Note that this method requires that the :attr:`Material.volume` has
|
||||
already been set.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
nuclides : str, optional
|
||||
Nuclide for which mass is desired. If not specified, the density
|
||||
for the entire material is given.
|
||||
|
||||
Returns
|
||||
-------
|
||||
float
|
||||
Mass of the nuclide/material in [g]
|
||||
|
||||
"""
|
||||
if self.volume is None:
|
||||
raise ValueError("Volume must be set in order to determine mass.")
|
||||
return self.volume*self.get_mass_density(nuclide)
|
||||
|
||||
def clone(self, memo=None):
|
||||
"""Create a copy of this material with a new unique ID.
|
||||
|
||||
|
|
|
|||
|
|
@ -2,6 +2,10 @@ from collections.abc import Iterable
|
|||
|
||||
import openmc
|
||||
from openmc.checkvalue import check_type
|
||||
import openmc.deplete as dep
|
||||
|
||||
_DEPLETE_METHODS = {'predictor': dep.integrator.predictor,
|
||||
'cecm': dep.integrator.cecm}
|
||||
|
||||
|
||||
class Model(object):
|
||||
|
|
@ -136,9 +140,38 @@ class Model(object):
|
|||
for plot in plots:
|
||||
self._plots.append(plot)
|
||||
|
||||
def export_to_xml(self):
|
||||
"""Export model to XML files.
|
||||
def deplete(self, timesteps, power, chain_file=None, method='cecm', **kwargs):
|
||||
"""Deplete model using specified timesteps/power
|
||||
|
||||
Parameters
|
||||
----------
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
power : float or iterable of float
|
||||
Power of the reactor in [W]. A single value indicates that the power
|
||||
is constant over all timesteps. An iterable indicates potentially
|
||||
different power levels for each timestep. For a 2D problem, the
|
||||
power can be given in [W/cm] as long as the "volume" assigned to a
|
||||
depletion material is actually an area in [cm^2].
|
||||
chain_file : str, optional
|
||||
Path to the depletion chain XML file. Defaults to the
|
||||
:envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists.
|
||||
method : {'cecm', 'predictor'}
|
||||
Integration method used for depletion
|
||||
**kwargs
|
||||
Keyword arguments passed to integration function (e.g.,
|
||||
:func:`openmc.deplete.integrator.cecm`)
|
||||
|
||||
"""
|
||||
# Create OpenMC transport operator
|
||||
op = dep.Operator(self.geometry, self.settings, chain_file)
|
||||
|
||||
# Perform depletion
|
||||
_DEPLETE_METHODS[method](op, timesteps, power, **kwargs)
|
||||
|
||||
def export_to_xml(self):
|
||||
"""Export model to XML files."""
|
||||
|
||||
self.settings.export_to_xml()
|
||||
self.geometry.export_to_xml()
|
||||
|
|
@ -166,7 +199,7 @@ class Model(object):
|
|||
Parameters
|
||||
----------
|
||||
**kwargs
|
||||
All keyword arguments are passed to openmc.run
|
||||
All keyword arguments are passed to :func:`openmc.run`
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -176,9 +209,7 @@ class Model(object):
|
|||
"""
|
||||
self.export_to_xml()
|
||||
|
||||
return_code = openmc.run(**kwargs)
|
||||
|
||||
assert (return_code == 0), "OpenMC did not execute successfully"
|
||||
openmc.run(**kwargs)
|
||||
|
||||
n = self.settings.batches
|
||||
if self.settings.statepoint is not None:
|
||||
|
|
|
|||
33
scripts/openmc-make-depletion-chain
Executable file
33
scripts/openmc-make-depletion-chain
Executable file
|
|
@ -0,0 +1,33 @@
|
|||
#!/usr/bin/env python3
|
||||
|
||||
import glob
|
||||
import os
|
||||
from zipfile import ZipFile
|
||||
|
||||
from openmc._utils import download
|
||||
import openmc.deplete
|
||||
|
||||
|
||||
URLS = [
|
||||
'http://www.nndc.bnl.gov/endf/b7.1/zips/ENDF-B-VII.1-neutrons.zip',
|
||||
'http://www.nndc.bnl.gov/endf/b7.1/zips/ENDF-B-VII.1-decay.zip',
|
||||
'http://www.nndc.bnl.gov/endf/b7.1/zips/ENDF-B-VII.1-nfy.zip'
|
||||
]
|
||||
|
||||
def main():
|
||||
for url in URLS:
|
||||
basename = download(url)
|
||||
with ZipFile(basename, 'r') as zf:
|
||||
print('Extracting {}...'.format(basename))
|
||||
zf.extractall()
|
||||
|
||||
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'))
|
||||
|
||||
chain = openmc.deplete.Chain.from_endf(decay_files, nfy_files, neutron_files)
|
||||
chain.export_to_xml('chain_endfb71.xml')
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
main()
|
||||
|
|
@ -121,6 +121,8 @@ contains
|
|||
energy_min_neutron = ZERO
|
||||
entropy_on = .false.
|
||||
gen_per_batch = 1
|
||||
index_entropy_mesh = -1
|
||||
index_ufs_mesh = -1
|
||||
keff = ONE
|
||||
legendre_to_tabular = .true.
|
||||
legendre_to_tabular_points = 33
|
||||
|
|
|
|||
|
|
@ -732,8 +732,10 @@ contains
|
|||
integer :: MT
|
||||
character(C_CHAR), pointer :: string(:)
|
||||
character(len=:, kind=C_CHAR), allocatable :: score_
|
||||
logical :: depletion_rx
|
||||
|
||||
err = E_UNASSIGNED
|
||||
depletion_rx = .false.
|
||||
if (index >= 1 .and. index <= size(tallies)) then
|
||||
associate (t => tallies(index) % obj)
|
||||
if (allocated(t % score_bins)) deallocate(t % score_bins)
|
||||
|
|
@ -757,10 +759,13 @@ contains
|
|||
t % score_bins(i) = SCORE_NU_SCATTER
|
||||
case ('(n,2n)')
|
||||
t % score_bins(i) = N_2N
|
||||
depletion_rx = .true.
|
||||
case ('(n,3n)')
|
||||
t % score_bins(i) = N_3N
|
||||
depletion_rx = .true.
|
||||
case ('(n,4n)')
|
||||
t % score_bins(i) = N_4N
|
||||
depletion_rx = .true.
|
||||
case ('absorption')
|
||||
t % score_bins(i) = SCORE_ABSORPTION
|
||||
case ('fission', '18')
|
||||
|
|
@ -829,8 +834,10 @@ contains
|
|||
t % score_bins(i) = N_NC
|
||||
case ('(n,gamma)')
|
||||
t % score_bins(i) = N_GAMMA
|
||||
depletion_rx = .true.
|
||||
case ('(n,p)')
|
||||
t % score_bins(i) = N_P
|
||||
depletion_rx = .true.
|
||||
case ('(n,d)')
|
||||
t % score_bins(i) = N_D
|
||||
case ('(n,t)')
|
||||
|
|
@ -839,6 +846,7 @@ contains
|
|||
t % score_bins(i) = N_3HE
|
||||
case ('(n,a)')
|
||||
t % score_bins(i) = N_A
|
||||
depletion_rx = .true.
|
||||
case ('(n,2a)')
|
||||
t % score_bins(i) = N_2A
|
||||
case ('(n,3a)')
|
||||
|
|
@ -879,6 +887,7 @@ contains
|
|||
end do
|
||||
|
||||
err = 0
|
||||
t % depletion_rx = depletion_rx
|
||||
end associate
|
||||
else
|
||||
err = E_OUT_OF_BOUNDS
|
||||
|
|
|
|||
49
tests/chain_simple.xml
Normal file
49
tests/chain_simple.xml
Normal file
|
|
@ -0,0 +1,49 @@
|
|||
<?xml version="1.0"?>
|
||||
<depletion_chain>
|
||||
<nuclide name="I135" decay_modes="1" reactions="1" half_life="2.36520E+04">
|
||||
<decay type="beta" target="Xe135" branching_ratio="1.0" />
|
||||
<reaction type="(n,gamma)" Q="0.0" target="Xe136" /> <!-- Not precisely true, but whatever -->
|
||||
</nuclide>
|
||||
<nuclide name="Xe135" decay_modes="1" reactions="1" half_life="3.29040E+04">
|
||||
<decay type=" beta" target="Cs135" branching_ratio="1.0" />
|
||||
<reaction type="(n,gamma)" Q="0.0" target="Xe136" />
|
||||
</nuclide>
|
||||
<nuclide name="Xe136" decay_modes="0" reactions="0" />
|
||||
<nuclide name="Cs135" decay_modes="0" reactions="0" />
|
||||
<nuclide name="Gd157" decay_modes="0" reactions="1" >
|
||||
<reaction type="(n,gamma)" Q="0.0" target="Nothing" />
|
||||
</nuclide>
|
||||
<nuclide name="Gd156" decay_modes="0" reactions="1">
|
||||
<reaction type="(n,gamma)" Q="0.0" target="Gd157" />
|
||||
</nuclide>
|
||||
<nuclide name="U234" decay_modes="0" reactions="1">
|
||||
<reaction type="fission" Q="191840000."/>
|
||||
<neutron_fission_yields>
|
||||
<energies>2.53000e-02</energies>
|
||||
<fission_yields energy="2.53000e-02">
|
||||
<products>Gd157 Gd156 I135 Xe135 Xe136 Cs135</products>
|
||||
<data>1.093250e-04 2.087260e-04 2.780820e-02 6.759540e-03 2.392300e-02 4.356330e-05</data>
|
||||
</fission_yields>
|
||||
</neutron_fission_yields>
|
||||
</nuclide>
|
||||
<nuclide name="U235" decay_modes="0" reactions="1">
|
||||
<reaction type="fission" Q="193410000."/>
|
||||
<neutron_fission_yields>
|
||||
<energies>2.53000e-02</energies>
|
||||
<fission_yields energy="2.53000e-02">
|
||||
<products>Gd157 Gd156 I135 Xe135 Xe136 Cs135</products>
|
||||
<data>6.142710e-5 1.483250e-04 0.0292737 0.002566345 0.0219242 4.9097e-6</data>
|
||||
</fission_yields>
|
||||
</neutron_fission_yields>
|
||||
</nuclide>
|
||||
<nuclide name="U238" decay_modes="0" reactions="1">
|
||||
<reaction type="fission" Q="197790000."/>
|
||||
<neutron_fission_yields>
|
||||
<energies>2.53000e-02</energies>
|
||||
<fission_yields energy="2.53000e-02">
|
||||
<products>Gd157 Gd156 I135 Xe135 Xe136 Cs135</products>
|
||||
<data>4.141120e-04 7.605360e-04 0.0135457 0.00026864 0.0024432 3.7100E-07</data>
|
||||
</fission_yields>
|
||||
</neutron_fission_yields>
|
||||
</nuclide>
|
||||
</depletion_chain>
|
||||
|
|
@ -1,3 +1,5 @@
|
|||
import pytest
|
||||
|
||||
from tests.regression_tests import config as regression_config
|
||||
|
||||
|
||||
|
|
@ -15,3 +17,12 @@ def pytest_configure(config):
|
|||
for opt in opts:
|
||||
if config.getoption(opt) is not None:
|
||||
regression_config[opt] = config.getoption(opt)
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def run_in_tmpdir(tmpdir):
|
||||
orig = tmpdir.chdir()
|
||||
try:
|
||||
yield
|
||||
finally:
|
||||
orig.chdir()
|
||||
|
|
|
|||
156
tests/dummy_operator.py
Normal file
156
tests/dummy_operator.py
Normal file
|
|
@ -0,0 +1,156 @@
|
|||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
from openmc.deplete.reaction_rates import ReactionRates
|
||||
from openmc.deplete.abc import TransportOperator, OperatorResult
|
||||
|
||||
|
||||
class DummyOperator(TransportOperator):
|
||||
"""This is a dummy operator class with no statistical uncertainty.
|
||||
|
||||
y_1' = sin(y_2) y_1 + cos(y_1) y_2
|
||||
y_2' = -cos(y_2) y_1 + sin(y_1) y_2
|
||||
|
||||
y_1(0) = 1
|
||||
y_2(0) = 1
|
||||
|
||||
y_1(1.5) ~ 2.3197067076743316
|
||||
y_2(1.5) ~ 3.1726475740397628
|
||||
|
||||
"""
|
||||
def __init__(self):
|
||||
pass
|
||||
|
||||
def __call__(self, vec, power, print_out=False):
|
||||
"""Evaluates F(y)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
vec : list of numpy.array
|
||||
Total atoms to be used in function.
|
||||
power : float
|
||||
Power in [W]
|
||||
print_out : bool, optional, ignored
|
||||
Whether or not to print out time.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.deplete.OperatorResult
|
||||
Result of transport operator
|
||||
|
||||
"""
|
||||
mats = ["1"]
|
||||
nuclides = ["1", "2"]
|
||||
reactions = ["1"]
|
||||
|
||||
reaction_rates = ReactionRates(mats, nuclides, reactions)
|
||||
|
||||
reaction_rates[0, 0, 0] = vec[0][0]
|
||||
reaction_rates[0, 1, 0] = vec[0][1]
|
||||
|
||||
# Create a fake rates object
|
||||
return OperatorResult(0.0, reaction_rates)
|
||||
|
||||
@property
|
||||
def chain(self):
|
||||
return self
|
||||
|
||||
def form_matrix(self, rates):
|
||||
"""Forms the f(y) matrix in y' = f(y)y.
|
||||
|
||||
Nominally a depletion matrix, this is abstracted on the off chance
|
||||
that the function f has nothing to do with depletion at all.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
rates : numpy.ndarray
|
||||
Slice of reaction rates for a single material
|
||||
|
||||
Returns
|
||||
-------
|
||||
scipy.sparse.csr_matrix
|
||||
Sparse matrix representing f(y).
|
||||
"""
|
||||
|
||||
y_1 = rates[0, 0]
|
||||
y_2 = rates[1, 0]
|
||||
|
||||
mat = np.zeros((2, 2))
|
||||
a11 = np.sin(y_2)
|
||||
a12 = np.cos(y_1)
|
||||
a21 = -np.cos(y_2)
|
||||
a22 = np.sin(y_1)
|
||||
|
||||
return sp.csr_matrix(np.array([[a11, a12], [a21, a22]]))
|
||||
|
||||
@property
|
||||
def volume(self):
|
||||
"""
|
||||
volume : dict of str float
|
||||
Volumes of material
|
||||
"""
|
||||
|
||||
return {"1": 0.0}
|
||||
|
||||
@property
|
||||
def nuc_list(self):
|
||||
"""
|
||||
nuc_list : list of str
|
||||
A list of all nuclide names. Used for sorting the simulation.
|
||||
"""
|
||||
|
||||
return ["1", "2"]
|
||||
|
||||
@property
|
||||
def local_mats(self):
|
||||
"""
|
||||
local_mats : list of str
|
||||
A list of all material IDs to be burned. Used for sorting the simulation.
|
||||
"""
|
||||
|
||||
return ["1"]
|
||||
|
||||
@property
|
||||
def burnable_mats(self):
|
||||
"""Maps cell name to index in global geometry."""
|
||||
return self.local_mats
|
||||
|
||||
|
||||
@property
|
||||
def reaction_rates(self):
|
||||
"""
|
||||
reaction_rates : ReactionRates
|
||||
Reaction rates from the last operator step.
|
||||
"""
|
||||
mats = ["1"]
|
||||
nuclides = ["1", "2"]
|
||||
reactions = ["1"]
|
||||
|
||||
return ReactionRates(mats, nuclides, reactions)
|
||||
|
||||
def initial_condition(self):
|
||||
"""Returns initial vector.
|
||||
|
||||
Returns
|
||||
-------
|
||||
list of numpy.array
|
||||
Total density for initial conditions.
|
||||
"""
|
||||
|
||||
return [np.array((1.0, 1.0))]
|
||||
|
||||
def get_results_info(self):
|
||||
"""Returns volume list, cell lists, and nuc lists.
|
||||
|
||||
Returns
|
||||
-------
|
||||
volume : dict of str float
|
||||
Volumes corresponding to materials in full_burn_dict
|
||||
nuc_list : list of str
|
||||
A list of all nuclide names. Used for sorting the simulation.
|
||||
burn_list : list of int
|
||||
A list of all cell IDs to be burned. Used for sorting the simulation.
|
||||
full_burn_list : OrderedDict of str to int
|
||||
Maps cell name to index in global geometry.
|
||||
|
||||
"""
|
||||
return self.volume, self.nuc_list, self.local_mats, self.burnable_mats
|
||||
378
tests/regression_tests/example_geometry.py
Normal file
378
tests/regression_tests/example_geometry.py
Normal file
|
|
@ -0,0 +1,378 @@
|
|||
"""An example file showing how to make a geometry.
|
||||
|
||||
This particular example creates a 3x3 geometry, with 8 regular pins and one
|
||||
Gd-157 2 wt-percent enriched. All pins are segmented.
|
||||
"""
|
||||
|
||||
from collections import OrderedDict
|
||||
import math
|
||||
|
||||
import numpy as np
|
||||
import openmc
|
||||
|
||||
|
||||
def density_to_mat(dens_dict):
|
||||
"""Generates an OpenMC material from a cell ID and self.number_density.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
dens_dict : dict
|
||||
Dictionary mapping nuclide names to densities
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.Material
|
||||
The OpenMC material filled with nuclides.
|
||||
|
||||
"""
|
||||
mat = openmc.Material()
|
||||
for key in dens_dict:
|
||||
mat.add_nuclide(key, 1.0e-24*dens_dict[key])
|
||||
mat.set_density('sum')
|
||||
|
||||
return mat
|
||||
|
||||
|
||||
def generate_initial_number_density():
|
||||
""" Generates initial number density.
|
||||
|
||||
These results were from a CASMO5 run in which the gadolinium pin was
|
||||
loaded with 2 wt percent of Gd-157.
|
||||
"""
|
||||
|
||||
# Concentration to be used for all fuel pins
|
||||
fuel_dict = OrderedDict()
|
||||
fuel_dict['U235'] = 1.05692e21
|
||||
fuel_dict['U234'] = 1.00506e19
|
||||
fuel_dict['U238'] = 2.21371e22
|
||||
fuel_dict['O16'] = 4.62954e22
|
||||
fuel_dict['O17'] = 1.127684e20
|
||||
fuel_dict['I135'] = 1.0e10
|
||||
fuel_dict['Xe135'] = 1.0e10
|
||||
fuel_dict['Xe136'] = 1.0e10
|
||||
fuel_dict['Cs135'] = 1.0e10
|
||||
fuel_dict['Gd156'] = 1.0e10
|
||||
fuel_dict['Gd157'] = 1.0e10
|
||||
# fuel_dict['O18'] = 9.51352e19 # Does not exist in ENDF71, merged into 17
|
||||
|
||||
# Concentration to be used for the gadolinium fuel pin
|
||||
fuel_gd_dict = OrderedDict()
|
||||
fuel_gd_dict['U235'] = 1.03579e21
|
||||
fuel_gd_dict['U238'] = 2.16943e22
|
||||
fuel_gd_dict['Gd156'] = 3.95517E+10
|
||||
fuel_gd_dict['Gd157'] = 1.08156e20
|
||||
fuel_gd_dict['O16'] = 4.64035e22
|
||||
fuel_dict['I135'] = 1.0e10
|
||||
fuel_dict['Xe136'] = 1.0e10
|
||||
fuel_dict['Xe135'] = 1.0e10
|
||||
fuel_dict['Cs135'] = 1.0e10
|
||||
# There are a whole bunch of 1e-10 stuff here.
|
||||
|
||||
# Concentration to be used for cladding
|
||||
clad_dict = OrderedDict()
|
||||
clad_dict['O16'] = 3.07427e20
|
||||
clad_dict['O17'] = 7.48868e17
|
||||
clad_dict['Cr50'] = 3.29620e18
|
||||
clad_dict['Cr52'] = 6.35639e19
|
||||
clad_dict['Cr53'] = 7.20763e18
|
||||
clad_dict['Cr54'] = 1.79413e18
|
||||
clad_dict['Fe54'] = 5.57350e18
|
||||
clad_dict['Fe56'] = 8.74921e19
|
||||
clad_dict['Fe57'] = 2.02057e18
|
||||
clad_dict['Fe58'] = 2.68901e17
|
||||
clad_dict['Cr50'] = 3.29620e18
|
||||
clad_dict['Cr52'] = 6.35639e19
|
||||
clad_dict['Cr53'] = 7.20763e18
|
||||
clad_dict['Cr54'] = 1.79413e18
|
||||
clad_dict['Ni58'] = 2.51631e19
|
||||
clad_dict['Ni60'] = 9.69278e18
|
||||
clad_dict['Ni61'] = 4.21338e17
|
||||
clad_dict['Ni62'] = 1.34341e18
|
||||
clad_dict['Ni64'] = 3.43127e17
|
||||
clad_dict['Zr90'] = 2.18320e22
|
||||
clad_dict['Zr91'] = 4.76104e21
|
||||
clad_dict['Zr92'] = 7.27734e21
|
||||
clad_dict['Zr94'] = 7.37494e21
|
||||
clad_dict['Zr96'] = 1.18814e21
|
||||
clad_dict['Sn112'] = 4.67352e18
|
||||
clad_dict['Sn114'] = 3.17992e18
|
||||
clad_dict['Sn115'] = 1.63814e18
|
||||
clad_dict['Sn116'] = 7.00546e19
|
||||
clad_dict['Sn117'] = 3.70027e19
|
||||
clad_dict['Sn118'] = 1.16694e20
|
||||
clad_dict['Sn119'] = 4.13872e19
|
||||
clad_dict['Sn120'] = 1.56973e20
|
||||
clad_dict['Sn122'] = 2.23076e19
|
||||
clad_dict['Sn124'] = 2.78966e19
|
||||
|
||||
# Gap concentration
|
||||
# Funny enough, the example problem uses air.
|
||||
gap_dict = OrderedDict()
|
||||
gap_dict['O16'] = 7.86548e18
|
||||
gap_dict['O17'] = 2.99548e15
|
||||
gap_dict['N14'] = 3.38646e19
|
||||
gap_dict['N15'] = 1.23717e17
|
||||
|
||||
# Concentration to be used for coolant
|
||||
# No boron
|
||||
cool_dict = OrderedDict()
|
||||
cool_dict['H1'] = 4.68063e22
|
||||
cool_dict['O16'] = 2.33427e22
|
||||
cool_dict['O17'] = 8.89086e18
|
||||
|
||||
# Store these dictionaries in the initial conditions dictionary
|
||||
initial_density = OrderedDict()
|
||||
initial_density['fuel_gd'] = fuel_gd_dict
|
||||
initial_density['fuel'] = fuel_dict
|
||||
initial_density['gap'] = gap_dict
|
||||
initial_density['clad'] = clad_dict
|
||||
initial_density['cool'] = cool_dict
|
||||
|
||||
# Set up libraries to use
|
||||
temperature = OrderedDict()
|
||||
sab = OrderedDict()
|
||||
|
||||
# Toggle betweeen MCNP and NNDC data
|
||||
MCNP = False
|
||||
|
||||
if MCNP:
|
||||
temperature['fuel_gd'] = 900.0
|
||||
temperature['fuel'] = 900.0
|
||||
# We approximate temperature of everything as 600K, even though it was
|
||||
# actually 580K.
|
||||
temperature['gap'] = 600.0
|
||||
temperature['clad'] = 600.0
|
||||
temperature['cool'] = 600.0
|
||||
else:
|
||||
temperature['fuel_gd'] = 293.6
|
||||
temperature['fuel'] = 293.6
|
||||
temperature['gap'] = 293.6
|
||||
temperature['clad'] = 293.6
|
||||
temperature['cool'] = 293.6
|
||||
|
||||
sab['cool'] = 'c_H_in_H2O'
|
||||
|
||||
# Set up burnable materials
|
||||
burn = OrderedDict()
|
||||
burn['fuel_gd'] = True
|
||||
burn['fuel'] = True
|
||||
burn['gap'] = False
|
||||
burn['clad'] = False
|
||||
burn['cool'] = False
|
||||
|
||||
return temperature, sab, initial_density, burn
|
||||
|
||||
def segment_pin(n_rings, n_wedges, r_fuel, r_gap, r_clad):
|
||||
""" Calculates a segmented pin.
|
||||
|
||||
Separates a pin with n_rings and n_wedges. All cells have equal volume.
|
||||
Pin is centered at origin.
|
||||
"""
|
||||
|
||||
# Calculate all the volumes of interest
|
||||
v_fuel = math.pi * r_fuel**2
|
||||
v_gap = math.pi * r_gap**2 - v_fuel
|
||||
v_clad = math.pi * r_clad**2 - v_fuel - v_gap
|
||||
v_ring = v_fuel / n_rings
|
||||
v_segment = v_ring / n_wedges
|
||||
|
||||
# Compute ring radiuses
|
||||
r_rings = np.zeros(n_rings)
|
||||
|
||||
for i in range(n_rings):
|
||||
r_rings[i] = math.sqrt(1.0/(math.pi) * v_ring * (i+1))
|
||||
|
||||
# Compute thetas
|
||||
theta = np.linspace(0, 2*math.pi, n_wedges + 1)
|
||||
|
||||
# Compute surfaces
|
||||
fuel_rings = [openmc.ZCylinder(x0=0, y0=0, R=r_rings[i])
|
||||
for i in range(n_rings)]
|
||||
|
||||
fuel_wedges = [openmc.Plane(A=math.cos(theta[i]), B=math.sin(theta[i]))
|
||||
for i in range(n_wedges)]
|
||||
|
||||
gap_ring = openmc.ZCylinder(x0=0, y0=0, R=r_gap)
|
||||
clad_ring = openmc.ZCylinder(x0=0, y0=0, R=r_clad)
|
||||
|
||||
# Create cells
|
||||
fuel_cells = []
|
||||
if n_wedges == 1:
|
||||
for i in range(n_rings):
|
||||
cell = openmc.Cell(name='fuel')
|
||||
if i == 0:
|
||||
cell.region = -fuel_rings[0]
|
||||
else:
|
||||
cell.region = +fuel_rings[i-1] & -fuel_rings[i]
|
||||
fuel_cells.append(cell)
|
||||
else:
|
||||
for i in range(n_rings):
|
||||
for j in range(n_wedges):
|
||||
cell = openmc.Cell(name='fuel')
|
||||
if i == 0:
|
||||
if j != n_wedges-1:
|
||||
cell.region = (-fuel_rings[0]
|
||||
& +fuel_wedges[j]
|
||||
& -fuel_wedges[j+1])
|
||||
else:
|
||||
cell.region = (-fuel_rings[0]
|
||||
& +fuel_wedges[j]
|
||||
& -fuel_wedges[0])
|
||||
else:
|
||||
if j != n_wedges-1:
|
||||
cell.region = (+fuel_rings[i-1]
|
||||
& -fuel_rings[i]
|
||||
& +fuel_wedges[j]
|
||||
& -fuel_wedges[j+1])
|
||||
else:
|
||||
cell.region = (+fuel_rings[i-1]
|
||||
& -fuel_rings[i]
|
||||
& +fuel_wedges[j]
|
||||
& -fuel_wedges[0])
|
||||
fuel_cells.append(cell)
|
||||
|
||||
# Gap ring
|
||||
gap_cell = openmc.Cell(name='gap')
|
||||
gap_cell.region = +fuel_rings[-1] & -gap_ring
|
||||
fuel_cells.append(gap_cell)
|
||||
|
||||
# Clad ring
|
||||
clad_cell = openmc.Cell(name='clad')
|
||||
clad_cell.region = +gap_ring & -clad_ring
|
||||
fuel_cells.append(clad_cell)
|
||||
|
||||
# Moderator
|
||||
mod_cell = openmc.Cell(name='cool')
|
||||
mod_cell.region = +clad_ring
|
||||
fuel_cells.append(mod_cell)
|
||||
|
||||
# Form universe
|
||||
fuel_u = openmc.Universe()
|
||||
fuel_u.add_cells(fuel_cells)
|
||||
|
||||
return fuel_u, v_segment, v_gap, v_clad
|
||||
|
||||
def generate_geometry(n_rings, n_wedges):
|
||||
""" Generates example geometry.
|
||||
|
||||
This function creates the initial geometry, a 9 pin reflective problem.
|
||||
One pin, containing gadolinium, is discretized into sectors.
|
||||
|
||||
In addition to what one would do with the general OpenMC geometry code, it
|
||||
is necessary to create a dictionary, volume, that maps a cell ID to a
|
||||
volume. Further, by naming cells the same as the above materials, the code
|
||||
can automatically handle the mapping.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n_rings : int
|
||||
Number of rings to generate for the geometry
|
||||
n_wedges : int
|
||||
Number of wedges to generate for the geometry
|
||||
"""
|
||||
|
||||
pitch = 1.26197
|
||||
r_fuel = 0.412275
|
||||
r_gap = 0.418987
|
||||
r_clad = 0.476121
|
||||
|
||||
n_pin = 3
|
||||
|
||||
# This table describes the 'fuel' to actual type mapping
|
||||
# It's not necessary to do it this way. Just adjust the initial conditions
|
||||
# below.
|
||||
mapping = ['fuel', 'fuel', 'fuel',
|
||||
'fuel', 'fuel_gd', 'fuel',
|
||||
'fuel', 'fuel', 'fuel']
|
||||
|
||||
# Form pin cell
|
||||
fuel_u, v_segment, v_gap, v_clad = segment_pin(n_rings, n_wedges, r_fuel, r_gap, r_clad)
|
||||
|
||||
# Form lattice
|
||||
all_water_c = openmc.Cell(name='cool')
|
||||
all_water_u = openmc.Universe(cells=(all_water_c, ))
|
||||
|
||||
lattice = openmc.RectLattice()
|
||||
lattice.pitch = [pitch]*2
|
||||
lattice.lower_left = [-pitch*n_pin/2, -pitch*n_pin/2]
|
||||
lattice_array = [[fuel_u for i in range(n_pin)] for j in range(n_pin)]
|
||||
lattice.universes = lattice_array
|
||||
lattice.outer = all_water_u
|
||||
|
||||
# Bound universe
|
||||
x_low = openmc.XPlane(x0=-pitch*n_pin/2, boundary_type='reflective')
|
||||
x_high = openmc.XPlane(x0=pitch*n_pin/2, boundary_type='reflective')
|
||||
y_low = openmc.YPlane(y0=-pitch*n_pin/2, boundary_type='reflective')
|
||||
y_high = openmc.YPlane(y0=pitch*n_pin/2, boundary_type='reflective')
|
||||
z_low = openmc.ZPlane(z0=-10, boundary_type='reflective')
|
||||
z_high = openmc.ZPlane(z0=10, boundary_type='reflective')
|
||||
|
||||
# Compute bounding box
|
||||
lower_left = [-pitch*n_pin/2, -pitch*n_pin/2, -10]
|
||||
upper_right = [pitch*n_pin/2, pitch*n_pin/2, 10]
|
||||
|
||||
root_c = openmc.Cell(fill=lattice)
|
||||
root_c.region = (+x_low & -x_high
|
||||
& +y_low & -y_high
|
||||
& +z_low & -z_high)
|
||||
root_u = openmc.Universe(universe_id=0, cells=(root_c, ))
|
||||
geometry = openmc.Geometry(root_u)
|
||||
|
||||
v_cool = pitch**2 - (v_gap + v_clad + n_rings * n_wedges * v_segment)
|
||||
|
||||
# Store volumes for later usage
|
||||
volume = {'fuel': v_segment, 'gap':v_gap, 'clad':v_clad, 'cool':v_cool}
|
||||
|
||||
return geometry, volume, mapping, lower_left, upper_right
|
||||
|
||||
def generate_problem(n_rings=5, n_wedges=8):
|
||||
""" Merges geometry and materials.
|
||||
|
||||
This function initializes the materials for each cell using the dictionaries
|
||||
provided by generate_initial_number_density. It is assumed a cell named
|
||||
'fuel' will have further region differentiation (see mapping).
|
||||
|
||||
Parameters
|
||||
----------
|
||||
n_rings : int, optional
|
||||
Number of rings to generate for the geometry
|
||||
n_wedges : int, optional
|
||||
Number of wedges to generate for the geometry
|
||||
"""
|
||||
|
||||
# Get materials dictionary, geometry, and volumes
|
||||
temperature, sab, initial_density, burn = generate_initial_number_density()
|
||||
geometry, volume, mapping, lower_left, upper_right = generate_geometry(n_rings, n_wedges)
|
||||
|
||||
# Apply distribmats, fill geometry
|
||||
cells = geometry.root_universe.get_all_cells()
|
||||
for cell_id in cells:
|
||||
cell = cells[cell_id]
|
||||
if cell.name == 'fuel':
|
||||
|
||||
omc_mats = []
|
||||
|
||||
for cell_type in mapping:
|
||||
omc_mat = density_to_mat(initial_density[cell_type])
|
||||
|
||||
if cell_type in sab:
|
||||
omc_mat.add_s_alpha_beta(sab[cell_type])
|
||||
omc_mat.temperature = temperature[cell_type]
|
||||
omc_mat.depletable = burn[cell_type]
|
||||
omc_mat.volume = volume['fuel']
|
||||
|
||||
omc_mats.append(omc_mat)
|
||||
|
||||
cell.fill = omc_mats
|
||||
elif cell.name != '':
|
||||
omc_mat = density_to_mat(initial_density[cell.name])
|
||||
|
||||
if cell.name in sab:
|
||||
omc_mat.add_s_alpha_beta(sab[cell.name])
|
||||
omc_mat.temperature = temperature[cell.name]
|
||||
omc_mat.depletable = burn[cell.name]
|
||||
omc_mat.volume = volume[cell.name]
|
||||
|
||||
cell.fill = omc_mat
|
||||
|
||||
return geometry, lower_left, upper_right
|
||||
102
tests/regression_tests/test_deplete_full.py
Normal file
102
tests/regression_tests/test_deplete_full.py
Normal file
|
|
@ -0,0 +1,102 @@
|
|||
""" Full system test suite. """
|
||||
|
||||
from math import floor
|
||||
import shutil
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import openmc
|
||||
from openmc.data import JOULE_PER_EV
|
||||
import openmc.deplete
|
||||
|
||||
from tests.regression_tests import config
|
||||
from .example_geometry import generate_problem
|
||||
|
||||
|
||||
def test_full(run_in_tmpdir):
|
||||
"""Full system test suite.
|
||||
|
||||
Runs an entire OpenMC simulation with depletion coupling and verifies
|
||||
that the outputs match a reference file. Sensitive to changes in
|
||||
OpenMC.
|
||||
|
||||
This test runs a complete OpenMC simulation and tests the outputs.
|
||||
It will take a while.
|
||||
"""
|
||||
|
||||
n_rings = 2
|
||||
n_wedges = 4
|
||||
|
||||
# Load geometry from example
|
||||
geometry, lower_left, upper_right = generate_problem(n_rings, n_wedges)
|
||||
|
||||
# OpenMC-specific settings
|
||||
settings = openmc.Settings()
|
||||
settings.particles = 100
|
||||
settings.batches = 100
|
||||
settings.inactive = 40
|
||||
space = openmc.stats.Box(lower_left, upper_right)
|
||||
settings.source = openmc.Source(space=space)
|
||||
settings.seed = 1
|
||||
settings.verbosity = 3
|
||||
|
||||
# Create operator
|
||||
chain_file = Path(__file__).parents[1] / 'chain_simple.xml'
|
||||
op = openmc.deplete.Operator(geometry, settings, chain_file)
|
||||
op.round_number = True
|
||||
|
||||
# Power and timesteps
|
||||
dt1 = 15.*24*60*60 # 15 days
|
||||
dt2 = 1.5*30*24*60*60 # 1.5 months
|
||||
N = floor(dt2/dt1)
|
||||
dt = np.full(N, dt1)
|
||||
power = 2.337e15*4*JOULE_PER_EV*1e6 # MeV/second cm from CASMO
|
||||
|
||||
# Perform simulation using the predictor algorithm
|
||||
openmc.deplete.integrator.predictor(op, dt, power)
|
||||
|
||||
# Get path to test and reference results
|
||||
path_test = op.output_dir / 'depletion_results.h5'
|
||||
path_reference = Path(__file__).with_name('test_reference.h5')
|
||||
|
||||
# If updating results, do so and return
|
||||
if config['update']:
|
||||
shutil.copyfile(str(path_test), str(path_reference))
|
||||
return
|
||||
|
||||
# Load the reference/test results
|
||||
res_test = openmc.deplete.ResultsList(path_test)
|
||||
res_ref = openmc.deplete.ResultsList(path_reference)
|
||||
|
||||
# Assert same mats
|
||||
for mat in res_ref[0].mat_to_ind:
|
||||
assert mat in res_test[0].mat_to_ind, \
|
||||
"Material {} not in new results.".format(mat)
|
||||
for nuc in res_ref[0].nuc_to_ind:
|
||||
assert nuc in res_test[0].nuc_to_ind, \
|
||||
"Nuclide {} not in new results.".format(nuc)
|
||||
|
||||
for mat in res_test[0].mat_to_ind:
|
||||
assert mat in res_ref[0].mat_to_ind, \
|
||||
"Material {} not in old results.".format(mat)
|
||||
for nuc in res_test[0].nuc_to_ind:
|
||||
assert nuc in res_ref[0].nuc_to_ind, \
|
||||
"Nuclide {} not in old results.".format(nuc)
|
||||
|
||||
tol = 1.0e-6
|
||||
for mat in res_test[0].mat_to_ind:
|
||||
for nuc in res_test[0].nuc_to_ind:
|
||||
_, y_test = res_test.get_atoms(mat, nuc)
|
||||
_, y_old = res_ref.get_atoms(mat, nuc)
|
||||
|
||||
# Test each point
|
||||
correct = True
|
||||
for i, ref in enumerate(y_old):
|
||||
if ref != y_test[i]:
|
||||
if ref != 0.0:
|
||||
correct = np.abs(y_test[i] - ref) / ref <= tol
|
||||
else:
|
||||
correct = False
|
||||
|
||||
assert correct, "Discrepancy in mat {} and nuc {}\n{}\n{}".format(
|
||||
mat, nuc, y_old, y_test)
|
||||
BIN
tests/regression_tests/test_reference.h5
Normal file
BIN
tests/regression_tests/test_reference.h5
Normal file
Binary file not shown.
|
|
@ -2,15 +2,6 @@ import openmc
|
|||
import pytest
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def run_in_tmpdir(tmpdir):
|
||||
orig = tmpdir.chdir()
|
||||
try:
|
||||
yield
|
||||
finally:
|
||||
orig.chdir()
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def uo2():
|
||||
m = openmc.Material(material_id=100, name='UO2')
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
from collections import Mapping
|
||||
from collections.abc import Mapping
|
||||
import os
|
||||
|
||||
import numpy as np
|
||||
|
|
|
|||
|
|
@ -58,3 +58,21 @@ def test_water_density():
|
|||
assert dens(300.0, 3.0) == pytest.approx(1e-3/0.100215168e-2, 1e-6)
|
||||
assert dens(300.0, 80.0) == pytest.approx(1e-3/0.971180894e-3, 1e-6)
|
||||
assert dens(500.0, 3.0) == pytest.approx(1e-3/0.120241800e-2, 1e-6)
|
||||
|
||||
|
||||
def test_gnd_name():
|
||||
assert openmc.data.gnd_name(1, 1) == 'H1'
|
||||
assert openmc.data.gnd_name(40, 90) == ('Zr90')
|
||||
assert openmc.data.gnd_name(95, 242, 0) == ('Am242')
|
||||
assert openmc.data.gnd_name(95, 242, 1) == ('Am242_m1')
|
||||
assert openmc.data.gnd_name(95, 242, 10) == ('Am242_m10')
|
||||
|
||||
|
||||
def test_zam():
|
||||
assert openmc.data.zam('H1') == (1, 1, 0)
|
||||
assert openmc.data.zam('Zr90') == (40, 90, 0)
|
||||
assert openmc.data.zam('Am242') == (95, 242, 0)
|
||||
assert openmc.data.zam('Am242_m1') == (95, 242, 1)
|
||||
assert openmc.data.zam('Am242_m10') == (95, 242, 10)
|
||||
with pytest.raises(ValueError):
|
||||
openmc.data.zam('garbage')
|
||||
|
|
|
|||
147
tests/unit_tests/test_deplete_atom_number.py
Normal file
147
tests/unit_tests/test_deplete_atom_number.py
Normal file
|
|
@ -0,0 +1,147 @@
|
|||
""" Tests for the AtomNumber class """
|
||||
|
||||
import numpy as np
|
||||
from openmc.deplete import atom_number
|
||||
|
||||
|
||||
def test_indexing():
|
||||
"""Tests the __getitem__ and __setitem__ routines simultaneously."""
|
||||
|
||||
local_mats = ["10000", "10001"]
|
||||
nuclides = ["U238", "U235", "U234"]
|
||||
volume = {"10000" : 0.38, "10001" : 0.21}
|
||||
|
||||
number = atom_number.AtomNumber(local_mats, nuclides, volume, 2)
|
||||
|
||||
number["10000", "U238"] = 1.0
|
||||
number["10001", "U238"] = 2.0
|
||||
number["10000", "U235"] = 3.0
|
||||
number["10001", "U235"] = 4.0
|
||||
|
||||
# String indexing
|
||||
assert number["10000", "U238"] == 1.0
|
||||
assert number["10001", "U238"] == 2.0
|
||||
assert number["10000", "U235"] == 3.0
|
||||
assert number["10001", "U235"] == 4.0
|
||||
|
||||
# Int indexing
|
||||
assert number[0, 0] == 1.0
|
||||
assert number[1, 0] == 2.0
|
||||
assert number[0, 1] == 3.0
|
||||
assert number[1, 1] == 4.0
|
||||
|
||||
number[0, 0] = 5.0
|
||||
|
||||
assert number[0, 0] == 5.0
|
||||
assert number["10000", "U238"] == 5.0
|
||||
|
||||
|
||||
def test_properties():
|
||||
"""Test properties. """
|
||||
local_mats = ["10000", "10001"]
|
||||
nuclides = ["U238", "U235", "Gd157"]
|
||||
volume = {"10000" : 0.38, "10001" : 0.21}
|
||||
|
||||
number = atom_number.AtomNumber(local_mats, nuclides, volume, 2)
|
||||
|
||||
assert list(number.materials) == ["10000", "10001"]
|
||||
assert number.n_nuc == 3
|
||||
assert list(number.nuclides) == ["U238", "U235", "Gd157"]
|
||||
assert number.burnable_nuclides == ["U238", "U235"]
|
||||
|
||||
|
||||
def test_density_indexing():
|
||||
"""Tests the get and set_atom_density routines simultaneously."""
|
||||
|
||||
local_mats = ["10000", "10001", "10002"]
|
||||
nuclides = ["U238", "U235", "U234"]
|
||||
volume = {"10000" : 0.38, "10001" : 0.21}
|
||||
|
||||
number = atom_number.AtomNumber(local_mats, nuclides, volume, 2)
|
||||
|
||||
number.set_atom_density("10000", "U238", 1.0)
|
||||
number.set_atom_density("10001", "U238", 2.0)
|
||||
number.set_atom_density("10002", "U238", 3.0)
|
||||
number.set_atom_density("10000", "U235", 4.0)
|
||||
number.set_atom_density("10001", "U235", 5.0)
|
||||
number.set_atom_density("10002", "U235", 6.0)
|
||||
number.set_atom_density("10000", "U234", 7.0)
|
||||
number.set_atom_density("10001", "U234", 8.0)
|
||||
number.set_atom_density("10002", "U234", 9.0)
|
||||
|
||||
# String indexing
|
||||
assert number.get_atom_density("10000", "U238") == 1.0
|
||||
assert number.get_atom_density("10001", "U238") == 2.0
|
||||
assert number.get_atom_density("10002", "U238") == 3.0
|
||||
assert number.get_atom_density("10000", "U235") == 4.0
|
||||
assert number.get_atom_density("10001", "U235") == 5.0
|
||||
assert number.get_atom_density("10002", "U235") == 6.0
|
||||
assert number.get_atom_density("10000", "U234") == 7.0
|
||||
assert number.get_atom_density("10001", "U234") == 8.0
|
||||
assert number.get_atom_density("10002", "U234") == 9.0
|
||||
|
||||
# Int indexing
|
||||
assert number.get_atom_density(0, 0) == 1.0
|
||||
assert number.get_atom_density(1, 0) == 2.0
|
||||
assert number.get_atom_density(2, 0) == 3.0
|
||||
assert number.get_atom_density(0, 1) == 4.0
|
||||
assert number.get_atom_density(1, 1) == 5.0
|
||||
assert number.get_atom_density(2, 1) == 6.0
|
||||
assert number.get_atom_density(0, 2) == 7.0
|
||||
assert number.get_atom_density(1, 2) == 8.0
|
||||
assert number.get_atom_density(2, 2) == 9.0
|
||||
|
||||
|
||||
number.set_atom_density(0, 0, 5.0)
|
||||
assert number.get_atom_density(0, 0) == 5.0
|
||||
|
||||
# Verify volume is used correctly
|
||||
assert number[0, 0] == 5.0 * 0.38
|
||||
assert number[1, 0] == 2.0 * 0.21
|
||||
assert number[2, 0] == 3.0 * 1.0
|
||||
assert number[0, 1] == 4.0 * 0.38
|
||||
assert number[1, 1] == 5.0 * 0.21
|
||||
assert number[2, 1] == 6.0 * 1.0
|
||||
assert number[0, 2] == 7.0 * 0.38
|
||||
assert number[1, 2] == 8.0 * 0.21
|
||||
assert number[2, 2] == 9.0 * 1.0
|
||||
|
||||
|
||||
def test_get_mat_slice():
|
||||
"""Tests getting slices."""
|
||||
|
||||
local_mats = ["10000", "10001", "10002"]
|
||||
nuclides = ["U238", "U235", "U234"]
|
||||
volume = {"10000" : 0.38, "10001" : 0.21}
|
||||
|
||||
number = atom_number.AtomNumber(local_mats, nuclides, volume, 2)
|
||||
|
||||
number.number = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])
|
||||
|
||||
sl = number.get_mat_slice(0)
|
||||
|
||||
np.testing.assert_array_equal(sl, np.array([1.0, 2.0]))
|
||||
|
||||
sl = number.get_mat_slice("10000")
|
||||
|
||||
np.testing.assert_array_equal(sl, np.array([1.0, 2.0]))
|
||||
|
||||
|
||||
def test_set_mat_slice():
|
||||
"""Tests getting slices."""
|
||||
|
||||
local_mats = ["10000", "10001", "10002"]
|
||||
nuclides = ["U238", "U235", "U234"]
|
||||
volume = {"10000" : 0.38, "10001" : 0.21}
|
||||
|
||||
number = atom_number.AtomNumber(local_mats, nuclides, volume, 2)
|
||||
|
||||
number.set_mat_slice(0, [1.0, 2.0])
|
||||
|
||||
assert number[0, 0] == 1.0
|
||||
assert number[0, 1] == 2.0
|
||||
|
||||
number.set_mat_slice("10000", [3.0, 4.0])
|
||||
|
||||
assert number[0, 0] == 3.0
|
||||
assert number[0, 1] == 4.0
|
||||
37
tests/unit_tests/test_deplete_cecm.py
Normal file
37
tests/unit_tests/test_deplete_cecm.py
Normal file
|
|
@ -0,0 +1,37 @@
|
|||
"""Regression tests for openmc.deplete.integrator.cecm algorithm.
|
||||
|
||||
These tests integrate a simple test problem described in dummy_geometry.py.
|
||||
"""
|
||||
|
||||
from pytest import approx
|
||||
import openmc.deplete
|
||||
|
||||
from tests import dummy_operator
|
||||
|
||||
|
||||
def test_cecm(run_in_tmpdir):
|
||||
"""Integral regression test of integrator algorithm using CE/CM."""
|
||||
|
||||
op = dummy_operator.DummyOperator()
|
||||
op.output_dir = "test_integrator_regression"
|
||||
|
||||
# Perform simulation using the MCNPX/MCNP6 algorithm
|
||||
dt = [0.75, 0.75]
|
||||
power = 1.0
|
||||
openmc.deplete.cecm(op, dt, power, print_out=False)
|
||||
|
||||
# Load the files
|
||||
res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
|
||||
|
||||
_, y1 = res.get_atoms("1", "1")
|
||||
_, y2 = res.get_atoms("1", "2")
|
||||
|
||||
# Mathematica solution
|
||||
s1 = [1.86872629872102, 1.395525772416039]
|
||||
s2 = [2.18097439443550, 2.69429754646747]
|
||||
|
||||
assert y1[1] == approx(s1[0])
|
||||
assert y2[1] == approx(s1[1])
|
||||
|
||||
assert y1[2] == approx(s2[0])
|
||||
assert y2[2] == approx(s2[1])
|
||||
233
tests/unit_tests/test_deplete_chain.py
Normal file
233
tests/unit_tests/test_deplete_chain.py
Normal file
|
|
@ -0,0 +1,233 @@
|
|||
"""Tests for openmc.deplete.Chain class."""
|
||||
|
||||
from collections.abc import Mapping
|
||||
import os
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
from openmc.data import zam, ATOMIC_SYMBOL
|
||||
from openmc.deplete import comm, Chain, reaction_rates, nuclide
|
||||
import pytest
|
||||
|
||||
from tests import cdtemp
|
||||
|
||||
_ENDF_DATA = Path(os.environ['OPENMC_ENDF_DATA'])
|
||||
|
||||
_TEST_CHAIN = """\
|
||||
<depletion_chain>
|
||||
<nuclide name="A" half_life="23652.0" decay_modes="2" decay_energy="0.0" reactions="1">
|
||||
<decay type="beta1" target="B" branching_ratio="0.6"/>
|
||||
<decay type="beta2" target="C" branching_ratio="0.4"/>
|
||||
<reaction type="(n,gamma)" Q="0.0" target="C"/>
|
||||
</nuclide>
|
||||
<nuclide name="B" half_life="32904.0" decay_modes="1" decay_energy="0.0" reactions="1">
|
||||
<decay type="beta" target="A" branching_ratio="1.0"/>
|
||||
<reaction type="(n,gamma)" Q="0.0" target="C"/>
|
||||
</nuclide>
|
||||
<nuclide name="C" reactions="3">
|
||||
<reaction type="fission" Q="200000000.0"/>
|
||||
<reaction type="(n,gamma)" Q="0.0" target="A" branching_ratio="0.7"/>
|
||||
<reaction type="(n,gamma)" Q="0.0" target="B" branching_ratio="0.3"/>
|
||||
<neutron_fission_yields>
|
||||
<energies>0.0253</energies>
|
||||
<fission_yields energy="0.0253">
|
||||
<products>A B</products>
|
||||
<data>0.0292737 0.002566345</data>
|
||||
</fission_yields>
|
||||
</neutron_fission_yields>
|
||||
</nuclide>
|
||||
</depletion_chain>
|
||||
"""
|
||||
|
||||
|
||||
@pytest.fixture(scope='module')
|
||||
def simple_chain():
|
||||
with cdtemp():
|
||||
with open('chain_test.xml', 'w') as fh:
|
||||
fh.write(_TEST_CHAIN)
|
||||
yield Chain.from_xml('chain_test.xml')
|
||||
|
||||
|
||||
def test_init():
|
||||
"""Test depletion chain initialization."""
|
||||
chain = Chain()
|
||||
|
||||
assert isinstance(chain.nuclides, list)
|
||||
assert isinstance(chain.nuclide_dict, Mapping)
|
||||
|
||||
|
||||
def test_len():
|
||||
"""Test depletion chain length."""
|
||||
chain = Chain()
|
||||
chain.nuclides = ["NucA", "NucB", "NucC"]
|
||||
|
||||
assert len(chain) == 3
|
||||
|
||||
|
||||
def test_from_endf():
|
||||
"""Test depletion chain building from ENDF files"""
|
||||
decay_data = (_ENDF_DATA / 'decay').glob('*.endf')
|
||||
fpy_data = (_ENDF_DATA / 'nfy').glob('*.endf')
|
||||
neutron_data = (_ENDF_DATA / 'neutrons').glob('*.endf')
|
||||
chain = Chain.from_endf(decay_data, fpy_data, neutron_data)
|
||||
|
||||
assert len(chain) == len(chain.nuclides) == len(chain.nuclide_dict) == 3820
|
||||
for nuc in chain.nuclides:
|
||||
assert nuc == chain[nuc.name]
|
||||
|
||||
|
||||
def test_from_xml(simple_chain):
|
||||
"""Read chain_test.xml and ensure all values are correct."""
|
||||
# Unfortunately, this routine touches a lot of the code, but most of
|
||||
# the components external to depletion_chain.py are simple storage
|
||||
# types.
|
||||
|
||||
chain = simple_chain
|
||||
|
||||
# Basic checks
|
||||
assert len(chain) == 3
|
||||
|
||||
# A tests
|
||||
nuc = chain["A"]
|
||||
|
||||
assert nuc.name == "A"
|
||||
assert nuc.half_life == 2.36520E+04
|
||||
assert nuc.n_decay_modes == 2
|
||||
modes = nuc.decay_modes
|
||||
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]
|
||||
|
||||
# B tests
|
||||
nuc = chain["B"]
|
||||
|
||||
assert nuc.name == "B"
|
||||
assert nuc.half_life == 3.29040E+04
|
||||
assert nuc.n_decay_modes == 1
|
||||
modes = nuc.decay_modes
|
||||
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]
|
||||
|
||||
# C tests
|
||||
nuc = chain["C"]
|
||||
|
||||
assert nuc.name == "C"
|
||||
assert nuc.n_decay_modes == 0
|
||||
assert nuc.n_reaction_paths == 3
|
||||
assert [r.target for r in nuc.reactions] == [None, "A", "B"]
|
||||
assert [r.type for r in nuc.reactions] == ["fission", "(n,gamma)", "(n,gamma)"]
|
||||
assert [r.branching_ratio for r in nuc.reactions] == [1.0, 0.7, 0.3]
|
||||
|
||||
# Yield tests
|
||||
assert nuc.yield_energies == [0.0253]
|
||||
assert list(nuc.yield_data) == [0.0253]
|
||||
assert nuc.yield_data[0.0253] == [("A", 0.0292737), ("B", 0.002566345)]
|
||||
|
||||
|
||||
def test_export_to_xml(run_in_tmpdir):
|
||||
"""Test writing a depletion chain to XML."""
|
||||
|
||||
# Prevent different MPI ranks from conflicting
|
||||
filename = 'test{}.xml'.format(comm.rank)
|
||||
|
||||
A = nuclide.Nuclide()
|
||||
A.name = "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)]
|
||||
|
||||
B = nuclide.Nuclide()
|
||||
B.name = "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)]
|
||||
|
||||
C = nuclide.Nuclide()
|
||||
C.name = "C"
|
||||
C.reactions = [
|
||||
nuclide.ReactionTuple("fission", None, 2.0e8, 1.0),
|
||||
nuclide.ReactionTuple("(n,gamma)", "A", 0.0, 0.7),
|
||||
nuclide.ReactionTuple("(n,gamma)", "B", 0.0, 0.3)
|
||||
]
|
||||
C.yield_energies = [0.0253]
|
||||
C.yield_data = {0.0253: [("A", 0.0292737), ("B", 0.002566345)]}
|
||||
|
||||
chain = Chain()
|
||||
chain.nuclides = [A, B, C]
|
||||
chain.export_to_xml(filename)
|
||||
|
||||
chain_xml = open(filename, 'r').read()
|
||||
assert _TEST_CHAIN == chain_xml
|
||||
|
||||
|
||||
def test_form_matrix(simple_chain):
|
||||
""" Using chain_test, and a dummy reaction rate, compute the matrix. """
|
||||
# Relies on test_from_xml passing.
|
||||
|
||||
chain = simple_chain
|
||||
|
||||
mats = ["10000", "10001"]
|
||||
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", "B", "(n,gamma)", 3.0)
|
||||
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
|
||||
|
||||
# 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
|
||||
|
||||
# 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
|
||||
|
||||
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
|
||||
|
||||
|
||||
def test_getitem():
|
||||
""" Test nuc_by_ind converter function. """
|
||||
chain = Chain()
|
||||
chain.nuclides = ["NucA", "NucB", "NucC"]
|
||||
chain.nuclide_dict = {nuc: chain.nuclides.index(nuc)
|
||||
for nuc in chain.nuclides}
|
||||
|
||||
assert "NucA" == chain["NucA"]
|
||||
assert "NucB" == chain["NucB"]
|
||||
assert "NucC" == chain["NucC"]
|
||||
37
tests/unit_tests/test_deplete_cram.py
Normal file
37
tests/unit_tests/test_deplete_cram.py
Normal file
|
|
@ -0,0 +1,37 @@
|
|||
""" Tests for cram.py
|
||||
|
||||
Compares a few Mathematica matrix exponentials to CRAM16/CRAM48.
|
||||
"""
|
||||
|
||||
from pytest import approx
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
from openmc.deplete.integrator import CRAM16, CRAM48
|
||||
|
||||
|
||||
def test_CRAM16():
|
||||
"""Test 16-term CRAM."""
|
||||
x = np.array([1.0, 1.0])
|
||||
mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
|
||||
dt = 0.1
|
||||
|
||||
z = CRAM16(mat, x, dt)
|
||||
|
||||
# Solution from mathematica
|
||||
z0 = np.array((0.904837418035960, 0.576799023327476))
|
||||
|
||||
assert z == approx(z0)
|
||||
|
||||
|
||||
def test_CRAM48():
|
||||
"""Test 48-term CRAM."""
|
||||
x = np.array([1.0, 1.0])
|
||||
mat = sp.csr_matrix([[-1.0, 0.0], [-2.0, -3.0]])
|
||||
dt = 0.1
|
||||
|
||||
z = CRAM48(mat, x, dt)
|
||||
|
||||
# Solution from mathematica
|
||||
z0 = np.array((0.904837418035960, 0.576799023327476))
|
||||
|
||||
assert z == approx(z0)
|
||||
93
tests/unit_tests/test_deplete_integrator.py
Normal file
93
tests/unit_tests/test_deplete_integrator.py
Normal file
|
|
@ -0,0 +1,93 @@
|
|||
"""Tests for saving results
|
||||
|
||||
It is worth noting that openmc.deplete.integrate is extremely complex, to the
|
||||
point I am unsure if it can be reasonably unit-tested. For the time being, it
|
||||
will be left unimplemented and testing will be done via regression.
|
||||
|
||||
"""
|
||||
|
||||
import copy
|
||||
import os
|
||||
from unittest.mock import MagicMock
|
||||
|
||||
import numpy as np
|
||||
from openmc.deplete import (ReactionRates, Results, ResultsList, comm,
|
||||
OperatorResult)
|
||||
|
||||
|
||||
def test_results_save(run_in_tmpdir):
|
||||
"""Test data save module"""
|
||||
|
||||
stages = 3
|
||||
|
||||
np.random.seed(comm.rank)
|
||||
|
||||
# Mock geometry
|
||||
op = MagicMock()
|
||||
|
||||
vol_dict = {}
|
||||
full_burn_list = []
|
||||
|
||||
for i in range(comm.size):
|
||||
vol_dict[str(2*i)] = 1.2
|
||||
vol_dict[str(2*i + 1)] = 1.2
|
||||
full_burn_list.append(str(2*i))
|
||||
full_burn_list.append(str(2*i + 1))
|
||||
|
||||
burn_list = full_burn_list[2*comm.rank : 2*comm.rank + 2]
|
||||
nuc_list = ["na", "nb"]
|
||||
|
||||
op.get_results_info.return_value = vol_dict, nuc_list, burn_list, full_burn_list
|
||||
|
||||
# Construct x
|
||||
x1 = []
|
||||
x2 = []
|
||||
|
||||
for i in range(stages):
|
||||
x1.append([np.random.rand(2), np.random.rand(2)])
|
||||
x2.append([np.random.rand(2), np.random.rand(2)])
|
||||
|
||||
# Construct r
|
||||
r1 = ReactionRates(burn_list, ["na", "nb"], ["ra", "rb"])
|
||||
r1[:] = np.random.rand(2, 2, 2)
|
||||
|
||||
rate1 = []
|
||||
rate2 = []
|
||||
|
||||
for i in range(stages):
|
||||
rate1.append(copy.deepcopy(r1))
|
||||
r1[:] = np.random.rand(2, 2, 2)
|
||||
rate2.append(copy.deepcopy(r1))
|
||||
r1[:] = np.random.rand(2, 2, 2)
|
||||
|
||||
# Create global terms
|
||||
eigvl1 = np.random.rand(stages)
|
||||
eigvl2 = np.random.rand(stages)
|
||||
|
||||
eigvl1 = comm.bcast(eigvl1, root=0)
|
||||
eigvl2 = comm.bcast(eigvl2, root=0)
|
||||
|
||||
t1 = [0.0, 1.0]
|
||||
t2 = [1.0, 2.0]
|
||||
|
||||
op_result1 = [OperatorResult(k, rates) for k, rates in zip(eigvl1, rate1)]
|
||||
op_result2 = [OperatorResult(k, rates) for k, rates in zip(eigvl2, rate2)]
|
||||
Results.save(op, x1, op_result1, t1, 0)
|
||||
Results.save(op, x2, op_result2, t2, 1)
|
||||
|
||||
# Load the files
|
||||
res = ResultsList("depletion_results.h5")
|
||||
|
||||
for i in range(stages):
|
||||
for mat_i, mat in enumerate(burn_list):
|
||||
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], 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)
|
||||
|
||||
np.testing.assert_array_equal(res[1].k, eigvl2)
|
||||
np.testing.assert_array_equal(res[1].time, t2)
|
||||
116
tests/unit_tests/test_deplete_nuclide.py
Normal file
116
tests/unit_tests/test_deplete_nuclide.py
Normal file
|
|
@ -0,0 +1,116 @@
|
|||
"""Tests for the openmc.deplete.Nuclide class."""
|
||||
|
||||
import xml.etree.ElementTree as ET
|
||||
|
||||
from openmc.deplete import nuclide
|
||||
|
||||
|
||||
def test_n_decay_modes():
|
||||
""" Test the decay mode count parameter. """
|
||||
|
||||
nuc = nuclide.Nuclide()
|
||||
|
||||
nuc.decay_modes = [
|
||||
nuclide.DecayTuple("beta1", "a", 0.5),
|
||||
nuclide.DecayTuple("beta2", "b", 0.3),
|
||||
nuclide.DecayTuple("beta3", "c", 0.2)
|
||||
]
|
||||
|
||||
assert nuc.n_decay_modes == 3
|
||||
|
||||
|
||||
def test_n_reaction_paths():
|
||||
""" Test the reaction path count parameter. """
|
||||
|
||||
nuc = nuclide.Nuclide()
|
||||
|
||||
nuc.reactions = [
|
||||
nuclide.ReactionTuple("(n,2n)", "a", 0.0, 1.0),
|
||||
nuclide.ReactionTuple("(n,3n)", "b", 0.0, 1.0),
|
||||
nuclide.ReactionTuple("(n,4n)", "c", 0.0, 1.0)
|
||||
]
|
||||
|
||||
assert nuc.n_reaction_paths == 3
|
||||
|
||||
|
||||
def test_from_xml():
|
||||
"""Test reading nuclide data from an XML element."""
|
||||
|
||||
data = """
|
||||
<nuclide name="U235" reactions="2">
|
||||
<decay type="sf" target="U235" branching_ratio="7.2e-11"/>
|
||||
<decay type="alpha" target="Th231" branching_ratio="0.999999999928"/>
|
||||
<reaction type="(n,2n)" target="U234" Q="-5297781.0"/>
|
||||
<reaction type="(n,3n)" target="U233" Q="-12142300.0"/>
|
||||
<reaction type="(n,4n)" target="U232" Q="-17885600.0"/>
|
||||
<reaction type="(n,gamma)" target="U236" Q="6545200.0"/>
|
||||
<reaction type="fission" Q="193405400.0"/>
|
||||
<neutron_fission_yields>
|
||||
<energies>0.0253</energies>
|
||||
<fission_yields energy="0.0253">
|
||||
<products>Te134 Zr100 Xe138</products>
|
||||
<data>0.062155 0.0497641 0.0481413</data>
|
||||
</fission_yields>
|
||||
</neutron_fission_yields>
|
||||
</nuclide>
|
||||
"""
|
||||
|
||||
element = ET.fromstring(data)
|
||||
u235 = nuclide.Nuclide.from_xml(element)
|
||||
|
||||
assert u235.decay_modes == [
|
||||
nuclide.DecayTuple('sf', 'U235', 7.2e-11),
|
||||
nuclide.DecayTuple('alpha', 'Th231', 1 - 7.2e-11)
|
||||
]
|
||||
assert u235.reactions == [
|
||||
nuclide.ReactionTuple('(n,2n)', 'U234', -5297781.0, 1.0),
|
||||
nuclide.ReactionTuple('(n,3n)', 'U233', -12142300.0, 1.0),
|
||||
nuclide.ReactionTuple('(n,4n)', 'U232', -17885600.0, 1.0),
|
||||
nuclide.ReactionTuple('(n,gamma)', 'U236', 6545200.0, 1.0),
|
||||
nuclide.ReactionTuple('fission', None, 193405400.0, 1.0),
|
||||
]
|
||||
assert u235.yield_energies == [0.0253]
|
||||
assert u235.yield_data == {
|
||||
0.0253: [('Te134', 0.062155), ('Zr100', 0.0497641),
|
||||
('Xe138', 0.0481413)]
|
||||
}
|
||||
|
||||
|
||||
def test_to_xml_element():
|
||||
"""Test writing nuclide data to an XML element."""
|
||||
|
||||
C = nuclide.Nuclide()
|
||||
C.name = "C"
|
||||
C.half_life = 0.123
|
||||
C.decay_modes = [
|
||||
nuclide.DecayTuple('beta-', 'B', 0.99),
|
||||
nuclide.DecayTuple('alpha', 'D', 0.01)
|
||||
]
|
||||
C.reactions = [
|
||||
nuclide.ReactionTuple('fission', None, 2.0e8, 1.0),
|
||||
nuclide.ReactionTuple('(n,gamma)', 'A', 0.0, 1.0)
|
||||
]
|
||||
C.yield_energies = [0.0253]
|
||||
C.yield_data = {0.0253: [("A", 0.0292737), ("B", 0.002566345)]}
|
||||
element = C.to_xml_element()
|
||||
|
||||
assert element.get("half_life") == "0.123"
|
||||
|
||||
decay_elems = element.findall("decay")
|
||||
assert len(decay_elems) == 2
|
||||
assert decay_elems[0].get("type") == "beta-"
|
||||
assert decay_elems[0].get("target") == "B"
|
||||
assert decay_elems[0].get("branching_ratio") == "0.99"
|
||||
assert decay_elems[1].get("type") == "alpha"
|
||||
assert decay_elems[1].get("target") == "D"
|
||||
assert decay_elems[1].get("branching_ratio") == "0.01"
|
||||
|
||||
rx_elems = element.findall("reaction")
|
||||
assert len(rx_elems) == 2
|
||||
assert rx_elems[0].get("type") == "fission"
|
||||
assert float(rx_elems[0].get("Q")) == 2.0e8
|
||||
assert rx_elems[1].get("type") == "(n,gamma)"
|
||||
assert rx_elems[1].get("target") == "A"
|
||||
assert float(rx_elems[1].get("Q")) == 0.0
|
||||
|
||||
assert element.find('neutron_fission_yields') is not None
|
||||
37
tests/unit_tests/test_deplete_predictor.py
Normal file
37
tests/unit_tests/test_deplete_predictor.py
Normal file
|
|
@ -0,0 +1,37 @@
|
|||
"""Regression tests for openmc.deplete.integrator.predictor algorithm.
|
||||
|
||||
These tests integrate a simple test problem described in dummy_geometry.py.
|
||||
"""
|
||||
|
||||
from pytest import approx
|
||||
import openmc.deplete
|
||||
|
||||
from tests import dummy_operator
|
||||
|
||||
|
||||
def test_predictor(run_in_tmpdir):
|
||||
"""Integral regression test of integrator algorithm using predictor/corrector"""
|
||||
|
||||
op = dummy_operator.DummyOperator()
|
||||
op.output_dir = "test_integrator_regression"
|
||||
|
||||
# Perform simulation using the predictor algorithm
|
||||
dt = [0.75, 0.75]
|
||||
power = 1.0
|
||||
openmc.deplete.predictor(op, dt, power, print_out=False)
|
||||
|
||||
# Load the files
|
||||
res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
|
||||
|
||||
_, y1 = res.get_atoms("1", "1")
|
||||
_, y2 = res.get_atoms("1", "2")
|
||||
|
||||
# Mathematica solution
|
||||
s1 = [2.46847546272295, 0.986431226850467]
|
||||
s2 = [4.11525874568034, -0.0581692232513460]
|
||||
|
||||
assert y1[1] == approx(s1[0])
|
||||
assert y2[1] == approx(s1[1])
|
||||
|
||||
assert y1[2] == approx(s2[0])
|
||||
assert y2[2] == approx(s2[1])
|
||||
63
tests/unit_tests/test_deplete_reaction.py
Normal file
63
tests/unit_tests/test_deplete_reaction.py
Normal file
|
|
@ -0,0 +1,63 @@
|
|||
"""Tests for the openmc.deplete.ReactionRates class."""
|
||||
|
||||
import numpy as np
|
||||
from openmc.deplete import ReactionRates
|
||||
|
||||
|
||||
def test_get_set():
|
||||
"""Tests the get/set methods."""
|
||||
|
||||
local_mats = ["10000", "10001"]
|
||||
nuclides = ["U238", "U235"]
|
||||
reactions = ["fission", "(n,gamma)"]
|
||||
|
||||
rates = ReactionRates(local_mats, nuclides, reactions)
|
||||
assert rates.shape == (2, 2, 2)
|
||||
assert np.all(rates == 0.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.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
|
||||
assert rates[1, 0, 0] == 2.0
|
||||
assert rates[0, 1, 0] == 3.0
|
||||
assert rates[1, 1, 0] == 4.0
|
||||
assert rates[0, 0, 1] == 5.0
|
||||
assert rates[1, 0, 1] == 6.0
|
||||
assert rates[0, 1, 1] == 7.0
|
||||
assert rates[1, 1, 1] == 8.0
|
||||
|
||||
rates[0, 0, 0] = 5.0
|
||||
|
||||
assert rates[0, 0, 0] == 5.0
|
||||
assert rates.get("10000", "U238", "fission") == 5.0
|
||||
|
||||
|
||||
def test_properties():
|
||||
"""Test number of materials property."""
|
||||
local_mats = ["10000", "10001"]
|
||||
nuclides = ["U238", "U235", "Gd157"]
|
||||
reactions = ["fission", "(n,gamma)", "(n,2n)", "(n,3n)"]
|
||||
|
||||
rates = ReactionRates(local_mats, nuclides, reactions)
|
||||
|
||||
assert rates.n_mat == 2
|
||||
assert rates.n_nuc == 3
|
||||
assert rates.n_react == 4
|
||||
51
tests/unit_tests/test_deplete_resultslist.py
Normal file
51
tests/unit_tests/test_deplete_resultslist.py
Normal file
|
|
@ -0,0 +1,51 @@
|
|||
"""Tests the ResultsList class"""
|
||||
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import pytest
|
||||
import openmc.deplete
|
||||
|
||||
|
||||
@pytest.fixture
|
||||
def res():
|
||||
"""Load the reference results"""
|
||||
filename = Path(__file__).parents[1] / 'regression_tests' / 'test_reference.h5'
|
||||
return openmc.deplete.ResultsList(filename)
|
||||
|
||||
|
||||
def test_get_atoms(res):
|
||||
"""Tests evaluating single nuclide concentration."""
|
||||
t, n = res.get_atoms("1", "Xe135")
|
||||
|
||||
t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
|
||||
n_ref = [6.6747328233649218e+08, 3.5421791038348462e+14,
|
||||
3.6208592242443462e+14, 3.3799758969347038e+14]
|
||||
|
||||
np.testing.assert_array_equal(t, t_ref)
|
||||
np.testing.assert_array_equal(n, n_ref)
|
||||
|
||||
def test_get_reaction_rate(res):
|
||||
"""Tests evaluating reaction rate."""
|
||||
t, r = res.get_reaction_rate("1", "Xe135", "(n,gamma)")
|
||||
|
||||
t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
|
||||
n_ref = np.array([6.6747328233649218e+08, 3.5421791038348462e+14,
|
||||
3.6208592242443462e+14, 3.3799758969347038e+14])
|
||||
xs_ref = np.array([4.0594392323131994e-05, 3.9249546927524987e-05,
|
||||
3.8394587728581798e-05, 4.1521845978371697e-05])
|
||||
|
||||
np.testing.assert_array_equal(t, t_ref)
|
||||
np.testing.assert_array_equal(r, n_ref * xs_ref)
|
||||
|
||||
|
||||
def test_get_eigenvalue(res):
|
||||
"""Tests evaluating eigenvalue."""
|
||||
t, k = res.get_eigenvalue()
|
||||
|
||||
t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
|
||||
k_ref = [1.181281798790367, 1.1798750921988739, 1.1965943696058159,
|
||||
1.2207119847790813]
|
||||
|
||||
np.testing.assert_array_equal(t, t_ref)
|
||||
np.testing.assert_array_equal(k, k_ref)
|
||||
|
|
@ -121,6 +121,23 @@ def test_get_nuclide_atom_densities(uo2):
|
|||
assert density > 0
|
||||
|
||||
|
||||
def test_mass():
|
||||
m = openmc.Material()
|
||||
m.add_nuclide('Zr90', 1.0, 'wo')
|
||||
m.add_nuclide('U235', 1.0, 'wo')
|
||||
m.set_density('g/cm3', 2.0)
|
||||
m.volume = 10.0
|
||||
|
||||
assert m.get_mass_density('Zr90') == pytest.approx(1.0)
|
||||
assert m.get_mass_density('U235') == pytest.approx(1.0)
|
||||
assert m.get_mass_density() == pytest.approx(2.0)
|
||||
|
||||
assert m.get_mass('Zr90') == pytest.approx(10.0)
|
||||
assert m.get_mass('U235') == pytest.approx(10.0)
|
||||
assert m.get_mass() == pytest.approx(20.0)
|
||||
assert m.fissionable_mass == pytest.approx(10.0)
|
||||
|
||||
|
||||
def test_materials(run_in_tmpdir):
|
||||
m1 = openmc.Material()
|
||||
m1.add_nuclide('U235', 1.0, 'wo')
|
||||
|
|
|
|||
|
|
@ -14,10 +14,15 @@ pip install cython
|
|||
pip install --upgrade pytest
|
||||
|
||||
# Pandas stopped supporting Python 3.4 with version 0.21
|
||||
if [[ "$TRAVIS_PYTHON_VERSION" == "3.4" ]]; then
|
||||
if [[ $TRAVIS_PYTHON_VERSION == "3.4" ]]; then
|
||||
pip install pandas==0.20.3
|
||||
fi
|
||||
|
||||
# Install mpi4py for MPI configurations
|
||||
if [[ $MPI == 'y' ]]; then
|
||||
pip install --no-binary=mpi4py mpi4py
|
||||
fi
|
||||
|
||||
# Build and install OpenMC executable
|
||||
python tools/ci/travis-install.py
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue