Merge pull request #2037 from paulromano/depletion-docstring-fixes

Collection of depletion-related fixes and enhancements
This commit is contained in:
Paul Romano 2022-04-24 12:24:01 -05:00 committed by GitHub
commit b5b928e2e0
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
6 changed files with 88 additions and 39 deletions

View file

@ -851,7 +851,7 @@ class Integrator(ABC):
return (self.operator.prev_res[-1].time[-1],
len(self.operator.prev_res) - 1)
def integrate(self, final_step=True):
def integrate(self, final_step=True, output=True):
"""Perform the entire depletion process across all steps
Parameters
@ -861,12 +861,18 @@ class Integrator(ABC):
of the last timestep.
.. versionadded:: 0.12.1
output : bool, optional
Indicate whether to display information about progress
.. versionadded:: 0.13.1
"""
with self.operator as conc:
t, self._i_res = self._get_start_data()
for i, (dt, source_rate) in enumerate(self):
if output:
print(f"[openmc.deplete] t={t} s, dt={dt} s, source={source_rate}")
# Solve transport equation (or obtain result from restart)
if i > 0 or self.operator.prev_res is None:
conc, res = self._get_bos_data_from_operator(i, source_rate, conc)
@ -892,6 +898,8 @@ class Integrator(ABC):
# source rate is passed to the transport operator (which knows to
# just return zero reaction rates without actually doing a transport
# solve)
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],
source_rate, self._i_res + len(self), proc_time)
@ -1004,12 +1012,23 @@ class SIIntegrator(Integrator):
self.operator.settings.particles //= self.n_steps
return inherited
def integrate(self):
"""Perform the entire depletion process across all steps"""
def integrate(self, output=True):
"""Perform the entire depletion process across all steps
Parameters
----------
output : bool, optional
Indicate whether to display information about progress
.. versionadded:: 0.13.1
"""
with self.operator as conc:
t, self._i_res = self._get_start_data()
for i, (dt, p) in enumerate(self):
if output:
print(f"[openmc.deplete] t={t} s, dt={dt} s, source={p}")
if i == 0:
if self.operator.prev_res is None:
conc, res = self._get_bos_data_from_operator(i, p, conc)

View file

