Change StatePoint.k_combined -> keff, ResultsList.get_eigenvalue -> get_keff

This commit is contained in:
Paul Romano 2022-04-28 15:43:42 -05:00
parent 99b128c421
commit 45a6b31c8f
20 changed files with 67 additions and 48 deletions

View file

@ -54,7 +54,7 @@ easy retrieval of k-effective, nuclide concentrations, and reaction rates over
time::
results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
time, keff = results.get_eigenvalue()
time, keff = results.get_keff()
Note that the coupling between the transport solver and the transmutation solver
happens in-memory rather than by reading/writing files on disk.

View file

@ -81,7 +81,7 @@ it is recommended that data be extracted from statepoint files in a context mana
.. code-block:: python
with openmc.StatePoint('statepoint.10.h5') as sp:
k_eff = sp.k_combined
k_eff = sp.keff
or that the :meth:`StatePoint.close` method is called before executing a subsequent
OpenMC run.

View file

@ -29,5 +29,5 @@ openmc.run()
# Get the resulting k-effective value
n = settings.batches
with openmc.StatePoint(f'statepoint.{n}.h5') as sp:
keff = sp.k_combined
keff = sp.keff
print(f'Final k-effective = {keff}')

View file

@ -59,33 +59,33 @@ integrator.integrate()
results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
# Obtain K_eff as a function of time
time, keff = results.get_eigenvalue()
time, keff = results.get_keff(time_units='d')
# Obtain U235 concentration as a function of time
time, n_U235 = results.get_atoms(uo2, 'U235')
uo2 = geometry.get_all_material_cells()[1]
_, n_U235 = results.get_atoms(uo2, 'U235')
# Obtain Xe135 capture reaction rate as a function of time
time, Xe_capture = results.get_reaction_rate(uo2, 'Xe135', '(n,gamma)')
_, Xe_capture = results.get_reaction_rate(uo2, 'Xe135', '(n,gamma)')
###############################################################################
# Generate plots
###############################################################################
days = 24*60*60
fig, ax = plt.subplots()
ax.errorbar(time/days, keff[:, 0], keff[:, 1], label="K-effective")
ax.errorbar(time, keff[:, 0], keff[:, 1], label="K-effective")
ax.set_xlabel("Time [d]")
ax.set_ylabel("Keff")
plt.show()
fig, ax = plt.subplots()
ax.plot(time/days, n_U235, label="U235")
ax.plot(time, n_U235, label="U235")
ax.set_xlabel("Time [d]")
ax.set_ylabel("U235 atoms")
plt.show()
fig, ax = plt.subplots()
ax.plot(time/days, Xe_capture, label="Xe135 capture")
ax.plot(time, Xe_capture, label="Xe135 capture")
ax.set_xlabel("Time [d]")
ax.set_ylabel("Xe135 capture rate")
plt.show()

View file

@ -104,33 +104,32 @@ integrator.integrate()
results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")
# Obtain K_eff as a function of time
time, keff = results.get_eigenvalue()
time, keff = results.get_keff(time_units='d')
# Obtain U235 concentration as a function of time
time, n_U235 = results.get_atoms(uo2, 'U235')
_, n_U235 = results.get_atoms(uo2, 'U235')
# Obtain Xe135 capture reaction rate as a function of time
time, Xe_capture = results.get_reaction_rate(uo2, 'Xe135', '(n,gamma)')
_, Xe_capture = results.get_reaction_rate(uo2, 'Xe135', '(n,gamma)')
###############################################################################
# Generate plots
###############################################################################
days = 24*60*60
fig, ax = plt.subplots()
ax.errorbar(time/days, keff[:, 0], keff[:, 1], label="K-effective")
ax.errorbar(time, keff[:, 0], keff[:, 1], label="K-effective")
ax.set_xlabel("Time [d]")
ax.set_ylabel("Keff")
plt.show()
fig, ax = plt.subplots()
ax.plot(time/days, n_U235, label="U235")
ax.plot(time, n_U235, label="U235")
ax.set_xlabel("Time [d]")
ax.set_ylabel("U235 atoms")
plt.show()
fig, ax = plt.subplots()
ax.plot(time/days, Xe_capture, label="Xe135 capture")
ax.plot(time, Xe_capture, label="Xe135 capture")
ax.set_xlabel("Time [d]")
ax.set_ylabel("Xe135 capture rate")
plt.show()

View file

