crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Don't write reaction rates in depletion results by default, remove per-stage data for multistage integrators (#3609)

Browse source
This commit is contained in:
Paul Romano 2025-11-19 11:26:57 -06:00 • committed by GitHub
parent 028f440448
commit f544d02e49
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 336 additions and 273 deletions
Download patch file Download diff file Expand all files Collapse all files

26
docs/source/io_formats/depletion_results.rst
View file

@ -4,7 +4,7 @@
Depletion Results File Format
=============================
The current version of the depletion results file format is 1.1.
The current version of the depletion results file format is 1.2.
**/**
@ -12,22 +12,20 @@ The current version of the depletion results file format is 1.1.
- **version** (*int[2]*) -- Major and minor version of the
statepoint file format.
:Datasets: - **eigenvalues** (*double[][][2]*) -- k-eigenvalues at each
time/stage. This array has shape (number of timesteps, number of
stages, value). The last axis contains the eigenvalue and the
associated uncertainty
- **number** (*double[][][][]*) -- Total number of atoms. This array
has shape (number of timesteps, number of stages, number of
:Datasets: - **eigenvalues** (*double[][2]*) -- k-eigenvalues at each timestep.
This array has shape (number of timesteps, 2). The second axis
contains the eigenvalue and its associated uncertainty.
- **number** (*double[][][]*) -- Total number of atoms at each
timestep. This array has shape (number of timesteps, 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).
- **reaction rates** (*double[][][][]*) -- Reaction rates at each
timestep. This array has shape (number of timesteps, number of
materials, number of nuclides, number of reactions). Only stored if
write_rates=True.
- **time** (*double[][2]*) -- Time in [s] at beginning/end of each
step.
- **source_rate** (*double[][]*) -- Power in [W] or source rate in
[neutron/sec]. This array has shape (number of timesteps, number
of stages).
- **source_rate** (*double[]*) -- Power in [W] or source rate in
[neutron/sec] for each timestep.
- **depletion time** (*double[]*) -- Average process time in [s]
spent depleting a material across all burnable materials and,
if applicable, MPI processes.

153
openmc/deplete/abc.py
View file

@ -631,17 +631,7 @@ class Integrator(ABC):
solver: str = "cram48",
continue_timesteps: bool = False,
):
# Check number of stages previously used
if operator.prev_res is not None:
res = operator.prev_res[-1]
if res.data.shape[0] != self._num_stages:
raise ValueError(
"{} incompatible with previous restart calculation. "
"Previous scheme used {} intermediate solutions, while "
"this uses {}".format(
self.__class__.__name__, res.data.shape[0],
self._num_stages))
elif continue_timesteps:
if continue_timesteps and operator.prev_res is None:
raise ValueError("Continuation run requires passing prev_results.")
self.operator = operator
self.chain = operator.chain
@ -775,12 +765,8 @@ class Integrator(ABC):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
n_end : list of numpy.ndarray
Concentrations at end of timestep
"""
@property
@ -811,9 +797,9 @@ class Integrator(ABC):
"""Get beginning of step concentrations, reaction rates from restart"""
res = self.operator.prev_res[-1]
# Depletion methods expect list of arrays
bos_conc = list(res.data[0])
rates = res.rates[0]
k = ufloat(res.k[0, 0], res.k[0, 1])
bos_conc = list(res.data)
rates = res.rates
k = ufloat(res.k[0], res.k[1])
if res.source_rate != 0.0:
# Scale reaction rates by ratio of source rates
@ -855,7 +841,8 @@ class Integrator(ABC):
self,
final_step: bool = True,
output: bool = True,
path: PathLike = 'depletion_results.h5'
path: PathLike = 'depletion_results.h5',
write_rates: bool = False
):
"""Perform the entire depletion process across all steps
@ -874,6 +861,11 @@ class Integrator(ABC):
Path to file to write. Defaults to 'depletion_results.h5'.
.. versionadded:: 0.15.0
write_rates : bool, optional
Whether reaction rates should be written to the results file for
each step. Defaults to ``False`` to reduce file size.
.. versionadded:: 0.15.3
"""
with change_directory(self.operator.output_dir):
n = self.operator.initial_condition()
@ -890,18 +882,22 @@ class Integrator(ABC):
n, res = self._get_bos_data_from_restart(source_rate, n)
# Solve Bateman equations over time interval
proc_time, n_list, res_list = self(n, res.rates, dt, source_rate, i)
proc_time, n_end = self(n, res.rates, dt, source_rate, i)
# Insert BOS concentration, transport results
n_list.insert(0, n)
res_list.insert(0, res)
# Remove actual EOS concentration for next step
n = n_list.pop()
StepResult.save(self.operator, n_list, res_list, [t, t + dt],
source_rate, self._i_res + i, proc_time, path)
StepResult.save(
self.operator,
n,
res,
[t, t + dt],
source_rate,
self._i_res + i,
proc_time,
write_rates=write_rates,
path=path
)
# Update for next step
n = n_end
t += dt
# Final simulation -- in the case that final_step is False, a zero
@ -910,9 +906,18 @@ class Integrator(ABC):
# solve)
if output and final_step and comm.rank == 0:
print(f"[openmc.deplete] t={t} (final operator evaluation)")
res_list = [self.operator(n, source_rate if final_step else 0.0)]
StepResult.save(self.operator, [n], res_list, [t, t],
source_rate, self._i_res + len(self), proc_time, path)
res_final = self.operator(n, source_rate if final_step else 0.0)
StepResult.save(
self.operator,
n,
res_final,
[t, t],
source_rate,
self._i_res + len(self),
proc_time,
write_rates=write_rates,
path=path
)
self.operator.write_bos_data(len(self) + self._i_res)
self.operator.finalize()
@ -1171,10 +1176,40 @@ class SIIntegrator(Integrator):
self.operator.settings.particles //= self.n_steps
return inherited
@abstractmethod
def __call__(self, n, rates, dt, source_rate, i):
"""Perform the integration across one time step
Parameters
----------
n : list of numpy.ndarray
List of atom number arrays for each material. Each array has
shape ``(n_nucs,)`` where ``n_nucs`` is the number of nuclides
rates : openmc.deplete.ReactionRates
Reaction rates (from transport operator)
dt : float
Time step in [s]
source_rate : float
Power in [W] or source rate in [neutron/sec]
i : int
Current time step index
Returns
-------
proc_time : float
Time spent in transport simulation
n_end : list of numpy.ndarray
Updated atom number densities for each material
op_result : OperatorResult
Eigenvalue and reaction rates resulting from transport simulation
"""
def integrate(
self,
output: bool = True,
path: PathLike = "depletion_results.h5"
path: PathLike = "depletion_results.h5",
write_rates: bool = False
):
"""Perform the entire depletion process across all steps
@ -1186,11 +1221,17 @@ class SIIntegrator(Integrator):
Path to file to write. Defaults to 'depletion_results.h5'.
.. versionadded:: 0.15.0
write_rates : bool, optional
Whether reaction rates should be written to the results file for
each step. Defaults to ``False`` to reduce file size.
.. versionadded:: 0.15.3
"""
with change_directory(self.operator.output_dir):
n = self.operator.initial_condition()
t, self._i_res = self._get_start_data()
res_end = None # Will be set in first iteration
for i, (dt, p) in enumerate(self):
if output:
print(f"[openmc.deplete] t={t} s, dt={dt} s, source={p}")
@ -1200,28 +1241,38 @@ class SIIntegrator(Integrator):
n, res = self._get_bos_data_from_operator(i, p, n)
else:
n, res = self._get_bos_data_from_restart(p, n)
else:
# Pull rates, k from previous iteration w/o
# re-running transport
res = res_list[-1] # defined in previous i iteration
proc_time, n_list, res_list = self(n, res.rates, dt, p, i)
proc_time, n_end, res_end = self(n, res.rates, dt, p, i)
# Insert BOS concentration, transport results
n_list.insert(0, n)
res_list.insert(0, res)
# Remove actual EOS concentration for next step
n = n_list.pop()
StepResult.save(self.operator, n_list, res_list, [t, t + dt],
p, self._i_res + i, proc_time, path)
StepResult.save(
self.operator,
n,
res,
[t, t + dt],
p,
self._i_res + i,
proc_time,
write_rates=write_rates,
path=path
)
# Update for next step
n = n_end
res = res_end
t += dt
# No final simulation for SIE, use last iteration results
StepResult.save(self.operator, [n], [res_list[-1]], [t, t],
p, self._i_res + len(self), proc_time, path)
StepResult.save(
self.operator,
n,
res_end,
[t, t],
p,
self._i_res + len(self),
proc_time,
write_rates=write_rates,
path=path
)
self.operator.write_bos_data(self._i_res + len(self))
self.operator.finalize()

85
openmc/deplete/integrators.py
View file

@ -46,15 +46,12 @@ class PredictorIntegrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
n_end : list of numpy.ndarray
Concentrations at end of interval
op_results : empty list
Kept for consistency with API. No intermediate calls to operator
with predictor
"""
proc_time, n_end = self._timed_deplete(n, rates, dt, _i)
return proc_time, [n_end], []
return proc_time, n_end
@add_params
@ -98,11 +95,8 @@ class CECMIntegrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from transport simulations
n_end : list of numpy.ndarray
Concentrations at end of interval
"""
# deplete across first half of interval
time0, n_middle = self._timed_deplete(n, rates, dt / 2, _i)
@ -112,7 +106,7 @@ class CECMIntegrator(Integrator):
# MOS reaction rates
time1, n_end = self._timed_deplete(n, res_middle.rates, dt, _i)
return time0 + time1, [n_middle, n_end], [res_middle]
return time0 + time1, n_end
@add_params
@ -162,12 +156,8 @@ class CF4Integrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
n_end : list of numpy.ndarray
Concentrations at end of interval
"""
# Step 1: deplete with matrix 1/2*A(y0)
time1, n_eos1 = self._timed_deplete(
@ -192,9 +182,7 @@ class CF4Integrator(Integrator):
time5, n_eos5 = self._timed_deplete(
n_inter, list_rates, dt, _i, matrix_func=cf4_f4)
return (time1 + time2 + time3 + time4 + time5,
[n_eos1, n_eos2, n_eos3, n_eos5],
[res1, res2, res3])
return time1 + time2 + time3 + time4 + time5, n_eos5
@add_params
@ -240,12 +228,8 @@ class CELIIntegrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulation
n_end : list of numpy.ndarray
Concentrations at end of interval
"""
# deplete to end using BOS rates
proc_time, n_ce = self._timed_deplete(n_bos, rates, dt, _i)
@ -260,7 +244,7 @@ class CELIIntegrator(Integrator):
time_le2, n_end = self._timed_deplete(
n_inter, list_rates, dt, _i, matrix_func=celi_f2)
return proc_time + time_le1 + time_le1, [n_ce, n_end], [res_ce]
return proc_time + time_le1 + time_le2, n_end
@add_params
@ -306,12 +290,8 @@ class EPCRK4Integrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
n_end : list of numpy.ndarray
Concentrations at end of interval
"""
# Step 1: deplete with matrix A(y0) / 2
@ -330,7 +310,7 @@ class EPCRK4Integrator(Integrator):
list_rates = list(zip(rates, res1.rates, res2.rates, res3.rates))
time4, n4 = self._timed_deplete(n, list_rates, dt, _i, matrix_func=rk4_f4)
return (time1 + time2 + time3 + time4, [n1, n2, n3, n4], [res1, res2, res3])
return time1 + time2 + time3 + time4, n4
@add_params
@ -388,12 +368,8 @@ class LEQIIntegrator(Integrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulation
n_end : list of numpy.ndarray
Concentrations at end of interval
"""
if i == 0:
if self._i_res < 1: # need at least previous transport solution
@ -402,7 +378,7 @@ class LEQIIntegrator(Integrator):
self, n_bos, bos_rates, dt, source_rate, i)
prev_res = self.operator.prev_res[-2]
prev_dt = self.timesteps[i] - prev_res.time[0]
self._prev_rates = prev_res.rates[0]
self._prev_rates = prev_res.rates
else:
prev_dt = self.timesteps[i - 1]
@ -431,9 +407,7 @@ class LEQIIntegrator(Integrator):
# store updated rates
self._prev_rates = copy.deepcopy(bos_res.rates)
return (
time1 + time2 + time3 + time4, [n_eos0, n_eos1],
[bos_res, res_inter])
return time1 + time2 + time3 + time4, n_eos1
@add_params
@ -470,10 +444,9 @@ class SICELIIntegrator(SIIntegrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_bos_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
n_end : list of numpy.ndarray
Concentrations at end of interval
op_result : openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
"""
@ -499,7 +472,7 @@ class SICELIIntegrator(SIIntegrator):
proc_time += time1 + time2
# end iteration
return proc_time, [n_eos, n_inter], [res_bar]
return proc_time, n_inter, res_bar
@add_params
@ -536,10 +509,9 @@ class SILEQIIntegrator(SIIntegrator):
-------
proc_time : float
Time spent in CRAM routines for all materials in [s]
n_list : list of list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
n_end : list of numpy.ndarray
Concentrations at end of interval
op_result : openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulation
"""
@ -551,7 +523,7 @@ class SILEQIIntegrator(SIIntegrator):
self, n_bos, bos_rates, dt, source_rate, i)
prev_res = self.operator.prev_res[-2]
prev_dt = self.timesteps[i] - prev_res.time[0]
self._prev_rates = prev_res.rates[0]
self._prev_rates = prev_res.rates
else:
prev_dt = self.timesteps[i - 1]
@ -584,7 +556,10 @@ class SILEQIIntegrator(SIIntegrator):
n_inter, inputs, dt, i, matrix_func=leqi_f4)
proc_time += time1 + time2
return proc_time, [n_eos, n_inter], [res_bar]
# Store updated rates for next step
self._prev_rates = copy.deepcopy(bos_rates)
return proc_time, n_inter, res_bar
integrator_by_name = {

8
openmc/deplete/results.py
View file

@ -203,7 +203,7 @@ class Results(list):
# Evaluate value in each region
for i, result in enumerate(self):
times[i] = result.time[0]
concentrations[i] = result[0, mat_id, nuc]
concentrations[i] = result[mat_id, nuc]
# Unit conversions
times = _get_time_as(times, time_units)
@ -363,7 +363,7 @@ class Results(list):
# 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]
rates[i] = result.rates.get(mat_id, nuc, rx) * result[mat, nuc]
return times, rates
@ -397,7 +397,7 @@ class Results(list):
# Get time/eigenvalue at each point
for i, result in enumerate(self):
times[i] = result.time[0]
eigenvalues[i] = result.k[0]
eigenvalues[i] = result.k
# Convert time units if necessary
times = _get_time_as(times, time_units)
@ -630,7 +630,7 @@ class Results(list):
for nuc in result.index_nuc:
if nuc not in available_cross_sections:
continue
atoms = result[0, mat_id, nuc]
atoms = result[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

222
openmc/deplete/stepresult.py
View file

@ -16,7 +16,7 @@ from openmc.mpi import comm, MPI
from openmc.checkvalue import PathLike
from .reaction_rates import ReactionRates
VERSION_RESULTS = (1, 1)
VERSION_RESULTS = (1, 2)
__all__ = ["StepResult"]
@ -30,8 +30,8 @@ class StepResult:
Attributes
----------
k : list of (float, float)
Eigenvalue and uncertainty for each substep.
k : tuple of (float, float)
Eigenvalue and uncertainty at end of step.
time : list of float
Time at beginning, end of step, in seconds.
source_rate : float
@ -40,8 +40,8 @@ class StepResult:
Number of mats.
n_nuc : int
Number of nuclides.
rates : list of ReactionRates
The reaction rates for each substep.
rates : ReactionRates
The reaction rates at end of step.
volume : dict of str to float
Dictionary mapping mat id to volume.
index_mat : dict of str to int
@ -52,10 +52,8 @@ class StepResult:
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.
Atom quantity, stored by mat, then by nuclide.
proc_time : int
Average time spent depleting a material across all
materials and processes
@ -86,17 +84,17 @@ class StepResult:
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
A two-length tuple containing a mat index and a nuc
index. Both can be integers or slices, or can be
strings corresponding to their respective dictionary.
Returns
-------
float
The atoms for stage, mat, nuc
The atoms for mat, nuc
"""
stage, mat, nuc = pos
mat, nuc = pos
if isinstance(mat, openmc.Material):
mat = str(mat.id)
if isinstance(mat, str):
@ -104,7 +102,7 @@ class StepResult:
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
return self.data[stage, mat, nuc]
return self.data[mat, nuc]
def __setitem__(self, pos, val):
"""Sets an item from results.
@ -112,21 +110,21 @@ class StepResult:
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
A two-length tuple containing a mat index and a nuc
index. Both can be integers or slices, or can be
strings corresponding to their respective dictionary.
val : float
The value to set data to.
"""
stage, mat, nuc = pos
mat, nuc = pos
if isinstance(mat, str):
mat = self.index_mat[mat]
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
self.data[stage, mat, nuc] = val
self.data[mat, nuc] = val
@property
def n_mat(self):
@ -140,11 +138,7 @@ class StepResult:
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):
def allocate(self, volume, nuc_list, burn_list, full_burn_list):
"""Allocate memory for depletion step data
Parameters
@ -157,8 +151,6 @@ class StepResult:
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)
@ -167,7 +159,7 @@ class StepResult:
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))
self.data = np.zeros((self.n_mat, self.n_nuc))
def distribute(self, local_materials, ranges):
"""Create a new object containing data for distributed materials
@ -196,8 +188,8 @@ class StepResult:
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]
new.data = self.data[ranges]
new.rates = self.rates[ranges]
return new
def get_material(self, mat_id):
@ -232,7 +224,7 @@ class StepResult:
f'values are {list(self.volume.keys())}'
) from e
for nuc, _ in sorted(self.index_nuc.items(), key=lambda x: x[1]):
atoms = self[0, mat_id, nuc]
atoms = self[mat_id, nuc]
if atoms <= 0.0:
continue
atom_per_bcm = atoms / vol * 1e-24
@ -240,7 +232,7 @@ class StepResult:
material.volume = vol
return material
def export_to_hdf5(self, filename, step):
def export_to_hdf5(self, filename, step, write_rates: bool = False):
"""Export results to an HDF5 file
Parameters
@ -249,6 +241,8 @@ class StepResult:
The filename to write to
step : int
What step is this?
write_rates : bool, optional
Whether to include reaction rate datasets in the results file.
"""
# Write new file if first time step, else add to existing file
@ -259,7 +253,8 @@ class StepResult:
kwargs['driver'] = 'mpio'
kwargs['comm'] = comm
with h5py.File(filename, **kwargs) as handle:
self._to_hdf5(handle, step, parallel=True)
self._to_hdf5(handle, step, parallel=True,
write_rates=write_rates)
else:
# Gather results at root process
all_results = comm.gather(self)
@ -268,15 +263,18 @@ class StepResult:
if comm.rank == 0:
with h5py.File(filename, **kwargs) as handle:
for res in all_results:
res._to_hdf5(handle, step, parallel=False)
res._to_hdf5(handle, step, parallel=False,
write_rates=write_rates)
def _write_hdf5_metadata(self, handle):
def _write_hdf5_metadata(self, handle, write_rates):
"""Writes result metadata in HDF5 file
Parameters
----------
handle : h5py.File or h5py.Group
An hdf5 file or group type to store this in.
write_rates : bool
Whether reaction rate datasets are being written.
"""
# Create and save the 5 dictionaries:
@ -284,8 +282,8 @@ class StepResult:
# self.index_mat -> self.volume (TODO: support for changing volumes)
# self.index_nuc
# reactions
# self.rates[0].index_nuc (can be different from above, above is superset)
# self.rates[0].index_rx
# self.rates.index_nuc (can be different from above, above is superset)
# self.rates.index_rx
# these are shared by every step of the simulation, and should be deduplicated.
# Store concentration mat and nuclide dictionaries (along with volumes)
@ -295,13 +293,19 @@ class StepResult:
mat_list = sorted(self.mat_to_hdf5_ind, key=int)
nuc_list = sorted(self.index_nuc)
rxn_list = sorted(self.rates[0].index_rx)
include_rates = (
write_rates
and self.rates is not None
and bool(self.rates.index_nuc)
and bool(self.rates.index_rx)
)
rxn_list = sorted(self.rates.index_rx) if include_rates else []
n_mats = self.n_hdf5_mats
n_nuc_number = len(nuc_list)
n_nuc_rxn = len(self.rates[0].index_nuc)
n_nuc_rxn = len(self.rates.index_nuc) if include_rates else 0
n_rxn = len(rxn_list)
n_stages = self.n_stages
mat_group = handle.create_group("materials")
@ -315,41 +319,44 @@ class StepResult:
for nuc in nuc_list:
nuc_single_group = nuc_group.create_group(nuc)
nuc_single_group.attrs["atom number index"] = self.index_nuc[nuc]
if nuc in self.rates[0].index_nuc:
nuc_single_group.attrs["reaction rate index"] = self.rates[0].index_nuc[nuc]
if include_rates and nuc in self.rates.index_nuc:
nuc_single_group.attrs["reaction rate index"] = (
self.rates.index_nuc[nuc])
rxn_group = handle.create_group("reactions")
if include_rates:
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]
for rxn in rxn_list:
rxn_single_group = rxn_group.create_group(rxn)
rxn_single_group.attrs["index"] = (
self.rates.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),
handle.create_dataset("number", (1, n_mats, n_nuc_number),
maxshape=(None, n_mats, n_nuc_number),
chunks=True,
dtype='float64')
if n_nuc_rxn > 0 and n_rxn > 0:
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=True,
dtype='float64')
if include_rates and n_nuc_rxn > 0 and n_rxn > 0:
handle.create_dataset(
"reaction rates", (1, n_mats, n_nuc_rxn, n_rxn),
maxshape=(None, n_mats, n_nuc_rxn, n_rxn),
chunks=True, dtype='float64')
handle.create_dataset("eigenvalues", (1, n_stages, 2),
maxshape=(None, n_stages, 2), dtype='float64')
handle.create_dataset("eigenvalues", (1, 2),
maxshape=(None, 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),
handle.create_dataset("source_rate", (1,), maxshape=(None,),
dtype='float64')
handle.create_dataset(
"depletion time", (1,), maxshape=(None,),
dtype="float64")
def _to_hdf5(self, handle, index, parallel=False):
def _to_hdf5(self, handle, index, parallel=False, write_rates: bool = False):
"""Converts results object into an hdf5 object.
Parameters
@ -360,12 +367,14 @@ class StepResult:
What step is this?
parallel : bool
Being called with parallel HDF5?
write_rates : bool, optional
Whether reaction rate datasets are being written.
"""
if "/number" not in handle:
if parallel:
comm.barrier()
self._write_hdf5_metadata(handle)
self._write_hdf5_metadata(handle, write_rates)
if parallel:
comm.barrier()
@ -417,18 +426,14 @@ class StepResult:
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.index_mat]
low = min(inds)
high = max(inds)
for i in range(n_stages):
number_dset[index, i, low:high+1] = self.data[i]
if has_reactions:
rxn_dset[index, i, low:high+1] = self.rates[i]
if comm.rank == 0:
eigenvalues_dset[index, i] = self.k[i]
number_dset[index, low:high+1] = self.data
if has_reactions:
rxn_dset[index, low:high+1] = self.rates
if comm.rank == 0:
eigenvalues_dset[index] = self.k
time_dset[index] = self.time
source_rate_dset[index] = self.source_rate
if self.proc_time is not None:
@ -459,10 +464,24 @@ class StepResult:
# Older versions used "power" instead of "source_rate"
source_rate_dset = handle["/power"]
results.data = number_dset[step, :, :, :]
results.k = eigenvalues_dset[step, :]
# Check if this is an old format file (with stages dimension) or new format
# Old format: number has shape (n_steps, n_stages, n_mats, n_nucs)
# New format: number has shape (n_steps, n_mats, n_nucs)
has_stages = len(number_dset.shape) == 4
if has_stages:
# Old format - extract data from first stage (index 0)
results.data = number_dset[step, 0, :, :]
results.k = eigenvalues_dset[step, 0, :]
# source_rate had shape (n_steps, n_stages) in old format
results.source_rate = source_rate_dset[step, 0]
else:
# New format - no stages dimension
results.data = number_dset[step, :, :]
results.k = eigenvalues_dset[step, :]
results.source_rate = source_rate_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"]
@ -493,33 +512,45 @@ class StepResult:
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"]
if "reactions" in handle:
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.index_mat, rxn_nuc_to_ind, rxn_to_ind, True)
if "reaction rates" in handle:
rate[:] = handle["/reaction rates"][step, i, :, :, :]
results.rates.append(rate)
# Reconstruct reaction rates
rate = ReactionRates(results.index_mat, rxn_nuc_to_ind, rxn_to_ind, True)
if "reaction rates" in handle:
if has_stages:
# Old format: (n_steps, n_stages, n_mats, n_nucs, n_rxns)
rate[:] = handle["/reaction rates"][step, 0, :, :, :]
else:
# New format: (n_steps, n_mats, n_nucs, n_rxns)
rate[:] = handle["/reaction rates"][step, :, :, :]
results.rates = rate
return results
@staticmethod
def save(op, x, op_results, t, source_rate, step_ind, proc_time=None,
path: PathLike = "depletion_results.h5"):
def save(
op,
x,
op_results,
t,
source_rate,
step_ind,
proc_time=None,
write_rates: bool = False,
path: PathLike = "depletion_results.h5"
):
"""Creates and writes depletion results to disk
Parameters
----------
op : openmc.deplete.abc.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
x : numpy.array
End-of-step concentrations for each material
op_results : openmc.deplete.OperatorResult
Result of applying transport operator at end of step
t : list of float
Time indices.
source_rate : float
@ -530,7 +561,8 @@ class StepResult:
Total process time spent depleting materials. This may
be process-dependent and will be reduced across MPI
processes.
write_rates : bool, optional
Whether reaction rates should be written to the results file.
path : PathLike
Path to file to write. Defaults to 'depletion_results.h5'.
@ -539,26 +571,20 @@ class StepResult:
# 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)
results.allocate(vol_dict, nuc_list, burn_list, full_burn_list)
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]
for mat_i in range(n_mat):
results[mat_i, :] = x[mat_i]
ks = []
for r in op_results:
if isinstance(r.k, type(None)):
ks += [(None, None)]
else:
ks += [(r.k.nominal_value, r.k.std_dev)]
results.k = ks
results.rates = [r.rates for r in op_results]
if isinstance(op_results.k, type(None)):
results.k = (None, None)
else:
results.k = (op_results.k.nominal_value, op_results.k.std_dev)
results.rates = op_results.rates
results.time = t
results.source_rate = source_rate
results.proc_time = proc_time
@ -567,7 +593,7 @@ class StepResult:
if not Path(path).is_file():
Path(path).parent.mkdir(parents=True, exist_ok=True)
results.export_to_hdf5(path, step_ind)
results.export_to_hdf5(path, step_ind, write_rates)
def transfer_volumes(self, model):
"""Transfers volumes from depletion results to geometry

8
tests/unit_tests/test_deplete_continue.py
View file

@ -17,7 +17,7 @@ def test_continue(run_in_tmpdir):
operator = dummy_operator.DummyOperator()
# initial depletion
bundle.solver(operator, [1.0, 2.0], [1.0, 2.0]).integrate()
bundle.solver(operator, [1.0, 2.0], [1.0, 2.0]).integrate(write_rates=True)
# set up continue run
prev_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
@ -25,7 +25,7 @@ def test_continue(run_in_tmpdir):
# if continue run happens, test passes
bundle.solver(operator, [1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 4.0],
continue_timesteps=True).integrate()
continue_timesteps=True).integrate(write_rates=True)
final_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
@ -42,7 +42,7 @@ def test_continue_continue(run_in_tmpdir):
operator = dummy_operator.DummyOperator()
# initial depletion
bundle.solver(operator, [1.0, 2.0], [1.0, 2.0]).integrate()
bundle.solver(operator, [1.0, 2.0], [1.0, 2.0]).integrate(write_rates=True)
# set up continue run
prev_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
@ -50,7 +50,7 @@ def test_continue_continue(run_in_tmpdir):
# first continue run
bundle.solver(operator, [1.0, 2.0, 3.0, 4.0], [1.0, 2.0, 3.0, 4.0],
continue_timesteps=True).integrate()
continue_timesteps=True).integrate(write_rates=True)
prev_res = openmc.deplete.Results(operator.output_dir / "depletion_results.h5")
# second continue run

93
tests/unit_tests/test_deplete_integrator.py
View file

@ -10,6 +10,7 @@ import copy
from random import uniform
from unittest.mock import MagicMock
import h5py
import numpy as np
from uncertainties import ufloat
import pytest
@ -38,8 +39,6 @@ INTEGRATORS = [
def test_results_save(run_in_tmpdir):
"""Test data save module"""
stages = 3
rng = np.random.RandomState(comm.rank)
# Mock geometry
@ -63,31 +62,22 @@ def test_results_save(run_in_tmpdir):
op.get_results_info.return_value = (
vol_dict, nuc_list, burn_list, full_burn_list)
# Construct x
x1 = []
x2 = []
# Construct end-of-step concentrations
x1 = [rng.random(2), rng.random(2)]
x2 = [rng.random(2), rng.random(2)]
for i in range(stages):
x1.append([rng.random(2), rng.random(2)])
x2.append([rng.random(2), rng.random(2)])
# Construct r
# Construct reaction rates
r1 = ReactionRates(burn_list, ["na", "nb"], ["ra", "rb"])
r1[:] = rng.random((2, 2, 2))
rate1 = copy.deepcopy(r1)
rate1 = []
rate2 = []
r2 = ReactionRates(burn_list, ["na", "nb"], ["ra", "rb"])
r2[:] = rng.random((2, 2, 2))
rate2 = copy.deepcopy(r2)
for i in range(stages):
rate1.append(copy.deepcopy(r1))
r1[:] = rng.random((2, 2, 2))
rate2.append(copy.deepcopy(r1))
r1[:] = rng.random((2, 2, 2))
# Create global terms
# Col 0: eig, Col 1: uncertainty
eigvl1 = rng.random((stages, 2))
eigvl2 = rng.random((stages, 2))
# Create global terms (eigenvalue and uncertainty)
eigvl1 = rng.random(2)
eigvl2 = rng.random(2)
eigvl1 = comm.bcast(eigvl1, root=0)
eigvl2 = comm.bcast(eigvl2, root=0)
@ -95,29 +85,35 @@ def test_results_save(run_in_tmpdir):
t1 = [0.0, 1.0]
t2 = [1.0, 2.0]
op_result1 = [OperatorResult(ufloat(*k), rates)
for k, rates in zip(eigvl1, rate1)]
op_result2 = [OperatorResult(ufloat(*k), rates)
for k, rates in zip(eigvl2, rate2)]
op_result1 = OperatorResult(ufloat(*eigvl1), rate1)
op_result2 = OperatorResult(ufloat(*eigvl2), rate2)
# saves within a subdirectory
StepResult.save(op, x1, op_result1, t1, 0, 0, path='out/put/depletion.h5')
StepResult.save(
op,
x1,
op_result1,
t1,
0,
0,
write_rates=True,
path='out/put/depletion.h5'
)
res = Results('out/put/depletion.h5')
# saves with default filename
StepResult.save(op, x1, op_result1, t1, 0, 0)
StepResult.save(op, x2, op_result2, t2, 0, 1)
StepResult.save(op, x1, op_result1, t1, 0, 0, write_rates=True)
StepResult.save(op, x2, op_result2, t2, 0, 1, write_rates=True)
# Load the files
res = Results("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])
for mat_i, mat in enumerate(burn_list):
for nuc_i, nuc in enumerate(nuc_list):
assert res[0][mat, nuc] == x1[mat_i][nuc_i]
assert res[1][mat, nuc] == x2[mat_i][nuc_i]
np.testing.assert_array_equal(res[0].rates, rate1)
np.testing.assert_array_equal(res[1].rates, rate2)
np.testing.assert_array_equal(res[0].k, eigvl1)
np.testing.assert_array_equal(res[0].time, t1)
@ -126,6 +122,31 @@ def test_results_save(run_in_tmpdir):
np.testing.assert_array_equal(res[1].time, t2)
def test_results_save_without_rates(run_in_tmpdir):
"""StepResult.save skips reaction-rate datasets by default"""
op = MagicMock()
op.prev_res = None
vol_dict = {"0": 1.0}
nuc_list = ["na"]
burn_list = ["0"]
op.get_results_info.return_value = (vol_dict, nuc_list, burn_list, burn_list)
x = [np.array([1.0])]
rates = ReactionRates(burn_list, nuc_list, ["ra"])
rates[:] = np.array([[[2.0]]])
op_result = OperatorResult(ufloat(1.0, 0.1), rates)
StepResult.save(op, x, op_result, [0.0, 1.0], 0.0, 0)
with h5py.File('depletion_results.h5', 'r') as handle:
assert 'reaction rates' not in handle
assert 'reactions' not in handle
res = Results('depletion_results.h5')
assert res[0].rates.size == 0
def test_bad_integrator_inputs():
"""Test failure modes for Integrator inputs"""

14
tests/unit_tests/test_deplete_restart.py
View file

@ -21,7 +21,7 @@ def test_restart_predictor_cecm(run_in_tmpdir):
# Perform simulation using the predictor algorithm
dt = [0.75]
power = 1.0
openmc.deplete.PredictorIntegrator(op, dt, power).integrate()
openmc.deplete.PredictorIntegrator(op, dt, power).integrate(write_rates=True)
# Load the files
prev_res = openmc.deplete.Results(op.output_dir / "depletion_results.h5")
@ -30,10 +30,6 @@ def test_restart_predictor_cecm(run_in_tmpdir):
op = dummy_operator.DummyOperator(prev_res)
op.output_dir = output_dir
# check ValueError is raised, indicating previous and current stages
with pytest.raises(ValueError, match="incompatible.* 1.*2"):
openmc.deplete.CECMIntegrator(op, dt, power)
def test_restart_cecm_predictor(run_in_tmpdir):
"""Integral regression test of integrator algorithm using CE/CM for the
@ -47,7 +43,7 @@ def test_restart_cecm_predictor(run_in_tmpdir):
dt = [0.75]
power = 1.0
cecm = openmc.deplete.CECMIntegrator(op, dt, power)
cecm.integrate()
cecm.integrate(write_rates=True)
# Load the files
prev_res = openmc.deplete.Results(op.output_dir / "depletion_results.h5")
@ -56,10 +52,6 @@ def test_restart_cecm_predictor(run_in_tmpdir):
op = dummy_operator.DummyOperator(prev_res)
op.output_dir = output_dir
# check ValueError is raised, indicating previous and current stages
with pytest.raises(ValueError, match="incompatible.* 2.*1"):
openmc.deplete.PredictorIntegrator(op, dt, power)
@pytest.mark.parametrize("scheme", dummy_operator.SCHEMES)
def test_restart(run_in_tmpdir, scheme):
@ -70,7 +62,7 @@ def test_restart(run_in_tmpdir, scheme):
operator = dummy_operator.DummyOperator()
# take first step
bundle.solver(operator, [0.75], 1.0).integrate()
bundle.solver(operator, [0.75], 1.0).integrate(write_rates=True)
# restart
prev_res = openmc.deplete.Results(