@ -6,6 +6,8 @@ from collections import OrderedDict
import numpy as np
from openmc import Material
class AtomNumber:
"""Stores local material compositions (atoms of each nuclide).
@ -58,6 +60,12 @@ class AtomNumber:
self.number = np.zeros((len(local_mats), len(nuclides)))
def _get_mat_index(self, mat):
"""Helper method for getting material index"""
if isinstance(mat, Material):
mat = str(mat.id)
return self.index_mat[mat] if isinstance(mat, str) else mat
def __getitem__(self, pos):
"""Retrieves total atom number from AtomNumber.
@ -75,8 +83,7 @@ class AtomNumber:
"""
mat, nuc = pos
if isinstance(mat, str):
mat = self.index_mat[mat]
mat = self._get_mat_index(mat)
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
@ -96,8 +103,7 @@ class AtomNumber:
"""
mat, nuc = pos
if isinstance(mat, str):
mat = self.index_mat[mat]
mat = self._get_mat_index(mat)
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
@ -121,11 +127,11 @@ class AtomNumber:
if ind < self.n_nuc_burn]
def get_atom_density(self, mat, nuc):
"""Accesses atom density instead of total number.
"""Return atom density of given material and nuclide
Parameters
----------
mat : str, int or slice
mat : str, int, openmc.Material or slice
Material index.
nuc : str, int or slice
Nuclide index.
@ -136,19 +142,43 @@ class AtomNumber:
Density in [atom/cm^3]
"""
if isinstance(mat, str):
mat = self.index_mat[mat]
mat = self._get_mat_index(mat)
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
return self[mat, nuc] / self.volume[mat]
def get_atom_densities(self, mat, units='atom/b-cm'):
"""Return atom densities for a given material
Parameters
----------
mat : str, int, openmc.Material or slice
Material index.
units : {"atom/b-cm", "atom/cm3"}, optional
Units for the returned concentration. Default is ``"atom/b-cm"``
.. versionadded:: 0.13.1
Returns
-------
dict
Dictionary mapping nuclides to atom densities
"""
mat = self._get_mat_index(mat)
normalization = (1.0e-24 if units == 'atom/b-cm' else 1.0) / self.volume[mat]
return {
name: normalization * self[mat, nuc]
for name, nuc in self.index_nuc.items()
}
def set_atom_density(self, mat, nuc, val):
"""Sets atom density instead of total number.
Parameters
----------
mat : str, int or slice
mat : str, int, openmc.Material or slice
Material index.
nuc : str, int or slice
Nuclide index.
@ -156,8 +186,7 @@ class AtomNumber:
Array of densities to set in [atom/cm^3]
"""
if isinstance(mat, str):
mat = self.index_mat[mat]
mat = self._get_mat_index(mat)
if isinstance(nuc, str):
nuc = self.index_nuc[nuc]
@ -168,7 +197,7 @@ class AtomNumber:
Parameters
----------
mat : str, int or slice
mat : str, int, openmc.Material or slice
Material index.
Returns
@ -177,9 +206,7 @@ class AtomNumber:
The slice requested in [atom].
"""
if isinstance(mat, str):
mat = self.index_mat[mat]
mat = self._get_mat_index(mat)
return self[mat, :self.n_nuc_burn]
def set_mat_slice(self, mat, val):
@ -187,15 +214,13 @@ class AtomNumber:
Parameters
----------
mat : str, int or slice
mat : str, int, openmc.Material, or slice
Material index.
val : numpy.ndarray
The slice to set in [atom]
"""
if isinstance(mat, str):
mat = self.index_mat[mat]
mat = self._get_mat_index(mat)
self[mat, :self.n_nuc_burn] = val
def set_density(self, total_density):

View file

@ -124,13 +124,13 @@ class CF4Integrator(Integrator):
.. math::
\begin{aligned}
\mathbf{A}_1 &= h\mathbf{A}(\mathbf{n}_0) \\
\hat{\mathbf{n}}_1 &= \exp \left ( \frac{\mathbf{A}_1}{2} \right ) \\
\mathbf{A}_1 &= h\mathbf{A}(\mathbf{n}_i) \\
\hat{\mathbf{n}}_1 &= \exp \left ( \frac{\mathbf{A}_1}{2} \right ) \mathbf{n}_i \\
\mathbf{A}_2 &= h\mathbf{A}(\hat{\mathbf{n}}_1) \\
\hat{\mathbf{n}}_2 &= \exp \left ( \frac{\mathbf{A}_2}{2} \right ) \\
\hat{\mathbf{n}}_2 &= \exp \left ( \frac{\mathbf{A}_2}{2} \right ) \mathbf{n}_i \\
\mathbf{A}_3 &= h \mathbf{A}(\hat{\mathbf{n}}_2) \\
\hat{\mathbf{n}}_3 &= \exp \left ( -\frac{\mathbf{A}_1}{2} + \mathbf{A}_3
\right ) \\
\right ) \hat{\mathbf{n}}_1 \\
\mathbf{A}_4 &= h\mathbf{A}(\hat{\mathbf{n}}_3) \\
\mathbf{n}_{i+1} &= \exp \left ( \frac{\mathbf{A}_1}{4} + \frac{\mathbf{A}_2}{6}
+ \frac{\mathbf{A}_3}{6} - \frac{\mathbf{A}_4}{12} \right )
@ -207,7 +207,8 @@ class CELIIntegrator(Integrator):
.. math::
\begin{aligned}
\mathbf{n}_{i+1}^p &= \exp \left ( h \mathbf{A}(\mathbf{n}_i ) \right ) \\
\mathbf{n}_{i+1}^p &= \exp \left ( h \mathbf{A}(\mathbf{n}_i ) \right )
\mathbf{n}_i \\
\mathbf{n}_{i+1} &= \exp \left( \frac{h}{12} \mathbf{A}(\mathbf{n}_i) +
\frac{5h}{12} \mathbf{A}(\mathbf{n}_{i+1}^p) \right)
\exp \left( \frac{5h}{12} \mathbf{A}(\mathbf{n}_i) +
@ -268,12 +269,12 @@ class EPCRK4Integrator(Integrator):
.. math::
\begin{aligned}
\mathbf{A}_1 &= h\mathbf{A}(\mathbf{n}_0) \\
\hat{\mathbf{n}}_1 &= \exp \left ( \frac{\mathbf{A}_1}{2} \right ) \\
\mathbf{A}_1 &= h\mathbf{A}(\mathbf{n}_i) \\
\hat{\mathbf{n}}_1 &= \exp \left ( \frac{\mathbf{A}_1}{2} \right ) \mathbf{n}_i \\
\mathbf{A}_2 &= h\mathbf{A}(\hat{\mathbf{n}}_1) \\
\hat{\mathbf{n}}_2 &= \exp \left ( \frac{\mathbf{A}_2}{2} \right ) \\
\hat{\mathbf{n}}_2 &= \exp \left ( \frac{\mathbf{A}_2}{2} \right ) \mathbf{n}_i \\
\mathbf{A}_3 &= h \mathbf{A}(\hat{\mathbf{n}}_2) \\
\hat{\mathbf{n}}_3 &= \exp \left ( \mathbf{A}_3 \right ) \\
\hat{\mathbf{n}}_3 &= \exp \left ( \mathbf{A}_3 \right ) \mathbf{n}_i \\
\mathbf{A}_4 &= h\mathbf{A}(\hat{\mathbf{n}}_3) \\
\mathbf{n}_{i+1} &= \exp \left ( \frac{\mathbf{A}_1}{6} + \frac{\mathbf{A}_2}{3}
+ \frac{\mathbf{A}_3}{3} + \frac{\mathbf{A}_4}{6} \right ) \mathbf{n}_i.

View file

@ -123,6 +123,10 @@ class Nuclide:
# Neutron fission yields, if present
self._yield_data = None
def __repr__(self):
n_modes, n_rx = self.n_decay_modes, self.n_reaction_paths
return f"<Nuclide: {self.name} ({n_modes} modes, {n_rx} reactions)>"
@property
def n_decay_modes(self):
return len(self.decay_modes)

View file

@ -84,11 +84,11 @@ class Operator(TransportOperator):
.. versionchanged:: 0.13.0
The geometry and settings parameters have been replaced with a
model parameter that takes an openmc.Model object
model parameter that takes a :class:`~openmc.model.Model` object
Parameters
----------
model : openmc.Model
model : openmc.model.Model
OpenMC model object
chain_file : str, optional
Path to the depletion chain XML file. Defaults to the file
@ -163,7 +163,7 @@ class Operator(TransportOperator):
Attributes
----------
model : openmc.Model
model : openmc.model.Model
OpenMC model object
geometry : openmc.Geometry
OpenMC geometry object
@ -195,8 +195,6 @@ class Operator(TransportOperator):
prev_res : ResultsList or None
Results from a previous depletion calculation. ``None`` if no
results are to be used.
diff_burnable_mats : bool
Whether to differentiate burnable materials with multiple instances
cleanup_when_done : bool
Whether to finalize and clear the shared library memory when the
depletion operation is complete. Defaults to clearing the library.
@ -248,7 +246,6 @@ class Operator(TransportOperator):
)
self.materials = model.materials
self.diff_burnable_mats = diff_burnable_mats
self.cleanup_when_done = True
# Reduce the chain before we create more materials
@ -262,8 +259,11 @@ class Operator(TransportOperator):
self.chain = self.chain.reduce(all_isotopes, reduce_chain_level)
# Differentiate burnable materials with multiple instances
if self.diff_burnable_mats:
if diff_burnable_mats:
self._differentiate_burnable_mats()
self.materials = openmc.Materials(
model.geometry.get_all_materials().values()
)
# Clear out OpenMC, create task lists, distribute
openmc.reset_auto_ids()

View file

@ -545,7 +545,7 @@ void print_results()
fmt::print(" Leakage Fraction = {:.5f}\n",
gt(GlobalTally::LEAKAGE, TallyResult::SUM) / n);
}
fmt::print("\n");
std::cout << std::endl;
}
//==============================================================================