@ -725,7 +725,7 @@ class Operator(TransportOperator):
rates.fill(0.0)
# Get k and uncertainty
k_combined = ufloat(*openmc.lib.keff())
keff = ufloat(*openmc.lib.keff())
# Extract tally bins
nuclides = self._rate_helper.nuclides
@ -779,7 +779,7 @@ class Operator(TransportOperator):
# Store new fission yields on the chain
self.chain.fission_yields = fission_yields
return OperatorResult(k_combined, rates)
return OperatorResult(keff, rates)
def _get_nuclides_with_data(self, cross_sections):
"""Loads cross_sections.xml file to find nuclides with neutron data"""

View file

@ -50,7 +50,7 @@ class Results:
Number of stages in simulation.
data : numpy.ndarray
Atom quantity, stored by stage, mat, then by nuclide.
proc_time: int
proc_time : int
Average time spent depleting a material across all
materials and processes
@ -518,4 +518,4 @@ class Results:
for material in materials:
if material.depletable:
material.volume = self.volume[str(material.id)]
material.volume = self.volume[str(material.id)]

View file

@ -1,6 +1,7 @@
import numbers
import bisect
import math
from warnings import warn
import h5py
import numpy as np
@ -165,7 +166,7 @@ class ResultsList(list):
return times, rates
def get_eigenvalue(self, time_units='s'):
def get_keff(self, time_units='s'):
"""Evaluates the eigenvalue from a results list.
Parameters
@ -197,14 +198,19 @@ class ResultsList(list):
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::
Will have one fewer row than number of other methods,
like :meth:`get_eigenvalues`, because no depletion
is performed at the final transport stage
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
-------

View file

@ -593,7 +593,7 @@ class Library:
self._nuclides = statepoint.summary.nuclides
if statepoint.run_mode == 'eigenvalue':
self._keff = statepoint.k_combined.n
self._keff = statepoint.keff.n
# Load tallies for each MGXS for each domain and mgxs type
for domain in self.domains:

View file

@ -53,7 +53,7 @@ def _search_keff(guess, target, model_builder, model_args, print_iterations,
# Run the model and obtain keff
sp_filepath = model.run(output=print_output)
with openmc.StatePoint(sp_filepath) as sp:
keff = sp.k_combined
keff = sp.keff
# Record the history
guesses.append(guess)

View file

@ -63,6 +63,8 @@ class StatePoint:
datatype has fields 'name', 'sum', 'sum_sq', 'mean', and 'std_dev'.
k_combined : uncertainties.UFloat
Combined estimator for k-effective
.. deprecated:: 0.13.1
k_col_abs : float
Cross-product of collision and absorption estimates of k-effective
k_col_tra : float
@ -71,6 +73,10 @@ class StatePoint:
Cross-product of absorption and tracklength estimates of k-effective
k_generation : numpy.ndarray
Estimate of k-effective for each batch/generation
keff : uncertainties.UFloat
Combined estimator for k-effective
.. versionadded:: 0.13.1
meshes : dict
Dictionary whose keys are mesh IDs and whose values are MeshBase objects
n_batches : int
@ -260,12 +266,20 @@ class StatePoint:
return None
@property
def k_combined(self):
def keff(self):
if self.run_mode == 'eigenvalue':
return ufloat(*self._f['k_combined'][()])
else:
return None
@property
def k_combined(self):
warnings.warn(
"The 'k_combined' property has been renamed to 'keff' and will be "
"removed in a future version of OpenMC.", FutureWarning
)
return self.keff
@property
def k_col_abs(self):
if self.run_mode == 'eigenvalue':

View file

@ -115,16 +115,16 @@ def test_full(run_in_tmpdir, problem, multiproc):
# Compare statepoint files with depletion results
t_test, k_test = res_test.get_eigenvalue()
t_ref, k_ref = res_ref.get_eigenvalue()
t_test, k_test = res_test.get_keff()
t_ref, k_ref = res_ref.get_keff()
k_state = np.empty_like(k_ref)
n_tallies = np.empty(N + 1, dtype=int)
# Get statepoint files for all BOS points and EOL
for n in range(N + 1):
statepoint = openmc.StatePoint("openmc_simulation_n{}.h5".format(n))
k_n = statepoint.k_combined
statepoint = openmc.StatePoint(f"openmc_simulation_n{n}.h5")
k_n = statepoint.keff
k_state[n] = [k_n.nominal_value, k_n.std_dev]
n_tallies[n] = len(statepoint.tallies)
# Look for exact match pulling from statepoint and depletion_results

View file

@ -14,7 +14,7 @@ class EntropyTestHarness(TestHarness):
with StatePoint(statepoint) as sp:
# Write out k-combined.
outstr = 'k-combined:\n'
outstr += '{:12.6E} {:12.6E}\n'.format(sp.k_combined.n, sp.k_combined.s)
outstr += '{:12.6E} {:12.6E}\n'.format(sp.keff.n, sp.keff.s)
# Write out entropy data.
outstr += 'entropy:\n'

View file

@ -145,7 +145,7 @@ class MGXSTestHarness(PyAPITestHarness):
# Write out k-combined.
outstr += 'k-combined:\n'
form = '{:12.6E} {:12.6E}\n'
outstr += form.format(sp.k_combined.n, sp.k_combined.s)
outstr += form.format(sp.keff.n, sp.keff.s)
return outstr

View file

@ -90,7 +90,7 @@ class TriggerStatepointRestartTestHarness(PyAPITestHarness):
assert spfile
with openmc.StatePoint(spfile) as sp:
sp_batchno_1 = sp.current_batch
k_combined_1 = sp.k_combined
keff_1 = sp.keff
assert sp_batchno_1 > 5
print('Last batch no = %d' % sp_batchno_1)
self._write_inputs(self._get_inputs())
@ -108,13 +108,13 @@ class TriggerStatepointRestartTestHarness(PyAPITestHarness):
assert spfile
with openmc.StatePoint(spfile) as sp:
sp_batchno_2 = sp.current_batch
k_combined_2 = sp.k_combined
keff_2 = sp.keff
assert sp_batchno_2 > 5
assert sp_batchno_1 == sp_batchno_2, \
'Different final batch number after restart'
# need str() here as uncertainties.ufloat instances are always different
assert str(k_combined_1) == str(k_combined_2), \
'Different final k_combined after restart'
assert str(keff_1) == str(keff_2), \
'Different final keff after restart'
self._write_inputs(self._get_inputs())
self._compare_inputs()
self._test_output_created()

View file

@ -92,7 +92,7 @@ class TestHarness:
# Write out k-combined.
outstr += 'k-combined:\n'
form = '{0:12.6E} {1:12.6E}\n'
outstr += form.format(sp.k_combined.n, sp.k_combined.s)
outstr += form.format(sp.keff.n, sp.keff.s)
# Write out tally data.
for i, tally_ind in enumerate(sp.tallies):

View file

@ -55,10 +55,10 @@ def test_get_reaction_rate(res):
np.testing.assert_allclose(r, np.array(n_ref) * xs_ref)
def test_get_eigenvalue(res):
"""Tests evaluating eigenvalue."""
t, k = res.get_eigenvalue()
t_min, k = res.get_eigenvalue(time_units='min')
def test_get_keff(res):
"""Tests evaluating keff."""
t, k = res.get_keff()
t_min, k = res.get_keff(time_units='min')
t_ref = [0.0, 1296000.0, 2592000.0, 3888000.0]
k_ref = [1.21409662, 1.16518654, 1.25357797, 1.22611968]

View file

@ -293,14 +293,14 @@ def test_run(run_in_tmpdir, pin_model_attributes, mpi_intracomm):
# C API execution modes and ensuring they give the same result.
sp_path = test_model.run(output=False)
with openmc.StatePoint(sp_path) as sp:
cli_keff = sp.k_combined
cli_keff = sp.keff
cli_flux = sp.get_tally(id=1).get_values(scores=['flux'])[0, 0, 0]
cli_fiss = sp.get_tally(id=1).get_values(scores=['fission'])[0, 0, 0]
test_model.init_lib(output=False, intracomm=mpi_intracomm)
sp_path = test_model.run(output=False)
with openmc.StatePoint(sp_path) as sp:
lib_keff = sp.k_combined
lib_keff = sp.keff
lib_flux = sp.get_tally(id=1).get_values(scores=['flux'])[0, 0, 0]
lib_fiss = sp.get_tally(id=1).get_values(scores=['fission'])[0, 0, 0]

View file

@ -147,7 +147,7 @@ def test_interpolation(model, method, temperature, fission_expected):
assert abs(nu_fission_mean - nu*fission_expected) < 3*nu_fission_unc
# Check that k-effective value matches expected
k = sp.k_combined
k = sp.keff
if isnan(k.s):
assert k.n == pytest.approx(nu*fission_expected)
else:

View file

@ -24,7 +24,7 @@ def get_torus_keff(cls, R, r, center=(0, 0, 0)):
sp_path = model.run()
with openmc.StatePoint(sp_path) as sp:
return sp.k_combined
return sp.keff
@pytest.mark.parametrize("R,r", [(2.1, 2.0), (3.0, 1.0)])