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

Add SI_Integrator, SI_CELI_Integrator classes; remove si_celi

Browse source 
The si_celi depletion function has been removed. The
SI-CELI method can be implemented with
    >>> from openmc.deplete import SI_CELI_Integrator
    >>> SI_CELI_Integrator(op, dt, power, stages).integrate()

The stages parameter is optional, and defaults to 10 to
be consistent with the previous default for si_celi.

The SI_CELI_Integrator inherits from the new abstract base class
SI_Integrator. There are a few differences in how the SI-based
methods perform the integration.

- The initial, t=0.0, i=0, no restart transport solution is scaled
  by the number of stages.
- The SI methods also do not perform a new transport solution at
  each successive iteration [t>0, i>0]. Instead, the reaction
  rates are pulled from the last stage of the previous step.
- There is no final transport solution once all the depletion
  steps have been taken. The final reaction rates are saved as
  those from the last stage of the last depletion step.

The si_celi function has been removed from documentation and testing,
and replaced with the SI_CELI_Integrator.
This commit is contained in:
Andrew Johnson 2019-07-26 12:05:06 -05:00
parent 98957fb473
commit d28546bcd6
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
5 changed files with 143 additions and 86 deletions
Download patch file Download diff file Expand all files Collapse all files

2
docs/source/pythonapi/deplete.rst
View file

@ -16,7 +16,6 @@ transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
:nosignatures:
:template: myfunction.rst
integrator.si_celi
integrator.si_leqi
.. autosummary::
@ -30,6 +29,7 @@ transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
integrator.CELIIntegrator
integrator.EPC_RK4_Integrator
integrator.LEQIIntegrator
integrator.SI_CELI_Integrator
Each of these functions expects a "transport operator" to be passed. An operator
specific to OpenMC is available using the following class:

84
openmc/deplete/integrator/abc.py
View file

@ -149,3 +149,87 @@ class Integrator(ABC):
Results.save(
self.operator, conc_list, results_list, time_list,
power, index, proc_time)
class SI_Integrator(Integrator):
"""Abstract for the Stochastic Implicit Euler integrators
Does not provide a ``__call__`` method, but scales and resets
the number of particles used in initial transport calculation
"""
def __init__(self, operator, timesteps, power=None, power_density=None,
n_steps=10):
"""
Parameters
----------
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.
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
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.
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.
n_stages : int, optional
Number of stages. Default : 10
"""
def __init__(self, operator, timesteps, power=None, power_density=None,
n_stages=10):
super().__init__(operator, timesteps, power, power_density)
self.n_stages = n_stages
def _get_bos_data_from_openmc(self, step_index, step_power, bos_conc):
reset_particles = False
if step_index == 0 and hasattr(self.operator, "settings"):
reset_particles = True
self.operator.settings.particles *= self.n_stages
inherited = super()._get_bos_data_from_openmc(
step_index, step_power, bos_conc)
if reset_particles:
self.operator.settings.particles //= self.n_stages
return inherited
def integrate(self):
"""Perform the entire depletion process across all steps"""
with self.operator as conc:
t, self._ires = self._get_start_data()
for i, (dt, p) in enumerate(self):
if i == 0:
if self.operator.prev_res is None:
conc, res = self._get_bos_data_from_openmc(i, p, conc)
else:
conc, res = self._get_bos_data_from_restart(i, p, conc)
else:
# Pull rates, k from previous iteration w/o
# re-running transport
res = res_list[-1] # defined in previous i iteration
proc_time, conc_list, res_list = self(conc, res.rates, dt, p, i)
# Insert BOS concentration, transport results
conc_list.insert(0, conc)
res_list.insert(0, res)
# Remove actual EOS concentration for next step
conc = conc_list.pop()
self._save_results(
conc_list, res_list, [t, t + dt], p, self._ires + i,
proc_time)
t += dt
# No final simulation for SIE, use last iteration results
self._save_results(
[conc], [res_list[-1]], [t, t], p, self._ires + len(self))

131
openmc/deplete/integrator/si_celi.py
View file

