mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Merge pull request #1481 from paulromano/depletion-by-burnup
Support different timestep units for depletion
This commit is contained in:
Andrew Johnson
GitHub
commit
9d883cc29b
6 changed files with 295 additions and 170 deletions
|
|
@ -33,12 +33,11 @@ material compositions over time. Each method appears as a different class.
|
|||
For example, :class:`openmc.deplete.CECMIntegrator` runs a depletion calculation
|
||||
using the CE/CM algorithm (deplete over a timestep using the middle-of-step
|
||||
reaction rates). An instance of :class:`openmc.deplete.Operator` is passed to
|
||||
one of these functions along with the power level and timesteps::
|
||||
one of these functions along with the timesteps and power level::
|
||||
|
||||
power = 1200.0e6
|
||||
days = 24*60*60
|
||||
timesteps = [10.0*days, 10.0*days, 10.0*days]
|
||||
openmc.deplete.CECMIntegrator(op, power, timesteps).integrate()
|
||||
power = 1200.0e6 # watts
|
||||
timesteps = [10.0, 10.0, 10.0] # days
|
||||
openmc.deplete.CECMIntegrator(op, timesteps, power, timestep_units='d').integrate()
|
||||
|
||||
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
|
||||
|
|
@ -67,7 +66,7 @@ the energy deposited during a transport calculation will be lower than expected.
|
|||
This causes the reaction rates to be over-adjusted to hit the user-specific power,
|
||||
or power density, leading to an over-depletion of burnable materials.
|
||||
|
||||
There are some remedies. First, the fission Q values can be directly set in a
|
||||
There are some remedies. First, the fission Q values can be directly set in a
|
||||
variety of ways. This requires knowing what the total fission energy release should
|
||||
be, including indirect components. Some examples are provided below::
|
||||
|
||||
|
|
@ -99,11 +98,11 @@ Local Spectra and Repeated Materials
|
|||
------------------------------------
|
||||
|
||||
It is not uncommon to explicitly create a single burnable material across many locations.
|
||||
From a pure transport perspective, there is nothing wrong with creating a single
|
||||
From a pure transport perspective, there is nothing wrong with creating a single
|
||||
3.5 wt.% enriched fuel ``fuel_3``, and placing that fuel in every fuel pin in an assembly
|
||||
or even full core problem. This certainly expedites the model making process, but can pose
|
||||
issues with depletion.
|
||||
Under this setup, :mod:`openmc.deplete` will deplete a single ``fuel_3`` material using
|
||||
issues with depletion.
|
||||
Under this setup, :mod:`openmc.deplete` will deplete a single ``fuel_3`` material using
|
||||
a single set of reaction rates, and produce a single new composition for the next time
|
||||
step. This can be problematic if the same ``fuel_3`` is used in very different regions
|
||||
of the problem.
|
||||
|
|
|
|||
|
|
@ -1,25 +1,7 @@
|
|||
import openmc
|
||||
import openmc.deplete
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
###############################################################################
|
||||
# Simulation Input File Parameters
|
||||
###############################################################################
|
||||
|
||||
# OpenMC simulation parameters
|
||||
batches = 100
|
||||
inactive = 10
|
||||
particles = 1000
|
||||
|
||||
# Depletion simulation parameters
|
||||
time_step = 1*24*60*60 # s
|
||||
final_time = 5*24*60*60 # s
|
||||
time_steps = np.full(final_time // time_step, time_step)
|
||||
|
||||
chain_file = './chain_simple.xml'
|
||||
power = 174 # W/cm, for 2D simulations only (use W for 3D)
|
||||
|
||||
###############################################################################
|
||||
# Load previous simulation results
|
||||
###############################################################################
|
||||
|
|
@ -37,31 +19,34 @@ previous_results = openmc.deplete.ResultsList("depletion_results.h5")
|
|||
###############################################################################
|
||||
|
||||
# Instantiate a Settings object, set all runtime parameters
|
||||
settings_file = openmc.Settings()
|
||||
settings_file.batches = batches
|
||||
settings_file.inactive = inactive
|
||||
settings_file.particles = particles
|
||||
settings = openmc.Settings()
|
||||
settings.batches = 100
|
||||
settings.inactive = 10
|
||||
settings.particles = 10000
|
||||
|
||||
# Create an initial uniform spatial source distribution over fissionable zones
|
||||
bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
|
||||
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
|
||||
settings_file.source = openmc.source.Source(space=uniform_dist)
|
||||
settings.source = openmc.source.Source(space=uniform_dist)
|
||||
|
||||
entropy_mesh = openmc.RegularMesh()
|
||||
entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]
|
||||
entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50]
|
||||
entropy_mesh.dimension = [10, 10, 1]
|
||||
settings_file.entropy_mesh = entropy_mesh
|
||||
settings.entropy_mesh = entropy_mesh
|
||||
|
||||
###############################################################################
|
||||
# Initialize and run depletion calculation
|
||||
###############################################################################
|
||||
|
||||
op = openmc.deplete.Operator(geometry, settings_file, chain_file,
|
||||
previous_results)
|
||||
# Create depletion "operator"
|
||||
chain_file = './chain_simple.xml'
|
||||
op = openmc.deplete.Operator(geometry, settings, chain_file, previous_results)
|
||||
|
||||
# Perform simulation using the predictor algorithm
|
||||
integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)
|
||||
time_steps = [1.0, 1.0, 1.0, 1.0, 1.0] # days
|
||||
power = 174 # W/cm, for 2D simulations only (use W for 3D)
|
||||
integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d')
|
||||
integrator.integrate()
|
||||
|
||||
###############################################################################
|
||||
|
|
@ -77,27 +62,28 @@ time, keff = results.get_eigenvalue()
|
|||
# Obtain U235 concentration as a function of time
|
||||
time, n_U235 = results.get_atoms('1', 'U235')
|
||||
|
||||
# Obtain Xe135 absorption as a function of time
|
||||
time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
|
||||
# Obtain Xe135 capture reaction rate as a function of time
|
||||
time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
|
||||
|
||||
###############################################################################
|
||||
# Generate plots
|
||||
###############################################################################
|
||||
|
||||
days = 24*60*60
|
||||
plt.figure()
|
||||
plt.plot(time/(24*60*60), keff, label="K-effective")
|
||||
plt.plot(time/days, keff, label="K-effective")
|
||||
plt.xlabel("Time (days)")
|
||||
plt.ylabel("Keff")
|
||||
plt.show()
|
||||
|
||||
plt.figure()
|
||||
plt.plot(time/(24*60*60), n_U235, label="U 235")
|
||||
plt.plot(time/days, n_U235, label="U 235")
|
||||
plt.xlabel("Time (days)")
|
||||
plt.ylabel("n U5 (-)")
|
||||
plt.show()
|
||||
|
||||
plt.figure()
|
||||
plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption")
|
||||
plt.plot(time/days, Xe_capture, label="Xe135 capture")
|
||||
plt.xlabel("Time (days)")
|
||||
plt.ylabel("RR (-)")
|
||||
plt.show()
|
||||
|
|
|
|||
|
|
@ -1,47 +1,31 @@
|
|||
from math import pi
|
||||
|
||||
import openmc
|
||||
import openmc.deplete
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
###############################################################################
|
||||
# Simulation Input File Parameters
|
||||
###############################################################################
|
||||
|
||||
# OpenMC simulation parameters
|
||||
batches = 100
|
||||
inactive = 10
|
||||
particles = 1000
|
||||
|
||||
# Depletion simulation parameters
|
||||
time_step = 1*24*60*60 # s
|
||||
final_time = 5*24*60*60 # s
|
||||
time_steps = np.full(final_time // time_step, time_step)
|
||||
chain_file = './chain_simple.xml'
|
||||
power = 174 # W/cm, for 2D simulations only (use W for 3D)
|
||||
|
||||
###############################################################################
|
||||
# Define materials
|
||||
###############################################################################
|
||||
|
||||
# Instantiate some Materials and register the appropriate Nuclides
|
||||
uo2 = openmc.Material(material_id=1, name='UO2 fuel at 2.4% wt enrichment')
|
||||
uo2 = openmc.Material(name='UO2 fuel at 2.4% wt enrichment')
|
||||
uo2.set_density('g/cm3', 10.29769)
|
||||
uo2.add_element('U', 1., enrichment=2.4)
|
||||
uo2.add_element('O', 2.)
|
||||
uo2.depletable = True
|
||||
|
||||
helium = openmc.Material(material_id=2, name='Helium for gap')
|
||||
helium = openmc.Material(name='Helium for gap')
|
||||
helium.set_density('g/cm3', 0.001598)
|
||||
helium.add_element('He', 2.4044e-4)
|
||||
|
||||
zircaloy = openmc.Material(material_id=3, name='Zircaloy 4')
|
||||
zircaloy = openmc.Material(name='Zircaloy 4')
|
||||
zircaloy.set_density('g/cm3', 6.55)
|
||||
zircaloy.add_element('Sn', 0.014 , 'wo')
|
||||
zircaloy.add_element('Sn', 0.014, 'wo')
|
||||
zircaloy.add_element('Fe', 0.00165, 'wo')
|
||||
zircaloy.add_element('Cr', 0.001 , 'wo')
|
||||
zircaloy.add_element('Cr', 0.001, 'wo')
|
||||
zircaloy.add_element('Zr', 0.98335, 'wo')
|
||||
|
||||
borated_water = openmc.Material(material_id=4, name='Borated water')
|
||||
borated_water = openmc.Material(name='Borated water')
|
||||
borated_water.set_density('g/cm3', 0.740582)
|
||||
borated_water.add_element('B', 4.0e-5)
|
||||
borated_water.add_element('H', 5.0e-2)
|
||||
|
|
@ -52,87 +36,62 @@ borated_water.add_s_alpha_beta('c_H_in_H2O')
|
|||
# Create geometry
|
||||
###############################################################################
|
||||
|
||||
# Instantiate ZCylinder surfaces
|
||||
fuel_or = openmc.ZCylinder(surface_id=1, x0=0, y0=0, r=0.39218, name='Fuel OR')
|
||||
clad_ir = openmc.ZCylinder(surface_id=2, x0=0, y0=0, r=0.40005, name='Clad IR')
|
||||
clad_or = openmc.ZCylinder(surface_id=3, x0=0, y0=0, r=0.45720, name='Clad OR')
|
||||
left = openmc.XPlane(surface_id=4, x0=-0.62992, name='left')
|
||||
right = openmc.XPlane(surface_id=5, x0=0.62992, name='right')
|
||||
bottom = openmc.YPlane(surface_id=6, y0=-0.62992, name='bottom')
|
||||
top = openmc.YPlane(surface_id=7, y0=0.62992, name='top')
|
||||
# Define surfaces
|
||||
pitch = 1.25984
|
||||
fuel_or = openmc.ZCylinder(r=0.39218, name='Fuel OR')
|
||||
clad_ir = openmc.ZCylinder(r=0.40005, name='Clad IR')
|
||||
clad_or = openmc.ZCylinder(r=0.45720, name='Clad OR')
|
||||
box = openmc.model.rectangular_prism(pitch, pitch, boundary_type='reflective')
|
||||
|
||||
left.boundary_type = 'reflective'
|
||||
right.boundary_type = 'reflective'
|
||||
top.boundary_type = 'reflective'
|
||||
bottom.boundary_type = 'reflective'
|
||||
# Define cells
|
||||
fuel = openmc.Cell(fill=uo2, region=-fuel_or)
|
||||
gap = openmc.Cell(fill=helium, region=+fuel_or & -clad_ir)
|
||||
clad = openmc.Cell(fill=zircaloy, region=+clad_ir & -clad_or)
|
||||
water = openmc.Cell(fill=borated_water, region=+clad_or & box)
|
||||
|
||||
# Instantiate Cells
|
||||
fuel = openmc.Cell(cell_id=1, name='cell 1')
|
||||
gap = openmc.Cell(cell_id=2, name='cell 2')
|
||||
clad = openmc.Cell(cell_id=3, name='cell 3')
|
||||
water = openmc.Cell(cell_id=4, name='cell 4')
|
||||
|
||||
# Use surface half-spaces to define regions
|
||||
fuel.region = -fuel_or
|
||||
gap.region = +fuel_or & -clad_ir
|
||||
clad.region = +clad_ir & -clad_or
|
||||
water.region = +clad_or & +left & -right & +bottom & -top
|
||||
|
||||
# Register Materials with Cells
|
||||
fuel.fill = uo2
|
||||
gap.fill = helium
|
||||
clad.fill = zircaloy
|
||||
water.fill = borated_water
|
||||
|
||||
# Instantiate Universe
|
||||
root = openmc.Universe(universe_id=0, name='root universe')
|
||||
|
||||
# Register Cells with Universe
|
||||
root.add_cells([fuel, gap, clad, water])
|
||||
|
||||
# Instantiate a Geometry, register the root Universe
|
||||
geometry = openmc.Geometry(root)
|
||||
# Define overall geometry
|
||||
geometry = openmc.Geometry([fuel, gap, clad, water])
|
||||
|
||||
###############################################################################
|
||||
# Set volumes of depletable materials
|
||||
###############################################################################
|
||||
|
||||
# Compute cell areas
|
||||
area = {}
|
||||
area[fuel] = np.pi * fuel_or.coefficients['r'] ** 2
|
||||
|
||||
# Set materials volume for depletion. Set to an area for 2D simulations
|
||||
uo2.volume = area[fuel]
|
||||
# Set material volume for depletion. For 2D simulations, this should be an area.
|
||||
uo2.volume = pi * fuel_or.r**2
|
||||
|
||||
###############################################################################
|
||||
# Transport calculation settings
|
||||
###############################################################################
|
||||
|
||||
# Instantiate a Settings object, set all runtime parameters, and export to XML
|
||||
settings_file = openmc.Settings()
|
||||
settings_file.batches = batches
|
||||
settings_file.inactive = inactive
|
||||
settings_file.particles = particles
|
||||
settings = openmc.Settings()
|
||||
settings.batches = 100
|
||||
settings.inactive = 10
|
||||
settings.particles = 1000
|
||||
|
||||
# Create an initial uniform spatial source distribution over fissionable zones
|
||||
bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
|
||||
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
|
||||
settings_file.source = openmc.source.Source(space=uniform_dist)
|
||||
settings.source = openmc.source.Source(space=uniform_dist)
|
||||
|
||||
entropy_mesh = openmc.RegularMesh()
|
||||
entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]
|
||||
entropy_mesh.upper_right = [0.39218, 0.39218, 1.e50]
|
||||
entropy_mesh.dimension = [10, 10, 1]
|
||||
settings_file.entropy_mesh = entropy_mesh
|
||||
settings.entropy_mesh = entropy_mesh
|
||||
|
||||
###############################################################################
|
||||
# Initialize and run depletion calculation
|
||||
###############################################################################
|
||||
|
||||
op = openmc.deplete.Operator(geometry, settings_file, chain_file)
|
||||
# Create depletion "operator"
|
||||
chain_file = './chain_simple.xml'
|
||||
op = openmc.deplete.Operator(geometry, settings, chain_file)
|
||||
|
||||
# Perform simulation using the predictor algorithm
|
||||
integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power)
|
||||
time_steps = [1.0, 1.0, 1.0, 1.0, 1.0] # days
|
||||
power = 174 # W/cm, for 2D simulations only (use W for 3D)
|
||||
integrator = openmc.deplete.PredictorIntegrator(op, time_steps, power, timestep_units='d')
|
||||
integrator.integrate()
|
||||
|
||||
###############################################################################
|
||||
|
|
@ -148,27 +107,28 @@ time, keff = results.get_eigenvalue()
|
|||
# Obtain U235 concentration as a function of time
|
||||
time, n_U235 = results.get_atoms('1', 'U235')
|
||||
|
||||
# Obtain Xe135 absorption as a function of time
|
||||
time, Xe_gam = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
|
||||
# Obtain Xe135 capture reaction rate as a function of time
|
||||
time, Xe_capture = results.get_reaction_rate('1', 'Xe135', '(n,gamma)')
|
||||
|
||||
###############################################################################
|
||||
# Generate plots
|
||||
###############################################################################
|
||||
|
||||
days = 24*60*60
|
||||
plt.figure()
|
||||
plt.plot(time/(24*60*60), keff, label="K-effective")
|
||||
plt.plot(time/days, keff, label="K-effective")
|
||||
plt.xlabel("Time (days)")
|
||||
plt.ylabel("Keff")
|
||||
plt.show()
|
||||
|
||||
plt.figure()
|
||||
plt.plot(time/(24*60*60), n_U235, label="U 235")
|
||||
plt.plot(time/days, n_U235, label="U235")
|
||||
plt.xlabel("Time (days)")
|
||||
plt.ylabel("n U5 (-)")
|
||||
plt.show()
|
||||
|
||||
plt.figure()
|
||||
plt.plot(time/(24*60*60), Xe_gam, label="Xe135 absorption")
|
||||
plt.plot(time/days, Xe_capture, label="Xe135 capture")
|
||||
plt.xlabel("Time (days)")
|
||||
plt.ylabel("RR (-)")
|
||||
plt.show()
|
||||
|
|
|
|||
|
|
@ -31,6 +31,10 @@ __all__ = [
|
|||
"Integrator", "SIIntegrator", "DepSystemSolver"]
|
||||
|
||||
|
||||
_SECONDS_PER_MINUTE = 60
|
||||
_SECONDS_PER_HOUR = 60*60
|
||||
_SECONDS_PER_DAY = 24*60*60
|
||||
|
||||
OperatorResult = namedtuple('OperatorResult', ['k', 'rates'])
|
||||
OperatorResult.__doc__ = """\
|
||||
Result of applying transport operator
|
||||
|
|
@ -597,9 +601,11 @@ class Integrator(ABC):
|
|||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
Operator to perform transport simulations
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
timesteps : iterable of float or iterable of tuple
|
||||
Array of timesteps. Note that values are not cumulative. The units are
|
||||
specified by the `timestep_units` argument when `timesteps` is an
|
||||
iterable of float. Alternatively, units can be specified for each step
|
||||
by passing an iterable of (value, unit) tuples.
|
||||
power : float or iterable of float, optional
|
||||
Power of the reactor in [W]. A single value indicates that
|
||||
the power is constant over all timesteps. An iterable
|
||||
|
|
@ -612,6 +618,11 @@ class Integrator(ABC):
|
|||
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.
|
||||
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
|
||||
that the values are given in burnup (MW-d of energy deposited per
|
||||
kilogram of initial heavy metal).
|
||||
|
||||
Attributes
|
||||
----------
|
||||
|
|
@ -625,7 +636,8 @@ class Integrator(ABC):
|
|||
Power of the reactor in [W] for each interval in :attr:`timesteps`
|
||||
"""
|
||||
|
||||
def __init__(self, operator, timesteps, power=None, power_density=None):
|
||||
def __init__(self, operator, timesteps, power=None, power_density=None,
|
||||
timestep_units='s'):
|
||||
# Check number of stages previously used
|
||||
if operator.prev_res is not None:
|
||||
res = operator.prev_res[-1]
|
||||
|
|
@ -638,27 +650,59 @@ class Integrator(ABC):
|
|||
self._num_stages))
|
||||
self.operator = operator
|
||||
self.chain = operator.chain
|
||||
if not isinstance(timesteps, Iterable):
|
||||
self.timesteps = [timesteps]
|
||||
else:
|
||||
self.timesteps = timesteps
|
||||
|
||||
# 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")
|
||||
if not isinstance(power_density, Iterable):
|
||||
power = power_density * operator.heavy_metal
|
||||
else:
|
||||
power = [p * operator.heavy_metal for p in power_density]
|
||||
|
||||
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(self.timesteps)
|
||||
elif len(power) != len(self.timesteps):
|
||||
raise ValueError(
|
||||
"Number of time steps != number of powers. {} vs {}".format(
|
||||
len(self.timesteps), len(power)))
|
||||
power = [power] * len(timesteps)
|
||||
|
||||
self.power = power
|
||||
if len(power) != len(timesteps):
|
||||
raise ValueError(
|
||||
"Number of time steps ({}) != number of powers ({})".format(
|
||||
len(timesteps), len(power)))
|
||||
|
||||
# Get list of times / units
|
||||
if isinstance(timesteps[0], Iterable):
|
||||
times, units = zip(*timesteps)
|
||||
else:
|
||||
times = timesteps
|
||||
units = [timestep_units] * len(timesteps)
|
||||
|
||||
# Determine number of seconds for each timestep
|
||||
seconds = []
|
||||
for time, unit, watts in zip(times, units, power):
|
||||
# Make sure values passed make sense
|
||||
check_type('timestep', time, Real)
|
||||
check_greater_than('timestep', time, 0.0, False)
|
||||
check_type('timestep units', unit, str)
|
||||
check_type('power', watts, Real)
|
||||
check_greater_than('power', watts, 0.0, True)
|
||||
|
||||
if unit in ('s', 'sec'):
|
||||
seconds.append(time)
|
||||
elif unit in ('min', 'minute'):
|
||||
seconds.append(time*_SECONDS_PER_MINUTE)
|
||||
elif unit in ('h', 'hr', 'hour'):
|
||||
seconds.append(time*_SECONDS_PER_HOUR)
|
||||
elif unit in ('d', 'day'):
|
||||
seconds.append(time*_SECONDS_PER_DAY)
|
||||
elif unit.lower() == 'mwd/kg':
|
||||
watt_days_per_kg = 1e6*time
|
||||
kilograms = 1e-3*operator.heavy_metal
|
||||
days = watt_days_per_kg * kilograms / watts
|
||||
seconds.append(days*_SECONDS_PER_DAY)
|
||||
else:
|
||||
raise ValueError("Invalid timestep unit '{}'".format(unit))
|
||||
|
||||
self.timesteps = asarray(seconds)
|
||||
self.power = asarray(power)
|
||||
|
||||
@abstractmethod
|
||||
def __call__(self, conc, rates, dt, power, i):
|
||||
|
|
@ -772,9 +816,11 @@ class SIIntegrator(Integrator):
|
|||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
The operator object to simulate on.
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
timesteps : iterable of float or iterable of tuple
|
||||
Array of timesteps. Note that values are not cumulative. The units are
|
||||
specified by the `timestep_units` argument when `timesteps` is an
|
||||
iterable of float. Alternatively, units can be specified for each step
|
||||
by passing an iterable of (value, unit) tuples.
|
||||
power : float or iterable of float, optional
|
||||
Power of the reactor in [W]. A single value indicates that
|
||||
the power is constant over all timesteps. An iterable
|
||||
|
|
@ -787,6 +833,11 @@ class SIIntegrator(Integrator):
|
|||
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.
|
||||
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
|
||||
that the values are given in burnup (MW-d of energy deposited per
|
||||
kilogram of initial heavy metal).
|
||||
n_steps : int, optional
|
||||
Number of stochastic iterations per depletion interval.
|
||||
Must be greater than zero. Default : 10
|
||||
|
|
@ -805,10 +856,10 @@ class SIIntegrator(Integrator):
|
|||
Number of stochastic iterations per depletion interval
|
||||
"""
|
||||
def __init__(self, operator, timesteps, power=None, power_density=None,
|
||||
n_steps=10):
|
||||
timestep_units='s', n_steps=10):
|
||||
check_type("n_steps", n_steps, Integral)
|
||||
check_greater_than("n_steps", n_steps, 0)
|
||||
super().__init__(operator, timesteps, power, power_density)
|
||||
super().__init__(operator, timesteps, power, power_density, timestep_units)
|
||||
self.n_steps = n_steps
|
||||
|
||||
def _get_bos_data_from_operator(self, step_index, step_power, bos_conc):
|
||||
|
|
|
|||
|
|
@ -31,9 +31,11 @@ class PredictorIntegrator(Integrator):
|
|||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
Operator to perform transport simulations
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
timesteps : iterable of float or iterable of tuple
|
||||
Array of timesteps. Note that values are not cumulative. The units are
|
||||
specified by the `timestep_units` argument when `timesteps` is an
|
||||
iterable of float. Alternatively, units can be specified for each step
|
||||
by passing an iterable of (value, unit) tuples.
|
||||
power : float or iterable of float, optional
|
||||
Power of the reactor in [W]. A single value indicates that
|
||||
the power is constant over all timesteps. An iterable
|
||||
|
|
@ -46,6 +48,11 @@ class PredictorIntegrator(Integrator):
|
|||
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.
|
||||
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
|
||||
that the values are given in burnup (MW-d of energy deposited per
|
||||
kilogram of initial heavy metal).
|
||||
|
||||
Attributes
|
||||
----------
|
||||
|
|
@ -113,9 +120,11 @@ class CECMIntegrator(Integrator):
|
|||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
Operator to perform transport simulations
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
timesteps : iterable of float or iterable of tuple
|
||||
Array of timesteps. Note that values are not cumulative. The units are
|
||||
specified by the `timestep_units` argument when `timesteps` is an
|
||||
iterable of float. Alternatively, units can be specified for each step
|
||||
by passing an iterable of (value, unit) tuples.
|
||||
power : float or iterable of float, optional
|
||||
Power of the reactor in [W]. A single value indicates that
|
||||
the power is constant over all timesteps. An iterable
|
||||
|
|
@ -128,6 +137,11 @@ class CECMIntegrator(Integrator):
|
|||
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.
|
||||
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
|
||||
that the values are given in burnup (MW-d of energy deposited per
|
||||
kilogram of initial heavy metal).
|
||||
|
||||
Attributes
|
||||
----------
|
||||
|
|
@ -203,9 +217,11 @@ class CF4Integrator(Integrator):
|
|||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
Operator to perform transport simulations
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
timesteps : iterable of float or iterable of tuple
|
||||
Array of timesteps. Note that values are not cumulative. The units are
|
||||
specified by the `timestep_units` argument when `timesteps` is an
|
||||
iterable of float. Alternatively, units can be specified for each step
|
||||
by passing an iterable of (value, unit) tuples.
|
||||
power : float or iterable of float, optional
|
||||
Power of the reactor in [W]. A single value indicates that
|
||||
the power is constant over all timesteps. An iterable
|
||||
|
|
@ -218,6 +234,11 @@ class CF4Integrator(Integrator):
|
|||
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.
|
||||
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
|
||||
that the values are given in burnup (MW-d of energy deposited per
|
||||
kilogram of initial heavy metal).
|
||||
|
||||
Attributes
|
||||
----------
|
||||
|
|
@ -310,9 +331,11 @@ class CELIIntegrator(Integrator):
|
|||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
Operator to perform transport simulations
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
timesteps : iterable of float or iterable of tuple
|
||||
Array of timesteps. Note that values are not cumulative. The units are
|
||||
specified by the `timestep_units` argument when `timesteps` is an
|
||||
iterable of float. Alternatively, units can be specified for each step
|
||||
by passing an iterable of (value, unit) tuples.
|
||||
power : float or iterable of float, optional
|
||||
Power of the reactor in [W]. A single value indicates that
|
||||
the power is constant over all timesteps. An iterable
|
||||
|
|
@ -325,6 +348,11 @@ class CELIIntegrator(Integrator):
|
|||
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.
|
||||
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
|
||||
that the values are given in burnup (MW-d of energy deposited per
|
||||
kilogram of initial heavy metal).
|
||||
|
||||
Attributes
|
||||
----------
|
||||
|
|
@ -404,9 +432,11 @@ class EPCRK4Integrator(Integrator):
|
|||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
Operator to perform transport simulations
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
timesteps : iterable of float or iterable of tuple
|
||||
Array of timesteps. Note that values are not cumulative. The units are
|
||||
specified by the `timestep_units` argument when `timesteps` is an
|
||||
iterable of float. Alternatively, units can be specified for each step
|
||||
by passing an iterable of (value, unit) tuples.
|
||||
power : float or iterable of float, optional
|
||||
Power of the reactor in [W]. A single value indicates that
|
||||
the power is constant over all timesteps. An iterable
|
||||
|
|
@ -419,6 +449,11 @@ class EPCRK4Integrator(Integrator):
|
|||
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.
|
||||
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
|
||||
that the values are given in burnup (MW-d of energy deposited per
|
||||
kilogram of initial heavy metal).
|
||||
|
||||
Attributes
|
||||
----------
|
||||
|
|
@ -518,9 +553,11 @@ class LEQIIntegrator(Integrator):
|
|||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
Operator to perform transport simulations
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
timesteps : iterable of float or iterable of tuple
|
||||
Array of timesteps. Note that values are not cumulative. The units are
|
||||
specified by the `timestep_units` argument when `timesteps` is an
|
||||
iterable of float. Alternatively, units can be specified for each step
|
||||
by passing an iterable of (value, unit) tuples.
|
||||
power : float or iterable of float, optional
|
||||
Power of the reactor in [W]. A single value indicates that
|
||||
the power is constant over all timesteps. An iterable
|
||||
|
|
@ -533,6 +570,11 @@ class LEQIIntegrator(Integrator):
|
|||
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.
|
||||
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
|
||||
that the values are given in burnup (MW-d of energy deposited per
|
||||
kilogram of initial heavy metal).
|
||||
|
||||
Attributes
|
||||
----------
|
||||
|
|
@ -629,9 +671,11 @@ class SICELIIntegrator(SIIntegrator):
|
|||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
The operator object to simulate on.
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
timesteps : iterable of float or iterable of tuple
|
||||
Array of timesteps. Note that values are not cumulative. The units are
|
||||
specified by the `timestep_units` argument when `timesteps` is an
|
||||
iterable of float. Alternatively, units can be specified for each step
|
||||
by passing an iterable of (value, unit) tuples.
|
||||
power : float or iterable of float, optional
|
||||
Power of the reactor in [W]. A single value indicates that
|
||||
the power is constant over all timesteps. An iterable
|
||||
|
|
@ -644,6 +688,11 @@ class SICELIIntegrator(SIIntegrator):
|
|||
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.
|
||||
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
|
||||
that the values are given in burnup (MW-d of energy deposited per
|
||||
kilogram of initial heavy metal).
|
||||
n_steps : int, optional
|
||||
Number of stochastic iterations per depletion interval.
|
||||
Must be greater than zero. Default : 10
|
||||
|
|
@ -730,9 +779,11 @@ class SILEQIIntegrator(SIIntegrator):
|
|||
----------
|
||||
operator : openmc.deplete.TransportOperator
|
||||
The operator object to simulate on.
|
||||
timesteps : iterable of float
|
||||
Array of timesteps in units of [s]. Note that values are not
|
||||
cumulative.
|
||||
timesteps : iterable of float or iterable of tuple
|
||||
Array of timesteps. Note that values are not cumulative. The units are
|
||||
specified by the `timestep_units` argument when `timesteps` is an
|
||||
iterable of float. Alternatively, units can be specified for each step
|
||||
by passing an iterable of (value, unit) tuples.
|
||||
power : float or iterable of float, optional
|
||||
Power of the reactor in [W]. A single value indicates that
|
||||
the power is constant over all timesteps. An iterable
|
||||
|
|
@ -745,6 +796,11 @@ class SILEQIIntegrator(SIIntegrator):
|
|||
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.
|
||||
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
|
||||
that the values are given in burnup (MW-d of energy deposited per
|
||||
kilogram of initial heavy metal).
|
||||
n_steps : int, optional
|
||||
Number of stochastic iterations per depletion interval.
|
||||
Must be greater than zero. Default : 10
|
||||
|
|
|
|||
|
|
@ -7,6 +7,7 @@ will be left unimplemented and testing will be done via regression.
|
|||
"""
|
||||
|
||||
import copy
|
||||
from random import uniform
|
||||
from unittest.mock import MagicMock
|
||||
|
||||
import numpy as np
|
||||
|
|
@ -15,11 +16,24 @@ import pytest
|
|||
|
||||
from openmc.deplete import (
|
||||
ReactionRates, Results, ResultsList, comm, OperatorResult,
|
||||
PredictorIntegrator, SICELIIntegrator)
|
||||
PredictorIntegrator, CECMIntegrator, CF4Integrator, CELIIntegrator,
|
||||
EPCRK4Integrator, LEQIIntegrator, SICELIIntegrator, SILEQIIntegrator)
|
||||
|
||||
from tests import dummy_operator
|
||||
|
||||
|
||||
INTEGRATORS = [
|
||||
PredictorIntegrator,
|
||||
CECMIntegrator,
|
||||
CF4Integrator,
|
||||
CELIIntegrator,
|
||||
EPCRK4Integrator,
|
||||
LEQIIntegrator,
|
||||
SICELIIntegrator,
|
||||
SILEQIIntegrator
|
||||
]
|
||||
|
||||
|
||||
def test_results_save(run_in_tmpdir):
|
||||
"""Test data save module"""
|
||||
|
||||
|
|
@ -105,14 +119,14 @@ def test_results_save(run_in_tmpdir):
|
|||
np.testing.assert_array_equal(res[1].time, t2)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("timesteps", (1, [1]))
|
||||
def test_bad_integrator_inputs(timesteps):
|
||||
def test_bad_integrator_inputs():
|
||||
"""Test failure modes for Integrator inputs"""
|
||||
|
||||
op = MagicMock()
|
||||
op.prev_res = None
|
||||
op.chain = None
|
||||
op.heavy_metal = 1.0
|
||||
timesteps = [1]
|
||||
|
||||
# No power nor power density given
|
||||
with pytest.raises(ValueError, match="Either power or power density"):
|
||||
|
|
@ -159,3 +173,62 @@ def test_integrator(run_in_tmpdir, scheme):
|
|||
dep_time = res.get_depletion_time()
|
||||
assert dep_time.shape == (2, )
|
||||
assert all(dep_time > 0)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("integrator", INTEGRATORS)
|
||||
def test_timesteps(integrator):
|
||||
# Crate fake operator
|
||||
op = MagicMock()
|
||||
op.prev_res = None
|
||||
op.chain = None
|
||||
|
||||
# Set heavy metal mass and power randomly
|
||||
op.heavy_metal = uniform(0, 10000)
|
||||
power = uniform(0, 1e6)
|
||||
|
||||
# Reference timesteps in seconds
|
||||
day = 86400.0
|
||||
ref_timesteps = [1*day, 2*day, 5*day, 10*day]
|
||||
|
||||
# Case 1, timesteps in seconds
|
||||
timesteps = ref_timesteps
|
||||
x = integrator(op, timesteps, power, timestep_units='s')
|
||||
assert np.allclose(x.timesteps, ref_timesteps)
|
||||
|
||||
# Case 2, timesteps in minutes
|
||||
minute = 60
|
||||
timesteps = [t / minute for t in ref_timesteps]
|
||||
x = integrator(op, timesteps, power, timestep_units='min')
|
||||
assert np.allclose(x.timesteps, ref_timesteps)
|
||||
|
||||
# Case 3, timesteps in hours
|
||||
hour = 60*60
|
||||
timesteps = [t / hour for t in ref_timesteps]
|
||||
x = integrator(op, timesteps, power, timestep_units='h')
|
||||
assert np.allclose(x.timesteps, ref_timesteps)
|
||||
|
||||
# Case 4, timesteps in days
|
||||
timesteps = [t / day for t in ref_timesteps]
|
||||
x = integrator(op, timesteps, power, timestep_units='d')
|
||||
assert np.allclose(x.timesteps, ref_timesteps)
|
||||
|
||||
# Case 5, timesteps in MWd/kg
|
||||
kilograms = op.heavy_metal / 1000.0
|
||||
days = [t/day for t in ref_timesteps]
|
||||
megawatts = power / 1000000.0
|
||||
burnup = [t * megawatts / kilograms for t in days]
|
||||
x = integrator(op, burnup, power, timestep_units='MWd/kg')
|
||||
assert np.allclose(x.timesteps, ref_timesteps)
|
||||
|
||||
# Case 6, mixed units
|
||||
burnup_per_day = (1e-6*power) / kilograms
|
||||
timesteps = [(burnup_per_day, 'MWd/kg'), (2*day, 's'), (5, 'd'),
|
||||
(10*burnup_per_day, 'MWd/kg')]
|
||||
x = integrator(op, timesteps, power)
|
||||
assert np.allclose(x.timesteps, ref_timesteps)
|
||||
|
||||
# Bad units should raise an exception
|
||||
with pytest.raises(ValueError, match="unit"):
|
||||
integrator(op, ref_timesteps, power, timestep_units='🐨')
|
||||
with pytest.raises(ValueError, match="unit"):
|
||||
integrator(op, [(800.0, 'gorillas')], power)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue