mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-26 13:15:39 -04:00
Merge pull request #2077 from paulromano/depletion-renaming
Rename deplete.ResultsList → deplete.Results
This commit is contained in:
Andrew Johnson
GitHub
commit
64becd2c70
16 changed files with 936 additions and 910 deletions
|
|
@ -134,7 +134,7 @@ data, such as number densities and reaction rates for each material.
|
|||
OperatorResult
|
||||
ReactionRates
|
||||
Results
|
||||
ResultsList
|
||||
StepResult
|
||||
|
||||
The following class and functions are used to solve the depletion equations,
|
||||
with :func:`cram.CRAM48` being the default.
|
||||
|
|
|
|||
|
|
@ -49,11 +49,11 @@ one of these functions along with the timesteps and power level::
|
|||
|
||||
The coupled transport-depletion problem is executed, and once it is done a
|
||||
``depletion_results.h5`` file is written. The results can be analyzed using the
|
||||
:class:`openmc.deplete.ResultsList` class. This class has methods that allow for
|
||||
:class:`openmc.deplete.Results` class. This class has methods that allow for
|
||||
easy retrieval of k-effective, nuclide concentrations, and reaction rates over
|
||||
time::
|
||||
|
||||
results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
|
||||
results = openmc.deplete.Results("depletion_results.h5")
|
||||
time, keff = results.get_keff()
|
||||
|
||||
Note that the coupling between the transport solver and the transmutation solver
|
||||
|
|
|
|||
|
|
@ -12,7 +12,7 @@ with openmc.StatePoint(statepoint) as sp:
|
|||
geometry = sp.summary.geometry
|
||||
|
||||
# Load previous depletion results
|
||||
previous_results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
|
||||
previous_results = openmc.deplete.Results("depletion_results.h5")
|
||||
|
||||
###############################################################################
|
||||
# Transport calculation settings
|
||||
|
|
@ -56,7 +56,7 @@ integrator.integrate()
|
|||
###############################################################################
|
||||
|
||||
# Open results file
|
||||
results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
|
||||
results = openmc.deplete.Results("depletion_results.h5")
|
||||
|
||||
# Obtain K_eff as a function of time
|
||||
time, keff = results.get_keff(time_units='d')
|
||||
|
|
|
|||
|
|
@ -101,7 +101,7 @@ integrator.integrate()
|
|||
###############################################################################
|
||||
|
||||
# Open results file
|
||||
results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
|
||||
results = openmc.deplete.Results("depletion_results.h5")
|
||||
|
||||
# Obtain K_eff as a function of time
|
||||
time, keff = results.get_keff(time_units='d')
|
||||
|
|
|
|||
|
|
@ -10,8 +10,8 @@ from .chain import *
|
|||
from .operator import *
|
||||
from .reaction_rates import *
|
||||
from .atom_number import *
|
||||
from .stepresult import *
|
||||
from .results import *
|
||||
from .results_list import *
|
||||
from .integrators import *
|
||||
from . import abc
|
||||
from . import cram
|
||||
|
|
|
|||
|
|
@ -22,9 +22,9 @@ from uncertainties import ufloat
|
|||
from openmc.lib import MaterialFilter, Tally
|
||||
from openmc.checkvalue import check_type, check_greater_than
|
||||
from openmc.mpi import comm
|
||||
from .results import Results
|
||||
from .stepresult import StepResult
|
||||
from .chain import Chain
|
||||
from .results_list import ResultsList
|
||||
from .results import Results
|
||||
from .pool import deplete
|
||||
|
||||
|
||||
|
|
@ -79,7 +79,7 @@ class TransportOperator(ABC):
|
|||
in initial condition to ensure they exist in the decay chain.
|
||||
Only done for nuclides with reaction rates.
|
||||
Defaults to 1.0e3.
|
||||
prev_results : ResultsList, optional
|
||||
prev_results : Results, optional
|
||||
Results from a previous depletion calculation.
|
||||
|
||||
Attributes
|
||||
|
|
@ -88,7 +88,7 @@ class TransportOperator(ABC):
|
|||
Initial atom density [atoms/cm^3] 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.
|
||||
prev_res : ResultsList or None
|
||||
prev_res : Results or None
|
||||
Results from a previous depletion calculation. ``None`` if no
|
||||
results are to be used.
|
||||
"""
|
||||
|
|
@ -102,7 +102,7 @@ class TransportOperator(ABC):
|
|||
if prev_results is None:
|
||||
self.prev_res = None
|
||||
else:
|
||||
check_type("previous results", prev_results, ResultsList)
|
||||
check_type("previous results", prev_results, Results)
|
||||
self.prev_res = prev_results
|
||||
|
||||
@property
|
||||
|
|
@ -719,6 +719,9 @@ class Integrator(ABC):
|
|||
elif unit.lower() == 'mwd/kg':
|
||||
watt_days_per_kg = 1e6*timestep
|
||||
kilograms = 1e-3*operator.heavy_metal
|
||||
if rate == 0.0:
|
||||
raise ValueError("Cannot specify a timestep in [MWd/kg] when"
|
||||
" the power is zero.")
|
||||
days = watt_days_per_kg * kilograms / rate
|
||||
seconds.append(days*_SECONDS_PER_DAY)
|
||||
else:
|
||||
|
|
@ -842,7 +845,7 @@ class Integrator(ABC):
|
|||
k = ufloat(res.k[0, 0], res.k[0, 1])
|
||||
|
||||
# Scale reaction rates by ratio of source rates
|
||||
rates *= source_rate / res.source_rate[0]
|
||||
rates *= source_rate / res.source_rate
|
||||
return bos_conc, OperatorResult(k, rates)
|
||||
|
||||
def _get_start_data(self):
|
||||
|
|
@ -889,7 +892,7 @@ class Integrator(ABC):
|
|||
# Remove actual EOS concentration for next step
|
||||
conc = conc_list.pop()
|
||||
|
||||
Results.save(self.operator, conc_list, res_list, [t, t + dt],
|
||||
StepResult.save(self.operator, conc_list, res_list, [t, t + dt],
|
||||
source_rate, self._i_res + i, proc_time)
|
||||
|
||||
t += dt
|
||||
|
|
@ -901,7 +904,7 @@ class Integrator(ABC):
|
|||
if output and final_step:
|
||||
print(f"[openmc.deplete] t={t} (final operator evaluation)")
|
||||
res_list = [self.operator(conc, source_rate if final_step else 0.0)]
|
||||
Results.save(self.operator, [conc], res_list, [t, t],
|
||||
StepResult.save(self.operator, [conc], res_list, [t, t],
|
||||
source_rate, self._i_res + len(self), proc_time)
|
||||
self.operator.write_bos_data(len(self) + self._i_res)
|
||||
|
||||
|
|
@ -1048,13 +1051,13 @@ class SIIntegrator(Integrator):
|
|||
# Remove actual EOS concentration for next step
|
||||
conc = conc_list.pop()
|
||||
|
||||
Results.save(self.operator, conc_list, res_list, [t, t + dt],
|
||||
StepResult.save(self.operator, conc_list, res_list, [t, t + dt],
|
||||
p, self._i_res + i, proc_time)
|
||||
|
||||
t += dt
|
||||
|
||||
# No final simulation for SIE, use last iteration results
|
||||
Results.save(self.operator, [conc], [res_list[-1]], [t, t],
|
||||
StepResult.save(self.operator, [conc], [res_list[-1]], [t, t],
|
||||
p, self._i_res + len(self), proc_time)
|
||||
self.operator.write_bos_data(self._i_res + len(self))
|
||||
|
||||
|
|
|
|||
|
|
@ -25,7 +25,7 @@ from .abc import TransportOperator, OperatorResult
|
|||
from .atom_number import AtomNumber
|
||||
from .chain import _find_chain_file
|
||||
from .reaction_rates import ReactionRates
|
||||
from .results_list import ResultsList
|
||||
from .results import Results
|
||||
from .helpers import (
|
||||
DirectReactionRateHelper, ChainFissionHelper, ConstantFissionYieldHelper,
|
||||
FissionYieldCutoffHelper, AveragedFissionYieldHelper, EnergyScoreHelper,
|
||||
|
|
@ -94,7 +94,7 @@ class Operator(TransportOperator):
|
|||
Path to the depletion chain XML file. Defaults to the file
|
||||
listed under ``depletion_chain`` in
|
||||
:envvar:`OPENMC_CROSS_SECTIONS` environment variable.
|
||||
prev_results : ResultsList, optional
|
||||
prev_results : Results, optional
|
||||
Results from a previous depletion calculation. If this argument is
|
||||
specified, the depletion calculation will start from the latest state
|
||||
in the previous results.
|
||||
|
|
@ -192,7 +192,7 @@ class Operator(TransportOperator):
|
|||
Initial heavy metal inventory [g]
|
||||
local_mats : list of str
|
||||
All burnable material IDs being managed by a single process
|
||||
prev_res : ResultsList or None
|
||||
prev_res : Results or None
|
||||
Results from a previous depletion calculation. ``None`` if no
|
||||
results are to be used.
|
||||
cleanup_when_done : bool
|
||||
|
|
@ -284,7 +284,7 @@ class Operator(TransportOperator):
|
|||
if comm.size == 1:
|
||||
self.prev_res = prev_results
|
||||
else:
|
||||
self.prev_res = ResultsList()
|
||||
self.prev_res = Results()
|
||||
mat_indexes = _distribute(range(len(self.burnable_mats)))
|
||||
for res_obj in prev_results:
|
||||
new_res = res_obj.distribute(self.local_mats, mat_indexes)
|
||||
|
|
@ -494,7 +494,7 @@ class Operator(TransportOperator):
|
|||
Volumes for the above materials in [cm^3]
|
||||
nuclides : list of str
|
||||
Nuclides to be used in the simulation.
|
||||
prev_res : ResultsList, optional
|
||||
prev_res : Results, optional
|
||||
Results from a previous depletion calculation
|
||||
|
||||
"""
|
||||
|
|
@ -543,7 +543,7 @@ class Operator(TransportOperator):
|
|||
----------
|
||||
mat : openmc.Material
|
||||
The material to read from
|
||||
prev_res : ResultsList
|
||||
prev_res : Results
|
||||
Results from a previous depletion calculation
|
||||
|
||||
"""
|
||||
|
|
|
|||
|
|
@ -1,521 +1,416 @@
|
|||
"""The results module.
|
||||
|
||||
Contains results generation and saving capabilities.
|
||||
"""
|
||||
|
||||
from collections import OrderedDict
|
||||
import copy
|
||||
import numbers
|
||||
import bisect
|
||||
import math
|
||||
from warnings import warn
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
|
||||
import openmc
|
||||
from openmc.mpi import comm, MPI
|
||||
from .reaction_rates import ReactionRates
|
||||
from .stepresult import StepResult, VERSION_RESULTS
|
||||
import openmc.checkvalue as cv
|
||||
from openmc.data.library import DataLibrary
|
||||
from openmc.material import Material, Materials
|
||||
from openmc.exceptions import DataError
|
||||
|
||||
VERSION_RESULTS = (1, 1)
|
||||
__all__ = ["Results", "ResultsList"]
|
||||
|
||||
|
||||
__all__ = ["Results"]
|
||||
def _get_time_as(seconds, units):
|
||||
if units == "d":
|
||||
return seconds / (60 * 60 * 24)
|
||||
elif units == "h":
|
||||
return seconds / (60 * 60)
|
||||
elif units == "min":
|
||||
return seconds / 60
|
||||
else:
|
||||
return seconds
|
||||
|
||||
|
||||
class Results:
|
||||
"""Output of a depletion run
|
||||
class Results(list):
|
||||
"""Results from a depletion simulation
|
||||
|
||||
Attributes
|
||||
The :class:`Results` class acts as a list that stores the results from
|
||||
each depletion step and provides extra methods for interrogating these
|
||||
results.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
k : list of (float, float)
|
||||
Eigenvalue and uncertainty for each substep.
|
||||
time : list of float
|
||||
Time at beginning, end of step, in seconds.
|
||||
source_rate : float
|
||||
Source rate during timestep in [W] or [neutron/sec]
|
||||
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 str 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.
|
||||
proc_time : int
|
||||
Average time spent depleting a material across all
|
||||
materials and processes
|
||||
filename : str
|
||||
Path to depletion result file
|
||||
|
||||
"""
|
||||
def __init__(self):
|
||||
self.k = None
|
||||
self.time = None
|
||||
self.source_rate = None
|
||||
self.rates = None
|
||||
self.volume = None
|
||||
self.proc_time = None
|
||||
def __init__(self, filename=None):
|
||||
data = []
|
||||
if filename is not None:
|
||||
with h5py.File(str(filename), "r") as fh:
|
||||
cv.check_filetype_version(fh, 'depletion results', VERSION_RESULTS[0])
|
||||
|
||||
self.mat_to_ind = None
|
||||
self.nuc_to_ind = None
|
||||
self.mat_to_hdf5_ind = None
|
||||
# Get number of results stored
|
||||
n = fh["number"][...].shape[0]
|
||||
|
||||
self.data = None
|
||||
for i in range(n):
|
||||
data.append(StepResult.from_hdf5(fh, i))
|
||||
super().__init__(data)
|
||||
|
||||
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, openmc.Material):
|
||||
mat = str(mat.id)
|
||||
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 distribute(self, local_materials, ranges):
|
||||
"""Create a new object containing data for distributed materials
|
||||
|
||||
Parameters
|
||||
----------
|
||||
local_materials : iterable of str
|
||||
Materials for this process
|
||||
ranges : iterable of int
|
||||
Slice-like object indicating indicies of ``local_materials``
|
||||
in the material dimension of :attr:`data` and each element
|
||||
in :attr:`rates`
|
||||
|
||||
Returns
|
||||
-------
|
||||
Results
|
||||
New results object
|
||||
"""
|
||||
new = Results()
|
||||
new.volume = {lm: self.volume[lm] for lm in local_materials}
|
||||
new.mat_to_ind = {mat: idx for (idx, mat) in enumerate(local_materials)}
|
||||
|
||||
# Direct transfer
|
||||
direct_attrs = ("time", "k", "source_rate", "nuc_to_ind",
|
||||
"mat_to_hdf5_ind", "proc_time")
|
||||
for attr in direct_attrs:
|
||||
setattr(new, attr, getattr(self, attr))
|
||||
# Get applicable slice of data
|
||||
new.data = self.data[:, ranges]
|
||||
new.rates = [r[ranges] for r in self.rates]
|
||||
return new
|
||||
|
||||
def export_to_hdf5(self, filename, step):
|
||||
"""Export results to an HDF5 file
|
||||
@classmethod
|
||||
def from_hdf5(cls, filename):
|
||||
"""Load in depletion results from a previous file
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str
|
||||
The filename to write to
|
||||
step : int
|
||||
What step is this?
|
||||
Path to depletion result file
|
||||
|
||||
Returns
|
||||
-------
|
||||
Results
|
||||
New instance of depletion results
|
||||
|
||||
"""
|
||||
# Write new file if first time step, else add to existing file
|
||||
kwargs = {'mode': "w" if step == 0 else "a"}
|
||||
warn(
|
||||
"The ResultsList.from_hdf5(...) method is no longer necessary and will "
|
||||
"be removed in a future version of OpenMC. Use Results(...) instead.",
|
||||
FutureWarning
|
||||
)
|
||||
return cls(filename)
|
||||
|
||||
if h5py.get_config().mpi and comm.size > 1:
|
||||
# Write results in parallel
|
||||
kwargs['driver'] = 'mpio'
|
||||
kwargs['comm'] = comm
|
||||
with h5py.File(filename, **kwargs) as handle:
|
||||
self._to_hdf5(handle, step, parallel=True)
|
||||
def get_atoms(self, mat, nuc, nuc_units="atoms", time_units="s"):
|
||||
"""Get number of nuclides over time from a single material
|
||||
|
||||
.. note::
|
||||
Initial values for some isotopes that do not appear in
|
||||
initial concentrations may be non-zero, depending on the
|
||||
value of the :attr:`openmc.deplete.Operator.dilute_initial`
|
||||
attribute. The :class:`openmc.deplete.Operator` class adds isotopes
|
||||
according to this setting, which can be set to zero.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mat : openmc.Material, str
|
||||
Material object or material id to evaluate
|
||||
nuc : str
|
||||
Nuclide name to evaluate
|
||||
nuc_units : {"atoms", "atom/b-cm", "atom/cm3"}, optional
|
||||
Units for the returned concentration. Default is ``"atoms"``
|
||||
|
||||
.. versionadded:: 0.12
|
||||
time_units : {"s", "min", "h", "d"}, optional
|
||||
Units for the returned time array. Default is ``"s"`` to
|
||||
return the value in seconds.
|
||||
|
||||
.. versionadded:: 0.12
|
||||
|
||||
Returns
|
||||
-------
|
||||
times : numpy.ndarray
|
||||
Array of times in units of ``time_units``
|
||||
concentrations : numpy.ndarray
|
||||
Concentration of specified nuclide in units of ``nuc_units``
|
||||
|
||||
"""
|
||||
cv.check_value("time_units", time_units, {"s", "d", "min", "h"})
|
||||
cv.check_value("nuc_units", nuc_units,
|
||||
{"atoms", "atom/b-cm", "atom/cm3"})
|
||||
|
||||
if isinstance(mat, Material):
|
||||
mat_id = str(mat.id)
|
||||
elif isinstance(mat, str):
|
||||
mat_id = mat
|
||||
else:
|
||||
# Gather results at root process
|
||||
all_results = comm.gather(self)
|
||||
raise TypeError('mat should be of type openmc.Material or str')
|
||||
times = np.empty_like(self, dtype=float)
|
||||
concentrations = np.empty_like(self, dtype=float)
|
||||
|
||||
# Only root process writes results
|
||||
if comm.rank == 0:
|
||||
with h5py.File(filename, **kwargs) as handle:
|
||||
for res in all_results:
|
||||
res._to_hdf5(handle, step, parallel=False)
|
||||
# Evaluate value in each region
|
||||
for i, result in enumerate(self):
|
||||
times[i] = result.time[0]
|
||||
concentrations[i] = result[0, mat_id, nuc]
|
||||
|
||||
def _write_hdf5_metadata(self, handle):
|
||||
"""Writes result metadata in HDF5 file
|
||||
# Unit conversions
|
||||
times = _get_time_as(times, time_units)
|
||||
if nuc_units != "atoms":
|
||||
# Divide by volume to get density
|
||||
concentrations /= self[0].volume[mat_id]
|
||||
if nuc_units == "atom/b-cm":
|
||||
# 1 barn = 1e-24 cm^2
|
||||
concentrations *= 1e-24
|
||||
|
||||
return times, concentrations
|
||||
|
||||
def get_reaction_rate(self, mat, nuc, rx):
|
||||
"""Get reaction rate in a single material/nuclide over time
|
||||
|
||||
.. note::
|
||||
|
||||
Initial values for some isotopes that do not appear in
|
||||
initial concentrations may be non-zero, depending on the
|
||||
value of :class:`openmc.deplete.Operator` ``dilute_initial``
|
||||
The :class:`openmc.deplete.Operator` adds isotopes according
|
||||
to this setting, which can be set to zero.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
handle : h5py.File or h5py.Group
|
||||
An hdf5 file or group type to store this in.
|
||||
mat : openmc.Material, str
|
||||
Material object or material id to evaluate
|
||||
nuc : str
|
||||
Nuclide name to evaluate
|
||||
rx : str
|
||||
Reaction rate to evaluate
|
||||
|
||||
Returns
|
||||
-------
|
||||
times : numpy.ndarray
|
||||
Array of times in [s]
|
||||
rates : numpy.ndarray
|
||||
Array of reaction rates
|
||||
|
||||
"""
|
||||
# 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.
|
||||
times = np.empty_like(self, dtype=float)
|
||||
rates = np.empty_like(self, dtype=float)
|
||||
|
||||
# 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, 2),
|
||||
maxshape=(None, n_stages, 2), dtype='float64')
|
||||
|
||||
handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
|
||||
|
||||
handle.create_dataset("source_rate", (1, n_stages), maxshape=(None, n_stages),
|
||||
dtype='float64')
|
||||
|
||||
handle.create_dataset(
|
||||
"depletion time", (1,), maxshape=(None,),
|
||||
dtype="float64")
|
||||
|
||||
def _to_hdf5(self, handle, index, parallel=False):
|
||||
"""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?
|
||||
parallel : bool
|
||||
Being called with parallel HDF5?
|
||||
|
||||
"""
|
||||
if "/number" not in handle:
|
||||
if parallel:
|
||||
comm.barrier()
|
||||
self._write_hdf5_metadata(handle)
|
||||
|
||||
if parallel:
|
||||
comm.barrier()
|
||||
|
||||
# Grab handles
|
||||
number_dset = handle["/number"]
|
||||
rxn_dset = handle["/reaction rates"]
|
||||
eigenvalues_dset = handle["/eigenvalues"]
|
||||
time_dset = handle["/time"]
|
||||
source_rate_dset = handle["/source_rate"]
|
||||
proc_time_dset = handle["/depletion 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)
|
||||
|
||||
source_rate_shape = list(source_rate_dset.shape)
|
||||
source_rate_shape[0] = new_shape
|
||||
source_rate_dset.resize(source_rate_shape)
|
||||
|
||||
proc_shape = list(proc_time_dset.shape)
|
||||
proc_shape[0] = new_shape
|
||||
proc_time_dset.resize(proc_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
|
||||
source_rate_dset[index] = self.source_rate
|
||||
if self.proc_time is not None:
|
||||
proc_time_dset[index] = (
|
||||
self.proc_time / (comm.size * self.n_hdf5_mats)
|
||||
)
|
||||
|
||||
@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
|
||||
Index for depletion step
|
||||
"""
|
||||
results = cls()
|
||||
|
||||
# Grab handles
|
||||
number_dset = handle["/number"]
|
||||
eigenvalues_dset = handle["/eigenvalues"]
|
||||
time_dset = handle["/time"]
|
||||
if "source_rate" in handle:
|
||||
source_rate_dset = handle["/source_rate"]
|
||||
if isinstance(mat, Material):
|
||||
mat_id = str(mat.id)
|
||||
elif isinstance(mat, str):
|
||||
mat_id = mat
|
||||
else:
|
||||
# Older versions used "power" instead of "source_rate"
|
||||
source_rate_dset = handle["/power"]
|
||||
raise TypeError('mat should be of type openmc.Material or str')
|
||||
|
||||
results.data = number_dset[step, :, :, :]
|
||||
results.k = eigenvalues_dset[step, :]
|
||||
results.time = time_dset[step, :]
|
||||
results.source_rate = source_rate_dset[step, :]
|
||||
# Evaluate value in each region
|
||||
for i, result in enumerate(self):
|
||||
times[i] = result.time[0]
|
||||
rates[i] = result.rates[0].get(mat_id, nuc, rx) * result[0, mat, nuc]
|
||||
|
||||
if "depletion time" in handle:
|
||||
proc_time_dset = handle["/depletion time"]
|
||||
if step < proc_time_dset.shape[0]:
|
||||
results.proc_time = proc_time_dset[step]
|
||||
return times, rates
|
||||
|
||||
if results.proc_time is None:
|
||||
results.proc_time = np.array([np.nan])
|
||||
def get_keff(self, time_units='s'):
|
||||
"""Evaluates the eigenvalue from a results list.
|
||||
|
||||
# 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, True)
|
||||
|
||||
rate[:] = handle["/reaction rates"][step, i, :, :, :]
|
||||
results.rates.append(rate)
|
||||
|
||||
return results
|
||||
|
||||
@staticmethod
|
||||
def save(op, x, op_results, t, source_rate, step_ind, proc_time=None):
|
||||
"""Creates and writes depletion results to disk
|
||||
.. versionadded:: 0.13.1
|
||||
|
||||
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.
|
||||
source_rate : float
|
||||
Source rate during time step in [W] or [neutron/sec]
|
||||
step_ind : int
|
||||
Step index.
|
||||
proc_time : float or None
|
||||
Total process time spent depleting materials. This may
|
||||
be process-dependent and will be reduced across MPI
|
||||
processes.
|
||||
time_units : {"s", "d", "h", "min"}, optional
|
||||
Desired units for the times array
|
||||
|
||||
Returns
|
||||
-------
|
||||
times : numpy.ndarray
|
||||
Array of times in specified units
|
||||
eigenvalues : numpy.ndarray
|
||||
k-eigenvalue at each time. Column 0
|
||||
contains the eigenvalue, while column
|
||||
1 contains the associated uncertainty
|
||||
|
||||
"""
|
||||
# Get indexing terms
|
||||
vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
|
||||
cv.check_value("time_units", time_units, {"s", "d", "min", "h"})
|
||||
|
||||
stages = len(x)
|
||||
times = np.empty_like(self, dtype=float)
|
||||
eigenvalues = np.empty((len(self), 2), dtype=float)
|
||||
|
||||
# Create results
|
||||
results = Results()
|
||||
results.allocate(vol_dict, nuc_list, burn_list, full_burn_list, stages)
|
||||
# Get time/eigenvalue at each point
|
||||
for i, result in enumerate(self):
|
||||
times[i] = result.time[0]
|
||||
eigenvalues[i] = result.k[0]
|
||||
|
||||
n_mat = len(burn_list)
|
||||
# Convert time units if necessary
|
||||
times = _get_time_as(times, time_units)
|
||||
return times, eigenvalues
|
||||
|
||||
for i in range(stages):
|
||||
for mat_i in range(n_mat):
|
||||
results[i, mat_i, :] = x[i][mat_i]
|
||||
def get_eigenvalue(self, time_units='s'):
|
||||
warn("The get_eigenvalue(...) function has been renamed get_keff and "
|
||||
"will be removed in a future version of OpenMC.", FutureWarning)
|
||||
return self.get_keff(time_units)
|
||||
|
||||
results.k = [(r.k.nominal_value, r.k.std_dev) for r in op_results]
|
||||
results.rates = [r.rates for r in op_results]
|
||||
results.time = t
|
||||
results.source_rate = source_rate
|
||||
results.proc_time = proc_time
|
||||
if results.proc_time is not None:
|
||||
results.proc_time = comm.reduce(proc_time, op=MPI.SUM)
|
||||
def get_depletion_time(self):
|
||||
"""Return an array of the average time to deplete a material
|
||||
|
||||
results.export_to_hdf5("depletion_results.h5", step_ind)
|
||||
.. note::
|
||||
The return value will have one fewer values than several other
|
||||
methods, such as :meth:`get_keff`, because no depletion is performed
|
||||
at the final transport stage.
|
||||
|
||||
def transfer_volumes(self, model):
|
||||
"""Transfers volumes from depletion results to geometry
|
||||
|
||||
Parameters
|
||||
----------
|
||||
model : OpenMC model to be used in a depletion restart
|
||||
calculation
|
||||
Returns
|
||||
-------
|
||||
times : numpy.ndarray
|
||||
Vector of average time to deplete a single material
|
||||
across all processes and materials.
|
||||
|
||||
"""
|
||||
|
||||
if not model.materials:
|
||||
materials = openmc.Materials(
|
||||
model.geometry.get_all_materials().values()
|
||||
)
|
||||
times = np.empty(len(self) - 1)
|
||||
# Need special logic because the predictor
|
||||
# writes EOS values for step i as BOS values
|
||||
# for step i+1
|
||||
# The first proc_time may be zero
|
||||
if self[0].proc_time > 0.0:
|
||||
items = self[:-1]
|
||||
else:
|
||||
materials = model.materials
|
||||
items = self[1:]
|
||||
for ix, res in enumerate(items):
|
||||
times[ix] = res.proc_time
|
||||
return times
|
||||
|
||||
for material in materials:
|
||||
if material.depletable:
|
||||
material.volume = self.volume[str(material.id)]
|
||||
def get_times(self, time_units="d") -> np.ndarray:
|
||||
"""Return the points in time that define the depletion schedule
|
||||
|
||||
.. versionadded:: 0.12.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
time_units : {"s", "d", "h", "min"}, optional
|
||||
Return the vector in these units. Default is to
|
||||
convert to days
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.ndarray
|
||||
1-D vector of time points
|
||||
|
||||
"""
|
||||
cv.check_value("time_units", time_units, {"s", "d", "min", "h"})
|
||||
|
||||
times = np.fromiter(
|
||||
(r.time[0] for r in self),
|
||||
dtype=self[0].time.dtype,
|
||||
count=len(self),
|
||||
)
|
||||
|
||||
return _get_time_as(times, time_units)
|
||||
|
||||
def get_step_where(
|
||||
self, time, time_units="d", atol=1e-6, rtol=1e-3
|
||||
) -> int:
|
||||
"""Return the index closest to a given point in time
|
||||
|
||||
In the event ``time`` lies exactly between two points, the
|
||||
lower index will be returned. It is possible that the index
|
||||
will be at most one past the point in time requested, but only
|
||||
according to tolerances requested.
|
||||
|
||||
Passing ``atol=math.inf`` and ``rtol=math.inf`` will return
|
||||
the closest index to the requested point.
|
||||
|
||||
.. versionadded:: 0.12.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
time : float
|
||||
Desired point in time
|
||||
time_units : {"s", "d", "min", "h"}, optional
|
||||
Units on ``time``. Default: days
|
||||
atol : float, optional
|
||||
Absolute tolerance (in ``time_units``) if ``time`` is not
|
||||
found.
|
||||
rtol : float, optional
|
||||
Relative tolerance if ``time`` is not found.
|
||||
|
||||
Returns
|
||||
-------
|
||||
int
|
||||
|
||||
"""
|
||||
cv.check_type("time", time, numbers.Real)
|
||||
cv.check_type("atol", atol, numbers.Real)
|
||||
cv.check_type("rtol", rtol, numbers.Real)
|
||||
|
||||
times = self.get_times(time_units)
|
||||
|
||||
if times[0] < time < times[-1]:
|
||||
ix = bisect.bisect_left(times, time)
|
||||
if ix == times.size:
|
||||
ix -= 1
|
||||
# Bisection will place us either directly on the point
|
||||
# or one-past the first value less than time
|
||||
elif time - times[ix - 1] <= times[ix] - time:
|
||||
ix -= 1
|
||||
elif times[0] >= time:
|
||||
ix = 0
|
||||
elif time >= times[-1]:
|
||||
ix = times.size - 1
|
||||
|
||||
if math.isclose(time, times[ix], rel_tol=rtol, abs_tol=atol):
|
||||
return ix
|
||||
|
||||
raise ValueError(
|
||||
"A value of {} {} was not found given absolute and "
|
||||
"relative tolerances {} and {}.".format(
|
||||
time, time_units, atol, rtol)
|
||||
)
|
||||
|
||||
def export_to_materials(self, burnup_index, nuc_with_data=None) -> Materials:
|
||||
"""Return openmc.Materials object based on results at a given step
|
||||
|
||||
.. versionadded:: 0.12.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
burn_index : int
|
||||
Index of burnup step to evaluate. See also: get_step_where for
|
||||
obtaining burnup step indices from other data such as the time.
|
||||
nuc_with_data : Iterable of str, optional
|
||||
Nuclides to include in resulting materials.
|
||||
This can be specified if not all nuclides appearing in
|
||||
depletion results have associated neutron cross sections, and
|
||||
as such cannot be used in subsequent transport calculations.
|
||||
If not provided, nuclides from the cross_sections element of
|
||||
materials.xml will be used. If that element is not present,
|
||||
nuclides from OPENMC_CROSS_SECTIONS will be used.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mat_file : Materials
|
||||
A modified Materials instance containing depleted material data
|
||||
and original isotopic compositions of non-depletable materials
|
||||
"""
|
||||
result = self[burnup_index]
|
||||
|
||||
# Only materials found in the original materials.xml file will be
|
||||
# updated. If for some reason you have modified OpenMC to produce
|
||||
# new materials as depletion takes place, this method will not
|
||||
# work as expected and leave out that material.
|
||||
mat_file = Materials.from_xml("materials.xml")
|
||||
|
||||
# Only nuclides with valid transport data will be written to
|
||||
# the new materials XML file. The precedence of nuclides to select
|
||||
# is first ones provided as a kwarg here, then ones specified
|
||||
# in the materials.xml file if provided, then finally from
|
||||
# the environment variable OPENMC_CROSS_SECTIONS.
|
||||
if nuc_with_data:
|
||||
cv.check_iterable_type('nuclide names', nuc_with_data, str)
|
||||
available_cross_sections = nuc_with_data
|
||||
else:
|
||||
# select cross_sections.xml file to use
|
||||
if mat_file.cross_sections:
|
||||
this_library = DataLibrary.from_xml(path=mat_file.cross_sections)
|
||||
else:
|
||||
this_library = DataLibrary.from_xml()
|
||||
|
||||
# Find neutron libraries we have access to
|
||||
available_cross_sections = set()
|
||||
for lib in this_library.libraries:
|
||||
if lib['type'] == 'neutron':
|
||||
available_cross_sections.update(lib['materials'])
|
||||
if not available_cross_sections:
|
||||
raise DataError('No neutron libraries found in cross_sections.xml')
|
||||
|
||||
# Overwrite material definitions, if they can be found in the depletion
|
||||
# results, and save them to the new depleted xml file.
|
||||
for mat in mat_file:
|
||||
mat_id = str(mat.id)
|
||||
if mat_id in result.mat_to_ind:
|
||||
mat.volume = result.volume[mat_id]
|
||||
mat.set_density('sum')
|
||||
for nuc in result.nuc_to_ind:
|
||||
if nuc not in available_cross_sections:
|
||||
continue
|
||||
atoms = result[0, mat_id, nuc]
|
||||
if atoms > 0.0:
|
||||
atoms_per_barn_cm = 1e-24 * atoms / mat.volume
|
||||
mat.remove_nuclide(nuc) # Replace if it's there
|
||||
mat.add_nuclide(nuc, atoms_per_barn_cm)
|
||||
|
||||
return mat_file
|
||||
|
||||
|
||||
# Retain deprecated name for the time being
|
||||
ResultsList = Results
|
||||
|
|
|
|||
|
|
@ -1,396 +0,0 @@
|
|||
import numbers
|
||||
import bisect
|
||||
import math
|
||||
from warnings import warn
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
|
||||
from .results import Results, VERSION_RESULTS
|
||||
import openmc.checkvalue as cv
|
||||
from openmc.data.library import DataLibrary
|
||||
from openmc.material import Material, Materials
|
||||
from openmc.exceptions import DataError, InvalidArgumentError
|
||||
|
||||
__all__ = ["ResultsList"]
|
||||
|
||||
|
||||
def _get_time_as(seconds, units):
|
||||
if units == "d":
|
||||
return seconds / (60 * 60 * 24)
|
||||
elif units == "h":
|
||||
return seconds / (60 * 60)
|
||||
elif units == "min":
|
||||
return seconds / 60
|
||||
else:
|
||||
return seconds
|
||||
|
||||
|
||||
class ResultsList(list):
|
||||
"""A list of openmc.deplete.Results objects
|
||||
|
||||
It is recommended to use :meth:`from_hdf5` over
|
||||
direct creation.
|
||||
"""
|
||||
|
||||
@classmethod
|
||||
def from_hdf5(cls, filename):
|
||||
"""Load in depletion results from a previous file
|
||||
|
||||
Parameters
|
||||
----------
|
||||
filename : str
|
||||
Path to depletion result file
|
||||
|
||||
Returns
|
||||
-------
|
||||
new : ResultsList
|
||||
New instance of depletion results
|
||||
"""
|
||||
with h5py.File(str(filename), "r") as fh:
|
||||
cv.check_filetype_version(fh, 'depletion results', VERSION_RESULTS[0])
|
||||
new = cls()
|
||||
|
||||
# Get number of results stored
|
||||
n = fh["number"][...].shape[0]
|
||||
|
||||
for i in range(n):
|
||||
new.append(Results.from_hdf5(fh, i))
|
||||
return new
|
||||
|
||||
def get_atoms(self, mat, nuc, nuc_units="atoms", time_units="s"):
|
||||
"""Get number of nuclides over time from a single material
|
||||
|
||||
.. note::
|
||||
Initial values for some isotopes that do not appear in
|
||||
initial concentrations may be non-zero, depending on the
|
||||
value of :class:`openmc.deplete.Operator` ``dilute_initial``.
|
||||
The :class:`openmc.deplete.Operator` adds isotopes according
|
||||
to this setting, which can be set to zero.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mat : openmc.Material, str
|
||||
Material object or material id to evaluate
|
||||
nuc : str
|
||||
Nuclide name to evaluate
|
||||
nuc_units : {"atoms", "atom/b-cm", "atom/cm3"}, optional
|
||||
Units for the returned concentration. Default is ``"atoms"``
|
||||
|
||||
.. versionadded:: 0.12
|
||||
time_units : {"s", "min", "h", "d"}, optional
|
||||
Units for the returned time array. Default is ``"s"`` to
|
||||
return the value in seconds.
|
||||
|
||||
.. versionadded:: 0.12
|
||||
|
||||
Returns
|
||||
-------
|
||||
times : numpy.ndarray
|
||||
Array of times in units of ``time_units``
|
||||
concentrations : numpy.ndarray
|
||||
Concentration of specified nuclide in units of ``nuc_units``
|
||||
|
||||
"""
|
||||
cv.check_value("time_units", time_units, {"s", "d", "min", "h"})
|
||||
cv.check_value("nuc_units", nuc_units,
|
||||
{"atoms", "atom/b-cm", "atom/cm3"})
|
||||
|
||||
if isinstance(mat, Material):
|
||||
mat_id = str(mat.id)
|
||||
elif isinstance(mat, str):
|
||||
mat_id = mat
|
||||
else:
|
||||
raise TypeError('mat should be of type openmc.Material or str')
|
||||
times = np.empty_like(self, dtype=float)
|
||||
concentrations = np.empty_like(self, dtype=float)
|
||||
|
||||
# Evaluate value in each region
|
||||
for i, result in enumerate(self):
|
||||
times[i] = result.time[0]
|
||||
concentrations[i] = result[0, mat_id, nuc]
|
||||
|
||||
# Unit conversions
|
||||
times = _get_time_as(times, time_units)
|
||||
if nuc_units != "atoms":
|
||||
# Divide by volume to get density
|
||||
concentrations /= self[0].volume[mat_id]
|
||||
if nuc_units == "atom/b-cm":
|
||||
# 1 barn = 1e-24 cm^2
|
||||
concentrations *= 1e-24
|
||||
|
||||
return times, concentrations
|
||||
|
||||
def get_reaction_rate(self, mat, nuc, rx):
|
||||
"""Get reaction rate in a single material/nuclide over time
|
||||
|
||||
.. note::
|
||||
|
||||
Initial values for some isotopes that do not appear in
|
||||
initial concentrations may be non-zero, depending on the
|
||||
value of :class:`openmc.deplete.Operator` ``dilute_initial``
|
||||
The :class:`openmc.deplete.Operator` adds isotopes according
|
||||
to this setting, which can be set to zero.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
mat : openmc.Material, str
|
||||
Material object or material id to evaluate
|
||||
nuc : str
|
||||
Nuclide name to evaluate
|
||||
rx : str
|
||||
Reaction rate to evaluate
|
||||
|
||||
Returns
|
||||
-------
|
||||
times : numpy.ndarray
|
||||
Array of times in [s]
|
||||
rates : numpy.ndarray
|
||||
Array of reaction rates
|
||||
|
||||
"""
|
||||
times = np.empty_like(self, dtype=float)
|
||||
rates = np.empty_like(self, dtype=float)
|
||||
|
||||
if isinstance(mat, Material):
|
||||
mat_id = str(mat.id)
|
||||
elif isinstance(mat, str):
|
||||
mat_id = mat
|
||||
else:
|
||||
raise TypeError('mat should be of type openmc.Material or str')
|
||||
|
||||
# Evaluate value in each region
|
||||
for i, result in enumerate(self):
|
||||
times[i] = result.time[0]
|
||||
rates[i] = result.rates[0].get(mat_id, nuc, rx) * result[0, mat, nuc]
|
||||
|
||||
return times, rates
|
||||
|
||||
def get_keff(self, time_units='s'):
|
||||
"""Evaluates the eigenvalue from a results list.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
time_units : {"s", "d", "h", "min"}, optional
|
||||
Desired units for the times array
|
||||
|
||||
Returns
|
||||
-------
|
||||
times : numpy.ndarray
|
||||
Array of times in specified units
|
||||
eigenvalues : numpy.ndarray
|
||||
k-eigenvalue at each time. Column 0
|
||||
contains the eigenvalue, while column
|
||||
1 contains the associated uncertainty
|
||||
|
||||
"""
|
||||
cv.check_value("time_units", time_units, {"s", "d", "min", "h"})
|
||||
|
||||
times = np.empty_like(self, dtype=float)
|
||||
eigenvalues = np.empty((len(self), 2), dtype=float)
|
||||
|
||||
# Get time/eigenvalue at each point
|
||||
for i, result in enumerate(self):
|
||||
times[i] = result.time[0]
|
||||
eigenvalues[i] = result.k[0]
|
||||
|
||||
# Convert time units if necessary
|
||||
times = _get_time_as(times, time_units)
|
||||
return times, eigenvalues
|
||||
|
||||
def get_eigenvalue(self, time_units='s'):
|
||||
warn("The get_eigenvalue(...) function has been renamed get_keff and "
|
||||
"will be removed in a future version of OpenMC.", FutureWarning)
|
||||
return self.get_keff(time_units)
|
||||
|
||||
|
||||
def get_depletion_time(self):
|
||||
"""Return an array of the average time to deplete a material
|
||||
|
||||
.. note::
|
||||
The return value will have one fewer values than several other
|
||||
methods, such as :meth:`get_keff`, because no depletion is performed
|
||||
at the final transport stage.
|
||||
|
||||
Returns
|
||||
-------
|
||||
times : numpy.ndarray
|
||||
Vector of average time to deplete a single material
|
||||
across all processes and materials.
|
||||
|
||||
"""
|
||||
times = np.empty(len(self) - 1)
|
||||
# Need special logic because the predictor
|
||||
# writes EOS values for step i as BOS values
|
||||
# for step i+1
|
||||
# The first proc_time may be zero
|
||||
if self[0].proc_time > 0.0:
|
||||
items = self[:-1]
|
||||
else:
|
||||
items = self[1:]
|
||||
for ix, res in enumerate(items):
|
||||
times[ix] = res.proc_time
|
||||
return times
|
||||
|
||||
def get_times(self, time_units="d") -> np.ndarray:
|
||||
"""Return the points in time that define the depletion schedule
|
||||
|
||||
|
||||
.. versionadded:: 0.12.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
time_units : {"s", "d", "h", "min"}, optional
|
||||
Return the vector in these units. Default is to
|
||||
convert to days
|
||||
|
||||
Returns
|
||||
-------
|
||||
numpy.ndarray
|
||||
1-D vector of time points
|
||||
|
||||
"""
|
||||
cv.check_value("time_units", time_units, {"s", "d", "min", "h"})
|
||||
|
||||
times = np.fromiter(
|
||||
(r.time[0] for r in self),
|
||||
dtype=self[0].time.dtype,
|
||||
count=len(self),
|
||||
)
|
||||
|
||||
return _get_time_as(times, time_units)
|
||||
|
||||
def get_step_where(
|
||||
self, time, time_units="d", atol=1e-6, rtol=1e-3
|
||||
) -> int:
|
||||
"""Return the index closest to a given point in time
|
||||
|
||||
In the event ``time`` lies exactly between two points, the
|
||||
lower index will be returned. It is possible that the index
|
||||
will be at most one past the point in time requested, but only
|
||||
according to tolerances requested.
|
||||
|
||||
Passing ``atol=math.inf`` and ``rtol=math.inf`` will return
|
||||
the closest index to the requested point.
|
||||
|
||||
|
||||
.. versionadded:: 0.12.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
time : float
|
||||
Desired point in time
|
||||
time_units : {"s", "d", "min", "h"}, optional
|
||||
Units on ``time``. Default: days
|
||||
atol : float, optional
|
||||
Absolute tolerance (in ``time_units``) if ``time`` is not
|
||||
found.
|
||||
rtol : float, optional
|
||||
Relative tolerance if ``time`` is not found.
|
||||
|
||||
Returns
|
||||
-------
|
||||
int
|
||||
|
||||
"""
|
||||
cv.check_type("time", time, numbers.Real)
|
||||
cv.check_type("atol", atol, numbers.Real)
|
||||
cv.check_type("rtol", rtol, numbers.Real)
|
||||
|
||||
times = self.get_times(time_units)
|
||||
|
||||
if times[0] < time < times[-1]:
|
||||
ix = bisect.bisect_left(times, time)
|
||||
if ix == times.size:
|
||||
ix -= 1
|
||||
# Bisection will place us either directly on the point
|
||||
# or one-past the first value less than time
|
||||
elif time - times[ix - 1] <= times[ix] - time:
|
||||
ix -= 1
|
||||
elif times[0] >= time:
|
||||
ix = 0
|
||||
elif time >= times[-1]:
|
||||
ix = times.size - 1
|
||||
|
||||
if math.isclose(time, times[ix], rel_tol=rtol, abs_tol=atol):
|
||||
return ix
|
||||
|
||||
raise ValueError(
|
||||
"A value of {} {} was not found given absolute and "
|
||||
"relative tolerances {} and {}.".format(
|
||||
time, time_units, atol, rtol)
|
||||
)
|
||||
|
||||
def export_to_materials(self, burnup_index, nuc_with_data=None) -> Materials:
|
||||
"""Return openmc.Materials object based on results at a given step
|
||||
|
||||
.. versionadded:: 0.12.1
|
||||
|
||||
Parameters
|
||||
----------
|
||||
burn_index : int
|
||||
Index of burnup step to evaluate. See also: get_step_where for
|
||||
obtaining burnup step indices from other data such as the time.
|
||||
nuc_with_data : Iterable of str, optional
|
||||
Nuclides to include in resulting materials.
|
||||
This can be specified if not all nuclides appearing in
|
||||
depletion results have associated neutron cross sections, and
|
||||
as such cannot be used in subsequent transport calculations.
|
||||
If not provided, nuclides from the cross_sections element of
|
||||
materials.xml will be used. If that element is not present,
|
||||
nuclides from OPENMC_CROSS_SECTIONS will be used.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mat_file : Materials
|
||||
A modified Materials instance containing depleted material data
|
||||
and original isotopic compositions of non-depletable materials
|
||||
"""
|
||||
result = self[burnup_index]
|
||||
|
||||
# Only materials found in the original materials.xml file will be
|
||||
# updated. If for some reason you have modified OpenMC to produce
|
||||
# new materials as depletion takes place, this method will not
|
||||
# work as expected and leave out that material.
|
||||
mat_file = Materials.from_xml("materials.xml")
|
||||
|
||||
# Only nuclides with valid transport data will be written to
|
||||
# the new materials XML file. The precedence of nuclides to select
|
||||
# is first ones provided as a kwarg here, then ones specified
|
||||
# in the materials.xml file if provided, then finally from
|
||||
# the environment variable OPENMC_CROSS_SECTIONS.
|
||||
if nuc_with_data:
|
||||
cv.check_iterable_type('nuclide names', nuc_with_data, str)
|
||||
available_cross_sections = nuc_with_data
|
||||
else:
|
||||
# select cross_sections.xml file to use
|
||||
if mat_file.cross_sections:
|
||||
this_library = DataLibrary.from_xml(path=mat_file.cross_sections)
|
||||
else:
|
||||
this_library = DataLibrary.from_xml()
|
||||
|
||||
# Find neutron libraries we have access to
|
||||
available_cross_sections = set()
|
||||
for lib in this_library.libraries:
|
||||
if lib['type'] == 'neutron':
|
||||
available_cross_sections.update(lib['materials'])
|
||||
if not available_cross_sections:
|
||||
raise DataError('No neutron libraries found in cross_sections.xml')
|
||||
|
||||
# Overwrite material definitions, if they can be found in the depletion
|
||||
# results, and save them to the new depleted xml file.
|
||||
for mat in mat_file:
|
||||
mat_id = str(mat.id)
|
||||
if mat_id in result.mat_to_ind:
|
||||
mat.volume = result.volume[mat_id]
|
||||
mat.set_density('sum')
|
||||
for nuc in result.nuc_to_ind:
|
||||
if nuc not in available_cross_sections:
|
||||
continue
|
||||
atoms = result[0, mat_id, nuc]
|
||||
if atoms > 0.0:
|
||||
atoms_per_barn_cm = 1e-24 * atoms / mat.volume
|
||||
mat.remove_nuclide(nuc) # Replace if it's there
|
||||
mat.add_nuclide(nuc, atoms_per_barn_cm)
|
||||
|
||||
return mat_file
|
||||
527
openmc/deplete/stepresult.py
Normal file
527
openmc/deplete/stepresult.py
Normal file
|
|
@ -0,0 +1,527 @@
|
|||
"""The stepresult module.
|
||||
|
||||
Contains capabilities for generating and saving results of a single depletion
|
||||
timestep.
|
||||
"""
|
||||
|
||||
from collections import OrderedDict
|
||||
import copy
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
|
||||
import openmc
|
||||
from openmc.mpi import comm, MPI
|
||||
from .reaction_rates import ReactionRates
|
||||
|
||||
VERSION_RESULTS = (1, 1)
|
||||
|
||||
|
||||
__all__ = ["StepResult"]
|
||||
|
||||
|
||||
class StepResult:
|
||||
"""Result of a single depletion timestep
|
||||
|
||||
Attributes
|
||||
----------
|
||||
k : list of (float, float)
|
||||
Eigenvalue and uncertainty for each substep.
|
||||
time : list of float
|
||||
Time at beginning, end of step, in seconds.
|
||||
source_rate : float
|
||||
Source rate during timestep in [W] or [neutron/sec]
|
||||
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 str 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.
|
||||
proc_time : int
|
||||
Average time spent depleting a material across all
|
||||
materials and processes
|
||||
|
||||
"""
|
||||
def __init__(self):
|
||||
self.k = None
|
||||
self.time = None
|
||||
self.source_rate = None
|
||||
self.rates = None
|
||||
self.volume = None
|
||||
self.proc_time = None
|
||||
|
||||
self.mat_to_ind = None
|
||||
self.nuc_to_ind = None
|
||||
self.mat_to_hdf5_ind = None
|
||||
|
||||
self.data = None
|
||||
|
||||
def __repr__(self):
|
||||
t = self.time[0]
|
||||
dt = self.time[1] - self.time[0]
|
||||
return f"<StepResult: t={t}, dt={dt}, source={self.source_rate}>"
|
||||
|
||||
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, openmc.Material):
|
||||
mat = str(mat.id)
|
||||
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):
|
||||
"""Allocate memory for depletion step data
|
||||
|
||||
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 distribute(self, local_materials, ranges):
|
||||
"""Create a new object containing data for distributed materials
|
||||
|
||||
Parameters
|
||||
----------
|
||||
local_materials : iterable of str
|
||||
Materials for this process
|
||||
ranges : iterable of int
|
||||
Slice-like object indicating indicies of ``local_materials``
|
||||
in the material dimension of :attr:`data` and each element
|
||||
in :attr:`rates`
|
||||
|
||||
Returns
|
||||
-------
|
||||
StepResult
|
||||
New results object
|
||||
"""
|
||||
new = StepResult()
|
||||
new.volume = {lm: self.volume[lm] for lm in local_materials}
|
||||
new.mat_to_ind = {mat: idx for (idx, mat) in enumerate(local_materials)}
|
||||
|
||||
# Direct transfer
|
||||
direct_attrs = ("time", "k", "source_rate", "nuc_to_ind",
|
||||
"mat_to_hdf5_ind", "proc_time")
|
||||
for attr in direct_attrs:
|
||||
setattr(new, attr, getattr(self, attr))
|
||||
# Get applicable slice of data
|
||||
new.data = self.data[:, ranges]
|
||||
new.rates = [r[ranges] for r in self.rates]
|
||||
return new
|
||||
|
||||
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?
|
||||
|
||||
"""
|
||||
# Write new file if first time step, else add to existing file
|
||||
kwargs = {'mode': "w" if step == 0 else "a"}
|
||||
|
||||
if h5py.get_config().mpi and comm.size > 1:
|
||||
# Write results in parallel
|
||||
kwargs['driver'] = 'mpio'
|
||||
kwargs['comm'] = comm
|
||||
with h5py.File(filename, **kwargs) as handle:
|
||||
self._to_hdf5(handle, step, parallel=True)
|
||||
else:
|
||||
# Gather results at root process
|
||||
all_results = comm.gather(self)
|
||||
|
||||
# Only root process writes results
|
||||
if comm.rank == 0:
|
||||
with h5py.File(filename, **kwargs) as handle:
|
||||
for res in all_results:
|
||||
res._to_hdf5(handle, step, parallel=False)
|
||||
|
||||
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, 2),
|
||||
maxshape=(None, n_stages, 2), dtype='float64')
|
||||
|
||||
handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
|
||||
|
||||
handle.create_dataset("source_rate", (1, n_stages), maxshape=(None, n_stages),
|
||||
dtype='float64')
|
||||
|
||||
handle.create_dataset(
|
||||
"depletion time", (1,), maxshape=(None,),
|
||||
dtype="float64")
|
||||
|
||||
def _to_hdf5(self, handle, index, parallel=False):
|
||||
"""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?
|
||||
parallel : bool
|
||||
Being called with parallel HDF5?
|
||||
|
||||
"""
|
||||
if "/number" not in handle:
|
||||
if parallel:
|
||||
comm.barrier()
|
||||
self._write_hdf5_metadata(handle)
|
||||
|
||||
if parallel:
|
||||
comm.barrier()
|
||||
|
||||
# Grab handles
|
||||
number_dset = handle["/number"]
|
||||
rxn_dset = handle["/reaction rates"]
|
||||
eigenvalues_dset = handle["/eigenvalues"]
|
||||
time_dset = handle["/time"]
|
||||
source_rate_dset = handle["/source_rate"]
|
||||
proc_time_dset = handle["/depletion 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)
|
||||
|
||||
source_rate_shape = list(source_rate_dset.shape)
|
||||
source_rate_shape[0] = new_shape
|
||||
source_rate_dset.resize(source_rate_shape)
|
||||
|
||||
proc_shape = list(proc_time_dset.shape)
|
||||
proc_shape[0] = new_shape
|
||||
proc_time_dset.resize(proc_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
|
||||
source_rate_dset[index] = self.source_rate
|
||||
if self.proc_time is not None:
|
||||
proc_time_dset[index] = (
|
||||
self.proc_time / (comm.size * self.n_hdf5_mats)
|
||||
)
|
||||
|
||||
@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
|
||||
Index for depletion step
|
||||
"""
|
||||
results = cls()
|
||||
|
||||
# Grab handles
|
||||
number_dset = handle["/number"]
|
||||
eigenvalues_dset = handle["/eigenvalues"]
|
||||
time_dset = handle["/time"]
|
||||
if "source_rate" in handle:
|
||||
source_rate_dset = handle["/source_rate"]
|
||||
else:
|
||||
# Older versions used "power" instead of "source_rate"
|
||||
source_rate_dset = handle["/power"]
|
||||
|
||||
results.data = number_dset[step, :, :, :]
|
||||
results.k = eigenvalues_dset[step, :]
|
||||
results.time = time_dset[step, :]
|
||||
results.source_rate = source_rate_dset[step, 0]
|
||||
|
||||
if "depletion time" in handle:
|
||||
proc_time_dset = handle["/depletion time"]
|
||||
if step < proc_time_dset.shape[0]:
|
||||
results.proc_time = proc_time_dset[step]
|
||||
|
||||
if results.proc_time is None:
|
||||
results.proc_time = np.array([np.nan])
|
||||
|
||||
# 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, True)
|
||||
|
||||
rate[:] = handle["/reaction rates"][step, i, :, :, :]
|
||||
results.rates.append(rate)
|
||||
|
||||
return results
|
||||
|
||||
@staticmethod
|
||||
def save(op, x, op_results, t, source_rate, step_ind, proc_time=None):
|
||||
"""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.
|
||||
source_rate : float
|
||||
Source rate during time step in [W] or [neutron/sec]
|
||||
step_ind : int
|
||||
Step index.
|
||||
proc_time : float or None
|
||||
Total process time spent depleting materials. This may
|
||||
be process-dependent and will be reduced across MPI
|
||||
processes.
|
||||
|
||||
"""
|
||||
# Get indexing terms
|
||||
vol_dict, nuc_list, burn_list, full_burn_list = op.get_results_info()
|
||||
|
||||
stages = len(x)
|
||||
|
||||
# Create results
|
||||
results = StepResult()
|
||||
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.nominal_value, r.k.std_dev) for r in op_results]
|
||||
results.rates = [r.rates for r in op_results]
|
||||
results.time = t
|
||||
results.source_rate = source_rate
|
||||
results.proc_time = proc_time
|
||||
if results.proc_time is not None:
|
||||
results.proc_time = comm.reduce(proc_time, op=MPI.SUM)
|
||||
|
||||
results.export_to_hdf5("depletion_results.h5", step_ind)
|
||||
|
||||
def transfer_volumes(self, model):
|
||||
"""Transfers volumes from depletion results to geometry
|
||||
|
||||
Parameters
|
||||
----------
|
||||
model : OpenMC model to be used in a depletion restart
|
||||
calculation
|
||||
|
||||
"""
|
||||
|
||||
if not model.materials:
|
||||
materials = openmc.Materials(
|
||||
model.geometry.get_all_materials().values()
|
||||
)
|
||||
else:
|
||||
materials = model.materials
|
||||
|
||||
for material in materials:
|
||||
if material.depletable:
|
||||
material.volume = self.volume[str(material.id)]
|
||||
|
|
@ -77,8 +77,8 @@ def test_full(run_in_tmpdir, problem, multiproc):
|
|||
return
|
||||
|
||||
# Load the reference/test results
|
||||
res_test = openmc.deplete.ResultsList.from_hdf5(path_test)
|
||||
res_ref = openmc.deplete.ResultsList.from_hdf5(path_reference)
|
||||
res_test = openmc.deplete.Results(path_test)
|
||||
res_ref = openmc.deplete.Results(path_reference)
|
||||
|
||||
# Assert same mats
|
||||
for mat in res_ref[0].mat_to_ind:
|
||||
|
|
@ -139,7 +139,7 @@ def test_depletion_results_to_material(run_in_tmpdir, problem):
|
|||
"""Checks openmc.Materials objects can be created from depletion results"""
|
||||
# Load the reference/test results
|
||||
path_reference = Path(__file__).with_name('test_reference.h5')
|
||||
res_ref = openmc.deplete.ResultsList.from_hdf5(path_reference)
|
||||
res_ref = openmc.deplete.Results(path_reference)
|
||||
|
||||
# Firstly need to export materials.xml file for the initial simulation state
|
||||
geometry, lower_left, upper_right = problem
|
||||
|
|
|
|||
|
|
@ -107,7 +107,7 @@ def test_activation(run_in_tmpdir, model, reaction_rate_mode, reaction_rate_opts
|
|||
integrator.integrate()
|
||||
|
||||
# Get resulting number of atoms
|
||||
results = openmc.deplete.ResultsList.from_hdf5('depletion_results.h5')
|
||||
results = openmc.deplete.Results('depletion_results.h5')
|
||||
_, atoms = results.get_atoms(w, "W186")
|
||||
|
||||
assert atoms[0] == pytest.approx(n0)
|
||||
|
|
@ -155,7 +155,7 @@ def test_decay(run_in_tmpdir):
|
|||
integrator.integrate()
|
||||
|
||||
# Get resulting number of atoms
|
||||
results = openmc.deplete.ResultsList.from_hdf5('depletion_results.h5')
|
||||
results = openmc.deplete.Results('depletion_results.h5')
|
||||
_, atoms = results.get_atoms(mat, "Sr89")
|
||||
|
||||
# Ensure density goes down by a factor of 2 after each half-life
|
||||
|
|
|
|||
|
|
@ -16,7 +16,7 @@ import pytest
|
|||
|
||||
from openmc.mpi import comm
|
||||
from openmc.deplete import (
|
||||
ReactionRates, Results, ResultsList, OperatorResult, PredictorIntegrator,
|
||||
ReactionRates, StepResult, Results, OperatorResult, PredictorIntegrator,
|
||||
CECMIntegrator, CF4Integrator, CELIIntegrator, EPCRK4Integrator,
|
||||
LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator, cram)
|
||||
|
||||
|
|
@ -99,11 +99,11 @@ def test_results_save(run_in_tmpdir):
|
|||
for k, rates in zip(eigvl1, rate1)]
|
||||
op_result2 = [OperatorResult(ufloat(*k), rates)
|
||||
for k, rates in zip(eigvl2, rate2)]
|
||||
Results.save(op, x1, op_result1, t1, 0, 0)
|
||||
Results.save(op, x2, op_result2, t2, 0, 1)
|
||||
StepResult.save(op, x1, op_result1, t1, 0, 0)
|
||||
StepResult.save(op, x2, op_result2, t2, 0, 1)
|
||||
|
||||
# Load the files
|
||||
res = ResultsList.from_hdf5("depletion_results.h5")
|
||||
res = Results("depletion_results.h5")
|
||||
|
||||
for i in range(stages):
|
||||
for mat_i, mat in enumerate(burn_list):
|
||||
|
|
@ -176,8 +176,7 @@ def test_integrator(run_in_tmpdir, scheme):
|
|||
|
||||
# get expected results
|
||||
|
||||
res = ResultsList.from_hdf5(
|
||||
operator.output_dir / "depletion_results.h5")
|
||||
res = Results(operator.output_dir / "depletion_results.h5")
|
||||
|
||||
t1, y1 = res.get_atoms("1", "1")
|
||||
t2, y2 = res.get_atoms("1", "2")
|
||||
|
|
|
|||
|
|
@ -24,8 +24,7 @@ def test_restart_predictor_cecm(run_in_tmpdir):
|
|||
openmc.deplete.PredictorIntegrator(op, dt, power).integrate()
|
||||
|
||||
# Load the files
|
||||
prev_res = openmc.deplete.ResultsList.from_hdf5(
|
||||
op.output_dir / "depletion_results.h5")
|
||||
prev_res = openmc.deplete.Results(op.output_dir / "depletion_results.h5")
|
||||
|
||||
# Re-create depletion operator and load previous results
|
||||
op = dummy_operator.DummyOperator(prev_res)
|
||||
|
|
@ -51,8 +50,7 @@ def test_restart_cecm_predictor(run_in_tmpdir):
|
|||
cecm.integrate()
|
||||
|
||||
# Load the files
|
||||
prev_res = openmc.deplete.ResultsList.from_hdf5(
|
||||
op.output_dir / "depletion_results.h5")
|
||||
prev_res = openmc.deplete.Results(op.output_dir / "depletion_results.h5")
|
||||
|
||||
# Re-create depletion operator and load previous results
|
||||
op = dummy_operator.DummyOperator(prev_res)
|
||||
|
|
@ -75,7 +73,7 @@ def test_restart(run_in_tmpdir, scheme):
|
|||
bundle.solver(operator, [0.75], 1.0).integrate()
|
||||
|
||||
# restart
|
||||
prev_res = openmc.deplete.ResultsList.from_hdf5(
|
||||
prev_res = openmc.deplete.Results(
|
||||
operator.output_dir / "depletion_results.h5")
|
||||
operator = dummy_operator.DummyOperator(prev_res)
|
||||
|
||||
|
|
@ -84,7 +82,7 @@ def test_restart(run_in_tmpdir, scheme):
|
|||
|
||||
# compare results
|
||||
|
||||
results = openmc.deplete.ResultsList.from_hdf5(
|
||||
results = openmc.deplete.Results(
|
||||
operator.output_dir / "depletion_results.h5")
|
||||
|
||||
_t, y1 = results.get_atoms("1", "1")
|
||||
|
|
|
|||
|
|
@ -1,4 +1,4 @@
|
|||
"""Tests the ResultsList class"""
|
||||
"""Tests the Results class"""
|
||||
|
||||
from pathlib import Path
|
||||
from math import inf
|
||||
|
|
@ -13,7 +13,7 @@ def res():
|
|||
"""Load the reference results"""
|
||||
filename = (Path(__file__).parents[1] / 'regression_tests' / 'deplete'
|
||||
/ 'test_reference.h5')
|
||||
return openmc.deplete.ResultsList.from_hdf5(filename)
|
||||
return openmc.deplete.Results(filename)
|
||||
|
||||
|
||||
def test_get_atoms(res):
|
||||
|
|
@ -72,9 +72,9 @@ def test_get_keff(res):
|
|||
|
||||
@pytest.mark.parametrize("unit", ("s", "d", "min", "h"))
|
||||
def test_get_steps(unit):
|
||||
# Make a ResultsList full of near-empty Result instances
|
||||
# Make a Results full of near-empty Result instances
|
||||
# Just fill out a time schedule
|
||||
results = openmc.deplete.ResultsList()
|
||||
results = openmc.deplete.Results()
|
||||
# Time in units of unit
|
||||
times = np.linspace(0, 100, num=5)
|
||||
if unit == "d":
|
||||
|
|
@ -87,7 +87,7 @@ def test_get_steps(unit):
|
|||
conversion_to_seconds = 1
|
||||
|
||||
for ix in range(times.size):
|
||||
res = openmc.deplete.Results()
|
||||
res = openmc.deplete.StepResult()
|
||||
res.time = times[ix:ix + 1] * conversion_to_seconds
|
||||
results.append(res)
|
||||
|
||||
|
|
|
|||
|
|
@ -2,7 +2,7 @@
|
|||
|
||||
from pytest import approx
|
||||
import openmc
|
||||
from openmc.deplete import PredictorIntegrator, ResultsList
|
||||
from openmc.deplete import PredictorIntegrator, Results
|
||||
|
||||
from tests import dummy_operator
|
||||
|
||||
|
|
@ -19,7 +19,7 @@ def test_transfer_volumes(run_in_tmpdir):
|
|||
PredictorIntegrator(op, dt, power).integrate()
|
||||
|
||||
# Load the files
|
||||
res = openmc.deplete.ResultsList.from_hdf5(op.output_dir / "depletion_results.h5")
|
||||
res = openmc.deplete.Results(op.output_dir / "depletion_results.h5")
|
||||
|
||||
# Create a dictionary of volumes to transfer
|
||||
res[0].volume['1'] = 1.5
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue