crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Remove si_leqi in favor of SI_LEQI_Integrator

Browse source 
The depletion function openmc.deplete.si_leqi has been
removed in favor of the SI_LEQI_Integrator class. The same
depletion scheme can be obtained with the following commands:

   >>> leqi = openmc.deplete.SI_LEQI_Integrator(operator, dt, power)
   >>> leqi.integrate()

The expression can be onlined for compactness.

The si_celi_inner function has been removed completely now,
as the SI_CELI iteration is performed by directly calling
SI_CELI_Integrator.__call__ through the SI_LEQI_Integrator.
This is similar to how the LEQIIntegrator handles the initial steps.

Tests have been updated to use this class, and the class has been
added to the documentation. No pure-function integration
schemes exist anymore.
This commit is contained in:
Andrew Johnson 2019-07-29 10:29:24 -05:00
parent 791f7c78b7
commit 1f3e7fdac1
No known key found for this signature in database
GPG key ID: 253418E91B7F6FEB
5 changed files with 75 additions and 211 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -11,13 +11,6 @@ algorithms for depletion calculations, which are described in detail in Colin
Josey's thesis, `Development and analysis of high order neutron
transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
.. autosummary::
:toctree: generated
:nosignatures:
:template: myfunction.rst
integrator.si_leqi
.. autosummary::
:toctree: generated
:nosignatures:
@ -30,6 +23,7 @@ transport-depletion coupling algorithms <http://hdl.handle.net/1721.1/113721>`_.
integrator.EPC_RK4_Integrator
integrator.LEQIIntegrator
integrator.SI_CELI_Integrator
integrator.SI_LEQI_Integrator
Each of these functions expects a "transport operator" to be passed. An operator
specific to OpenMC is available using the following class:

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

@ -71,71 +71,3 @@ class SI_CELI_Integrator(SI_Integrator):
# 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):
""" The inner loop of SI-CE/LI CFQ4.
Parameters
----------
operator : Operator
The operator object to simulate on.
x : list of nuclide vector
Nuclide vector, beginning of time.
op_results : list of OperatorResult
Operator result at BOS.
p : float
Power of the reactor in [W]
i : int
Current iteration number.
i_res : int
Starting index, for restart calculation.
t : float
Time at start of step.
dt : float
Time step.
print_out : bool
Whether or not to print out time.
m : int, optional
Number of stages.
Returns
-------
list of nuclide vector (numpy.array)
Nuclide vector, end of time.
float
Next time
list of OperatorResult
Operator result at end of time.
"""
chain = operator.chain
# Deplete to end
proc_time, x_new = timed_deplete(
chain, x[0], op_results[0].rates, dt, print_out)
x.append(x_new)
for j in range(m + 1):
op_res = operator(x_new, p)
if j <= 1:
op_res_bar = copy.deepcopy(op_res)
else:
rates = 1/j * op_res.rates + (1 - 1/j) * op_res_bar.rates
k = 1/j * op_res.k + (1 - 1/j) * op_res_bar.k
op_res_bar = OperatorResult(k, rates)
rates = list(zip(op_results[0].rates, op_res_bar.rates))
time_1, x_new = timed_deplete(
chain, x[0], rates, dt, print_out, matrix_func=_celi_f1)
time_2, x_new = timed_deplete(
chain, x_new, rates, dt, print_out, matrix_func=_celi_f2)
proc_time += time_1 + time_2
# Create results, write to disk
op_results.append(op_res_bar)
Results.save(operator, x, op_results, [t, t+dt], p, i_res+i, proc_time)
# return updated time and vectors
return [x_new], t + dt, [op_res_bar]

197
openmc/deplete/integrator/si_leqi.py
View file

