crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Change power to source_rate throughout deplete module

Browse source
This commit is contained in:
Paul Romano 2020-07-15 08:09:37 -05:00
parent 04501e3213
commit e2ae3c7401
5 changed files with 74 additions and 65 deletions
Download patch file Download diff file Expand all files Collapse all files

97
openmc/deplete/abc.py
View file

@ -138,15 +138,15 @@ class TransportOperator(ABC):
self._dilute_initial = value
@abstractmethod
def __call__(self, vec, power):
def __call__(self, vec, source_rate):
"""Runs a simulation.
Parameters
----------
vec : list of numpy.ndarray
Total atoms to be used in function.
power : float
Power of the reactor in [W]
source_rate : float
Power in [W] or source rate in [neutron/sec]
Returns
-------
@ -366,9 +366,9 @@ class NormalizationHelper(ABC):
check_type("nuclides", nuclides, list, str)
self._nuclides = nuclides
def factor(self, power):
def factor(self, source_rate):
# Reduce energy produced from all processes
# J / s / source neutron
# J / source neutron
energy = comm.allreduce(self._energy) * JOULE_PER_EV
# Guard against divide by zero
@ -380,8 +380,9 @@ class NormalizationHelper(ABC):
comm.barrier()
comm.Abort(1)
# Return normalization factor for scaling reaction rates
return power / energy
# Return normalization factor for scaling reaction rates. In this case,
# the source rate is the power in [W], so [W] / [J/src] = [src/s]
return source_rate / energy
class FissionYieldHelper(ABC):
@ -619,12 +620,16 @@ class Integrator(ABC):
indicates potentially different power levels for each timestep.
For a 2D problem, the power can be given in [W/cm] as long
as the "volume" assigned to a depletion material is actually
an area in [cm^2]. Either ``power`` or ``power_density`` must be
specified.
an area in [cm^2]. Either ``power``, ``power_density``, or
``source_rates`` must be specified.
power_density : float or iterable of float, optional
Power density of the reactor in [W/gHM]. It is multiplied by
initial heavy metal inventory to get total power if ``power``
is not speficied.
source_rates : iterable of float
Source rate in [neutrons/sec] for each interval in :attr:`timesteps`
.. versionadded:: 0.12.1
timestep_units : {'s', 'min', 'h', 'd', 'MWd/kg'}
Units for values specified in the `timesteps` argument. 's' means
seconds, 'min' means minutes, 'h' means hours, and 'MWd/kg' indicates
@ -650,8 +655,9 @@ class Integrator(ABC):
Depletion chain
timesteps : iterable of float
Size of each depletion interval in [s]
power : iterable of float
Power of the reactor in [W] for each interval in :attr:`timesteps`
source_rates : iterable of float
Source rate in [W] or [neutrons/sec] for each interval in
:attr:`timesteps`
solver : callable
Function that will solve the Bateman equations
:math:`\frac{\partial}{\partial t}\vec{n} = A_i\vec{n}_i` with a step
@ -672,7 +678,7 @@ class Integrator(ABC):
"""
def __init__(self, operator, timesteps, power=None, power_density=None,
timestep_units='s', solver="cram48"):
source_rates=None, timestep_units='s', solver="cram48"):
# Check number of stages previously used
if operator.prev_res is not None:
res = operator.prev_res[-1]
@ -686,22 +692,25 @@ class Integrator(ABC):
self.operator = operator
self.chain = operator.chain
# Determine power and normalize units to W
if power is None:
if power_density is None:
raise ValueError("Either power or power density must be set")
# Determine source rate and normalize units to W in using power
if power is not None:
source_rates = power
elif power_density is not None:
if not isinstance(power_density, Iterable):
power = power_density * operator.heavy_metal
source_rates = power_density * operator.heavy_metal
else:
power = [p*operator.heavy_metal for p in power_density]
if not isinstance(power, Iterable):
# Ensure that power is single value if that is the case
power = [power] * len(timesteps)
source_rates = [p*operator.heavy_metal for p in power_density]
elif source_rates is None:
raise ValueError("Either power, power_density, or source_rates must be set")
if len(power) != len(timesteps):
if not isinstance(source_rates, Iterable):
# Ensure that rate is single value if that is the case
source_rates = [source_rates] * len(timesteps)
if len(source_rates) != len(timesteps):
raise ValueError(
"Number of time steps ({}) != number of powers ({})".format(
len(timesteps), len(power)))
len(timesteps), len(source_rates)))
# Get list of times / units
if isinstance(timesteps[0], Iterable):
@ -712,13 +721,13 @@ class Integrator(ABC):
# Determine number of seconds for each timestep
seconds = []
for timestep, unit, watts in zip(times, units, power):
for timestep, unit, rate in zip(times, units, source_rates):
# Make sure values passed make sense
check_type('timestep', timestep, Real)
check_greater_than('timestep', timestep, 0.0, False)
check_type('timestep units', unit, str)
check_type('power', watts, Real)
check_greater_than('power', watts, 0.0, True)
check_type('source rate', rate, Real)
check_greater_than('source rate', rate, 0.0, True)
if unit in ('s', 'sec'):
seconds.append(timestep)
@ -731,13 +740,13 @@ class Integrator(ABC):
elif unit.lower() == 'mwd/kg':
watt_days_per_kg = 1e6*timestep
kilograms = 1e-3*operator.heavy_metal
days = watt_days_per_kg * kilograms / watts
days = watt_days_per_kg * kilograms / rate
seconds.append(days*_SECONDS_PER_DAY)
else:
raise ValueError("Invalid timestep unit '{}'".format(unit))
self.timesteps = asarray(seconds)
self.power = asarray(power)
self.source_rates = asarray(source_rates)
if isinstance(solver, str):
# Delay importing of cram module, which requires this file
@ -830,22 +839,22 @@ class Integrator(ABC):
"""
def __iter__(self):
"""Return pairs of time steps in [s] and powers in [W]"""
return zip(self.timesteps, self.power)
"""Return pair of time step in [s] and source rate in [W] or [neutrons/sec]"""
return zip(self.timesteps, self.source_rates)
def __len__(self):
"""Return integer number of depletion intervals"""
return len(self.timesteps)
def _get_bos_data_from_operator(self, step_index, step_power, bos_conc):
def _get_bos_data_from_operator(self, step_index, source_rate, bos_conc):
"""Get beginning of step concentrations, reaction rates from Operator
"""
x = deepcopy(bos_conc)
res = self.operator(x, step_power)
res = self.operator(x, source_rate)
self.operator.write_bos_data(step_index + self._i_res)
return x, res
def _get_bos_data_from_restart(self, step_index, step_power, bos_conc):
def _get_bos_data_from_restart(self, step_index, source_rate, bos_conc):
"""Get beginning of step concentrations, reaction rates from restart"""
res = self.operator.prev_res[-1]
# Depletion methods expect list of arrays
@ -853,8 +862,8 @@ class Integrator(ABC):
rates = res.rates[0]
k = ufloat(res.k[0, 0], res.k[0, 1])
# Scale rates by ratio of powers
rates *= step_power / res.power[0]
# Scale reaction rates by ratio of source rates
rates *= source_rate / res.source_rate[0]
return bos_conc, OperatorResult(k, rates)
def _get_start_data(self):
@ -868,12 +877,12 @@ class Integrator(ABC):
with self.operator as conc:
t, self._i_res = self._get_start_data()
for i, (dt, p) in enumerate(self):
for i, (dt, source_rate) in enumerate(self):
if i > 0 or self.operator.prev_res is None:
conc, res = self._get_bos_data_from_operator(i, p, conc)
conc, res = self._get_bos_data_from_operator(i, source_rate, conc)
else:
conc, res = self._get_bos_data_from_restart(i, p, conc)
proc_time, conc_list, res_list = self(conc, res.rates, dt, p, i)
conc, res = self._get_bos_data_from_restart(i, source_rate, conc)
proc_time, conc_list, res_list = self(conc, res.rates, dt, source_rate, i)
# Insert BOS concentration, transport results
conc_list.insert(0, conc)
@ -883,14 +892,14 @@ class Integrator(ABC):
conc = conc_list.pop()
Results.save(self.operator, conc_list, res_list, [t, t + dt],
p, self._i_res + i, proc_time)
source_rate, self._i_res + i, proc_time)
t += dt
# Final simulation
res_list = [self.operator(conc, p)]
res_list = [self.operator(conc, source_rate)]
Results.save(self.operator, [conc], res_list, [t, t],
p, self._i_res + len(self), proc_time)
source_rate, self._i_res + len(self), proc_time)
self.operator.write_bos_data(len(self) + self._i_res)
@ -976,8 +985,8 @@ class SIIntegrator(Integrator):
check_type("n_steps", n_steps, Integral)
check_greater_than("n_steps", n_steps, 0)
super().__init__(
operator, timesteps, power, power_density, timestep_units,
solver=solver)
operator, timesteps, power, power_density,
timestep_units=timestep_units, solver=solver)
self.n_steps = n_steps
def _get_bos_data_from_operator(self, step_index, step_power, bos_conc):

8
openmc/deplete/operator.py
View file

@ -604,7 +604,7 @@ class Operator(TransportOperator):
nuc_list = comm.bcast(nuc_list)
return [nuc for nuc in nuc_list if nuc in self.chain]
def _unpack_tallies_and_normalize(self, power):
def _unpack_tallies_and_normalize(self, source_rate):
"""Unpack tallies from OpenMC and return an operator result
This method uses OpenMC's C API bindings to determine the k-effective
@ -614,8 +614,8 @@ class Operator(TransportOperator):
Parameters
----------
power : float
Power of the reactor in [W]
source_rate : float
Power in [W] or source rate in [neutrons/sec]
Returns
-------
@ -677,7 +677,7 @@ class Operator(TransportOperator):
rates[i] = self._rate_helper.divide_by_adens(number)
# Scale reaction rates to obtain units of reactions/sec
rates *= self._normalization_helper.factor(power)
rates *= self._normalization_helper.factor(source_rate)
# Store new fission yields on the chain
self.chain.fission_yields = fission_yields