@ -1,103 +1,76 @@
"""The SI-CE/LI CFQ4 integrator."""
import copy
from collections.abc import Iterable
from uncertainties import ufloat
from .abc import SI_Integrator
from .cram import timed_deplete
from ..results import Results
from ..abc import OperatorResult
from .celi import _celi_f1, _celi_f2
from ..abc import OperatorResult
from ..results import Results
def si_celi(operator, timesteps, power=None, power_density=None,
print_out=True, m=10):
r"""Deplete using the SI-CE/LI CFQ4 algorithm.
class SI_CELI_Integrator(SI_Integrator):
r"""Deplete using the si-ce/li cfq4 algorithm.
Implements the Stochastic Implicit CE/LI Predictor-Corrector algorithm using
Implements the stochastic implicit ce/li predictor-corrector algorithm using
the `fourth order commutator-free integrator <https://doi.org/10.1137/05063042>`_.
Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis
detailed algorithm can be found in section 3.2 in `colin josey's thesis
<http://hdl.handle.net/1721.1/113721>`_.
Parameters
----------
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.
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 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.
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.
print_out : bool, optional
Whether or not to print out time.
m : int, optional
Number of stages.
"""
if power is None:
if power_density is None:
raise ValueError(
"Neither power nor power density was specified.")
if not isinstance(power_density, Iterable):
power = power_density*operator.heavy_metal
else:
power = [i*operator.heavy_metal for i in power_density]
def __call__(self, bos_conc, bos_rates, dt, power, _i):
"""Perform the integration across one time step
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
Parameters
----------
bos_conc : numpy.ndarray
Initial bos_concentrations for all nuclides in [atom]
bos_rates : openmc.deplete.ReactionRates
Reaction rates from operator
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system [W]
_i : int
Current depletion step index. Unused
# Generate initial conditions
with operator as vec:
# Initialize time and starting index
if operator.prev_res is None:
t = 0.0
i_res = 0
else:
t = operator.prev_res[-1].time[-1]
i_res = len(operator.prev_res)
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials
bos_conc_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final bos_concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulations
"""
proc_time, eos_conc = timed_deplete(
self.chain, bos_conc, bos_rates, dt)
inter_conc = copy.deepcopy(eos_conc)
# Get the concentrations and reaction rates for the first
# beginning-of-timestep (BOS). Compute with m (stage number) times as
# many neutrons as later simulations for statistics reasons if no
# previous calculation results present
if operator.prev_res is None:
x = [copy.deepcopy(vec)]
if hasattr(operator, "settings"):
operator.settings.particles *= m
op_results = [operator(x[0], power[0])]
if hasattr(operator, "settings"):
operator.settings.particles //= m
else:
# Get initial concentration
x = [operator.prev_res[-1].data[0]]
# Begin iteration
for j in range(self.n_stages + 1):
inter_res = self.operator(inter_conc, power)
# Get rates
op_results = [operator.prev_res[-1]]
op_results[0].rates = op_results[0].rates[0]
if j <= 1:
res_bar = copy.deepcopy(inter_res)
else:
rates = 1/j * inter_res.rates + (1 - 1 / j) * res_bar.rates
k = 1/j * inter_res.k + (1 - 1 / j) * res_bar.k
res_bar = OperatorResult(k, rates)
# Set first stage value of keff
op_results[0].k = ufloat(*op_results[0].k[0])
list_rates = list(zip(bos_rates, res_bar.rates))
time1, inter_conc = timed_deplete(
self.chain, bos_conc, list_rates, dt, matrix_func=_celi_f1)
time2, inter_conc = timed_deplete(
self.chain, inter_conc, list_rates, dt, matrix_func=_celi_f2)
proc_time += time1 + time2
# Scale reaction rates by ratio of powers
power_res = operator.prev_res[-1].power
ratio_power = power[0] / power_res
op_results[0].rates *= ratio_power[0]
for i, (dt, p) in enumerate(zip(timesteps, power)):
x, t, op_results = si_celi_inner(operator, x, op_results, p,
i, i_res, t, dt, print_out, m)
# Create results for last point, write to disk
Results.save(
operator, x, op_results, [t, t], p, i_res + len(timesteps))
# end iteration
return proc_time, [eos_conc, inter_conc], [res_bar]
def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out, m=10):

6
tests/unit_tests/test_deplete_restart.py
View file

@ -8,7 +8,7 @@ from pytest import approx
import openmc.deplete
from openmc.deplete import (
CECMIntegrator, PredictorIntegrator, CELIIntegrator, LEQIIntegrator,
EPC_RK4_Integrator, CF4Integrator
EPC_RK4_Integrator, CF4Integrator, SI_CELI_Integrator
)
from tests import dummy_operator
@ -341,7 +341,7 @@ def test_restart_si_celi(run_in_tmpdir):
# Perform simulation
dt = [0.75]
power = 1.0
openmc.deplete.si_celi(op, dt, power, print_out=False)
SI_CELI_Integrator(op, dt, power).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
@ -351,7 +351,7 @@ def test_restart_si_celi(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation
openmc.deplete.si_celi(op, dt, power, print_out=False)
SI_CELI_Integrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")

6
tests/unit_tests/test_deplete_si_celi.py
View file

@ -4,7 +4,7 @@ These tests integrate a simple test problem described in dummy_geometry.py.
"""
from pytest import approx
import openmc.deplete
from openmc.deplete import SI_CELI_Integrator, ResultsList
from tests import dummy_operator
@ -18,10 +18,10 @@ def test_si_celi(run_in_tmpdir):
# Perform simulation using the si_celi algorithm
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.si_celi(op, dt, power, print_out=False)
SI_CELI_Integrator(op, dt, power).integrate()
# Load the files
res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
res = ResultsList(op.output_dir / "depletion_results.h5")
_, y1 = res.get_atoms("1", "1")
_, y2 = res.get_atoms("1", "2")