@ -6,15 +6,15 @@ from itertools import repeat
from uncertainties import ufloat
from .si_celi import si_celi_inner
from .abc import SI_Integrator
from .si_celi import SI_CELI_Integrator
from .leqi import _leqi_f1, _leqi_f2, _leqi_f3, _leqi_f4
from .cram import timed_deplete
from ..results import Results
from ..abc import OperatorResult
def si_leqi(operator, timesteps, power=None, power_density=None,
print_out=True, m=10):
class SI_LEQI_Integrator(SI_Integrator):
r"""Deplete using the SI-LE/QI CFQ4 algorithm.
Implements the Stochastic Implicit LE/QI Predictor-Corrector algorithm using
@ -22,138 +22,75 @@ def si_leqi(operator, timesteps, power=None, power_density=None,
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
def __call__(self, bos_conc, bos_rates, dt, power, i):
"""Perform the integration across one time step
Parameters
----------
bos_conc : list of numpy.ndarray
Initial concentrations for all nuclides in [atom] for
all depletable materials
bos_rates : list of openmc.deplete.ReactionRates
Reaction rates from operator for all depletable materials
dt : float
Time in [s] for the entire depletion interval
power : float
Power of the system [W]
i : int
Current depletion step index
Returns
-------
proc_time : float
Time spent in CRAM routines for all materials
conc_list : list of numpy.ndarray
Concentrations at each of the intermediate points with
the final concentration as the last element
op_results : list of openmc.deplete.OperatorResult
Eigenvalue and reaction rates from intermediate transport
simulation
"""
if i == 0:
if self._ires <= 1:
self._prev_rates = bos_rates
# Perform CELI for initial steps
return SI_CELI_Integrator.__call__(
self, bos_conc, bos_rates, dt, power, i)
prev_res = self.operator.prev_res[-2]
prevdt = self.timesteps[i] - prev_res.time[0]
self._prev_rates = prev_res.rates[0]
else:
power = [i*operator.heavy_metal for i in power_density]
prevdt = self.timesteps[i - 1]
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Perform remaining LE/QI
inputs = list(zip(self._prev_rates, bos_rates,
repeat(prevdt), repeat(dt)))
proc_time, inter_conc = timed_deplete(
self.chain, bos_conc, inputs, dt, matrix_func=_leqi_f1)
time1, eos_conc = timed_deplete(
self.chain, inter_conc, inputs, dt, matrix_func=_leqi_f2)
# 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)
proc_time += time1
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]]
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])
inputs = list(zip(self._prev_rates, bos_rates, res_bar.rates,
repeat(prevdt), repeat(dt)))
time1, inter_conc = timed_deplete(
self.chain, bos_conc, inputs, dt, matrix_func=_leqi_f3)
time2, inter_conc = timed_deplete(
self.chain, inter_conc, inputs, dt, matrix_func=_leqi_f4)
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]
chain = operator.chain
for i, (dt, p) in enumerate(zip(timesteps, power)):
# LE/QI needs the last step results to start
# Perform SI-CE/LI CFQ4 or restore results for the first step
if i == 0:
dt_l = dt
if i_res <= 1:
op_res_last = copy.deepcopy(op_results[0])
x, t, op_results = si_celi_inner(operator, x, op_results, p,
i, i_res, t, dt, print_out)
continue
else:
dt_l = t - operator.prev_res[-2].time[0]
op_res_last = operator.prev_res[-2]
op_res_last.rates = op_res_last.rates[0]
x = [operator.prev_res[-1].data[0]]
# Perform remaining LE/QI
inputs = list(zip(op_res_last.rates, op_results[0].rates,
repeat(dt_l), repeat(dt)))
proc_time, x_new = timed_deplete(
chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f1)
time_1, x_new = timed_deplete(
chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f2)
x.append(x_new)
proc_time += time_1
# Loop on inner
for j in range(m + 1):
op_res = operator(x_new, p)
if j <= 1:
op_res_bar = copy.deepcopy(op_res)
else:
rates = 1/j * op_res.rates + (1 - 1/j) * op_res_bar.rates
k = 1/j * op_res.k + (1 - 1/j) * op_res_bar.k
op_res_bar = OperatorResult(k, rates)
inputs = list(zip(op_res_last.rates, op_results[0].rates,
op_res_bar.rates, repeat(dt_l), repeat(dt)))
time_1, x_new = timed_deplete(
chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f3)
time_2, x_new = timed_deplete(
chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f4)
proc_time += time_1 + time_2
# Create results, write to disk
op_results.append(op_res_bar)
Results.save(
operator, x, op_results, [t, t+dt], p, i_res+i, proc_time)
# update results
x = [x_new]
op_res_last = copy.deepcopy(op_results[0])
op_results = [op_res_bar]
t += dt
dt_l = dt
# Create results for last point, write to disk
Results.save(
operator, x, op_results, [t, t], p, i_res+len(timesteps))
return proc_time, [eos_conc, inter_conc], [res_bar]

7
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, SI_CELI_Integrator
EPC_RK4_Integrator, CF4Integrator, SI_CELI_Integrator, SI_LEQI_Integrator
)
from tests import dummy_operator
@ -380,7 +380,8 @@ def test_restart_si_leqi(run_in_tmpdir):
# Perform simulation
dt = [0.75]
power = 1.0
openmc.deplete.si_leqi(op, dt, power, print_out=False)
nstages = 10
SI_LEQI_Integrator(op, dt, power, nstages).integrate()
# Load the files
prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")
@ -390,7 +391,7 @@ def test_restart_si_leqi(run_in_tmpdir):
op.output_dir = output_dir
# Perform restarts simulation
openmc.deplete.si_leqi(op, dt, power, print_out=False)
SI_LEQI_Integrator(op, dt, power, nstages).integrate()
# Load the files
res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5")

6
tests/unit_tests/test_deplete_si_leqi.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_LEQI_Integrator, ResultsList
from tests import dummy_operator
@ -18,10 +18,10 @@ def test_si_leqi(run_in_tmpdir):
# Perform simulation using the si_leqi algorithm
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.si_leqi(op, dt, power, print_out=False)
SI_LEQI_Integrator(op, dt, power, 10).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")