32
openmc/deplete/results.py
View file

@ -27,8 +27,8 @@ class Results:
Eigenvalue and uncertainty for each substep.
time : list of float
Time at beginning, end of step, in seconds.
power : float
Power during time step, in Watts
source_rate : float
Source rate during timestep in [W] or [neutrons/sec]
n_mat : int
Number of mats.
n_nuc : int
@ -57,7 +57,7 @@ class Results:
def __init__(self):
self.k = None
self.time = None
self.power = None
self.source_rate = None
self.rates = None
self.volume = None
self.proc_time = None
@ -177,7 +177,7 @@ class Results:
new.mat_to_ind = {mat: idx for (idx, mat) in enumerate(local_materials)}
# Direct transfer
direct_attrs = ("time", "k", "power", "nuc_to_ind",
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))
@ -287,7 +287,7 @@ class Results:
handle.create_dataset("time", (1, 2), maxshape=(None, 2), dtype='float64')
handle.create_dataset("power", (1, n_stages), maxshape=(None, n_stages),
handle.create_dataset("source_rate", (1, n_stages), maxshape=(None, n_stages),
dtype='float64')
handle.create_dataset(
@ -320,7 +320,7 @@ class Results:
rxn_dset = handle["/reaction rates"]
eigenvalues_dset = handle["/eigenvalues"]
time_dset = handle["/time"]
power_dset = handle["/power"]
source_rate_dset = handle["/source_rate"]
proc_time_dset = handle["/depletion time"]
# Get number of results stored
@ -346,9 +346,9 @@ class Results:
time_shape[0] = new_shape
time_dset.resize(time_shape)
power_shape = list(power_dset.shape)
power_shape[0] = new_shape
power_dset.resize(power_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
@ -371,7 +371,7 @@ class Results:
eigenvalues_dset[index, i] = self.k[i]
if comm.rank == 0:
time_dset[index] = self.time
power_dset[index] = self.power
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)
@ -394,12 +394,12 @@ class Results:
number_dset = handle["/number"]
eigenvalues_dset = handle["/eigenvalues"]
time_dset = handle["/time"]
power_dset = handle["/power"]
source_rate_dset = handle["/source_rate"]
results.data = number_dset[step, :, :, :]
results.k = eigenvalues_dset[step, :]
results.time = time_dset[step, :]
results.power = power_dset[step, :]
results.source_rate = source_rate_dset[step, :]
if "depletion time" in handle:
proc_time_dset = handle["/depletion time"]
@ -444,7 +444,7 @@ class Results:
return results
@staticmethod
def save(op, x, op_results, t, power, step_ind, proc_time=None):
def save(op, x, op_results, t, source_rate, step_ind, proc_time=None):
"""Creates and writes depletion results to disk
Parameters
@ -457,8 +457,8 @@ class Results:
Results of applying transport operator
t : list of float
Time indices.
power : float
Power during time step
source_rate : float
Source rate during time step in [W] or [neutrons/sec]
step_ind : int
Step index.
proc_time : float or None
@ -485,7 +485,7 @@ class Results:
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.power = power
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)

BIN
tests/regression_tests/deplete/test_reference.h5
View file

Binary file not shown.

2
tests/unit_tests/test_deplete_integrator.py
View file

@ -130,7 +130,7 @@ def test_bad_integrator_inputs():
timesteps = [1]
# No power nor power density given
with pytest.raises(ValueError, match="Either power or power density"):
with pytest.raises(ValueError, match="Either power"):
PredictorIntegrator(op, timesteps)
# Length of power != length time