mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Demonstrate Integrator class with cecm
New abstract class openmc.deplete.integrate.Integrator and concrete CECMIntegrator for performing depletion analysis. Concrete classes only have to implement the __call__ method responsible for performing the time integration across a time interval. The abstract base class is responsible for iterating through all time steps, collecting and writing results to file, and executing intermediate transport solutions.
This commit is contained in:
Andrew Johnson
parent
e05dec4234
commit
cd444510cf
3 changed files with 194 additions and 91 deletions
|
|
@ -5,6 +5,7 @@ Integrator
|
|||
The integrator subcomponents.
|
||||
"""
|
||||
|
||||
from .abc import Integrator
|
||||
from .cf4 import *
|
||||
from .cecm import *
|
||||
from .celi import *
|
||||
|
|
|
|||
152
openmc/deplete/integrator/abc.py
Normal file
152
openmc/deplete/integrator/abc.py
Normal file
|
|
@ -0,0 +1,152 @@
|
|||
from copy import deepcopy
|
||||
from abc import ABC, abstractmethod
|
||||
from collections.abc import Iterable
|
||||
|
||||
from openmc.deplete import Results
|
||||
|
||||
|
||||
class Integrator(ABC):
|
||||
"""Abstract class for solving the time-integration for depletion
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, operator, timesteps, power=None, power_density=None):
|
||||
"""
|
||||
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.
|
||||
"""
|
||||
self.operator = operator
|
||||
self.chain = operator.chain
|
||||
if not isinstance(timesteps, Iterable):
|
||||
self.timesteps = [timesteps]
|
||||
else:
|
||||
self.timesteps = timesteps
|
||||
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]
|
||||
|
||||
if not isinstance(power, Iterable):
|
||||
# TODO Maybe use itertools.zip_longest?
|
||||
# Ensure that power is single value if that is the case
|
||||
power = [power] * len(self.timesteps)
|
||||
|
||||
self.power = power
|
||||
|
||||
@abstractmethod
|
||||
def __call__(self, conc, rates, dt, power):
|
||||
"""Perform the integration across one time step
|
||||
|
||||
Parameters
|
||||
----------
|
||||
conc : numpy.ndarray
|
||||
Initial concentrations for all nuclides in [atom]
|
||||
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]
|
||||
|
||||
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
|
||||
simulations
|
||||
"""
|
||||
|
||||
def __iter__(self):
|
||||
for dt, p in zip(self.timesteps, self.power):
|
||||
yield dt, p
|
||||
|
||||
def __len__(self):
|
||||
return len(self.timesteps)
|
||||
|
||||
def _get_bos_data(self, step_index, step_power, prev_conc):
|
||||
if step_index > 0 or self.operator.prev_res is None:
|
||||
x = deepcopy(prev_conc)
|
||||
res = self.operator(x, step_power)
|
||||
else:
|
||||
# Get previous concentration
|
||||
x = self.operator.prev_res[-1].data[0]
|
||||
|
||||
# Get reaction rates and keff
|
||||
res = self.operator.prev_res[-1]
|
||||
res.rates = res.rates[0]
|
||||
res.k = res.k[0]
|
||||
|
||||
# Scale rates by ratio of powers
|
||||
res.rates *= step_power / res.power[0]
|
||||
return x, res
|
||||
|
||||
def _get_start_data(self):
|
||||
if self.operator.prev_res is None:
|
||||
return 0.0, 0
|
||||
return self.operator.prev_res[-1].time[-1], len(self.operator.prev_res)
|
||||
|
||||
def integrate_all(self):
|
||||
"""Perform the entire depletion process across all steps"""
|
||||
with self.operator as conc:
|
||||
t, i_start = self._get_start_data()
|
||||
|
||||
for i, (dt, p) in enumerate(self):
|
||||
conc, res = self._get_bos_data(i, p, conc)
|
||||
proc_time, conc_list, res_list = self(conc, res.rates, dt, p)
|
||||
|
||||
# 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, i_start + i,
|
||||
proc_time)
|
||||
|
||||
t += dt
|
||||
|
||||
# Final simulation
|
||||
res_list = [self.operator(conc, p)]
|
||||
self._save_results(
|
||||
[conc], res_list, [t, t], p, i_start + len(self))
|
||||
|
||||
def _save_results(self, conc_list, results_list, time_list, power,
|
||||
index, proc_time=None):
|
||||
"""Save the results at the end of of one step
|
||||
|
||||
Abstracted to support the predictor's unique save location
|
||||
"""
|
||||
Results.save(
|
||||
self.operator, conc_list, results_list, time_list,
|
||||
power, index, proc_time)
|
||||
|
||||
@classmethod
|
||||
def integrate(cls, operator, timesteps, power=None, power_density=None):
|
||||
"""High-level interface for depleting with this integrator"""
|
||||
return cls(operator, timesteps, power, power_density).integrate_all()
|
||||
|
|
@ -1,13 +1,12 @@
|
|||
"""The CE/CM integrator."""
|
||||
|
||||
import copy
|
||||
from collections.abc import Iterable
|
||||
from textwrap import dedent
|
||||
|
||||
from .abc import Integrator
|
||||
from .cram import timed_deplete
|
||||
from ..results import Results
|
||||
|
||||
|
||||
def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
|
||||
class CECMIntegrator(Integrator):
|
||||
r"""Deplete using the CE/CM algorithm.
|
||||
|
||||
Implements the second order `CE/CM predictor-corrector algorithm
|
||||
|
|
@ -24,100 +23,51 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True):
|
|||
A_c &= A(y_m, t_n + h/2) \\
|
||||
y_{n+1} &= \text{expm}(A_c h) y_n
|
||||
\end{aligned}
|
||||
|
||||
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.
|
||||
|
||||
"""
|
||||
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]
|
||||
|
||||
if not isinstance(power, Iterable):
|
||||
power = [power]*len(timesteps)
|
||||
def __call__(self, conc, rates, dt, power):
|
||||
"""Integrate using CE/CM
|
||||
|
||||
# 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)
|
||||
Parameters
|
||||
----------
|
||||
conc : numpy.ndarray
|
||||
Initial concentrations for all nuclides in [atom]
|
||||
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]
|
||||
|
||||
chain = operator.chain
|
||||
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 transport simulations
|
||||
"""
|
||||
# deplete across first half of inteval
|
||||
time0, x_middle = timed_deplete(self.chain, conc, rates, dt / 2)
|
||||
res_middle = self.operator(x_middle, power)
|
||||
|
||||
for i, (dt, p) in enumerate(zip(timesteps, power)):
|
||||
# Get beginning-of-timestep concentrations and reaction rates
|
||||
# Avoid doing first transport run if already done in previous
|
||||
# calculation
|
||||
if i > 0 or operator.prev_res is None:
|
||||
x = [copy.deepcopy(vec)]
|
||||
op_results = [operator(x[0], p)]
|
||||
# deplete across entire interval with BOS concentrations,
|
||||
# MOS reaction rates
|
||||
time1, x_end = timed_deplete(self.chain, conc, res_middle.rates, dt)
|
||||
|
||||
else:
|
||||
# Get initial concentration
|
||||
x = [operator.prev_res[-1].data[0]]
|
||||
return time0 + time1, [x_middle, x_end], [res_middle]
|
||||
|
||||
# Get rates
|
||||
op_results = [operator.prev_res[-1]]
|
||||
op_results[0].rates = op_results[0].rates[0]
|
||||
|
||||
# Set first stage value of keff
|
||||
op_results[0].k = op_results[0].k[0]
|
||||
def cecm(operator, timesteps, power=None, power_density=None, print_out=False):
|
||||
# TODO Remove print_out since depletion timings are stored
|
||||
return CECMIntegrator(
|
||||
operator, timesteps, power, power_density).integrate_all()
|
||||
|
||||
# Scale reaction rates by ratio of powers
|
||||
power_res = operator.prev_res[-1].power
|
||||
ratio_power = p / power_res
|
||||
op_results[0].rates *= ratio_power[0]
|
||||
|
||||
# Deplete for first half of timestep
|
||||
proc_time, x_middle = timed_deplete(
|
||||
chain, x[0], op_results[0].rates, dt/2, print_out)
|
||||
|
||||
# Get middle-of-timestep reaction rates
|
||||
x.append(x_middle)
|
||||
op_results.append(operator(x_middle, p))
|
||||
|
||||
# Deplete for full timestep using beginning-of-step materials
|
||||
# and middle-of-timestep reaction rates
|
||||
pt_end, x_end = timed_deplete(
|
||||
chain, x[0], op_results[1].rates, dt, print_out)
|
||||
|
||||
# Create results, write to disk
|
||||
Results.save(
|
||||
operator, x, op_results, [t, t + dt], p, i_res + i,
|
||||
proc_time + pt_end)
|
||||
|
||||
# Advance time, update vector
|
||||
t += dt
|
||||
vec = copy.deepcopy(x_end)
|
||||
|
||||
# Perform one last simulation
|
||||
x = [copy.deepcopy(vec)]
|
||||
op_results = [operator(x[0], power[-1])]
|
||||
|
||||
# Create results, write to disk
|
||||
Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))
|
||||
try:
|
||||
cecm.__doc__ = (
|
||||
dedent(CECMIntegrator.__doc__) + dedent(Integrator.__init__.__doc__))
|
||||
except AttributeError:
|
||||
pass
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue