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

Make timesteps and power specified in integrator functions

Browse source
This commit is contained in:
Paul Romano 2018-02-19 17:26:44 -06:00
parent c19f396aa7
commit ccfee55115
8 changed files with 68 additions and 44 deletions
Download patch file Download diff file Expand all files Collapse all files

8
openmc/deplete/abc.py
View file

@ -18,8 +18,6 @@ class Settings(object):
Attributes
----------
dt_vec : numpy.array
Array of time steps to take.
output_dir : pathlib.Path
Path to output directory to save results.
chain_file : str
@ -29,10 +27,6 @@ class Settings(object):
Initial atom density to add for nuclides that are zero in initial
condition to ensure they exist in the decay chain. Only done for
nuclides with reaction rates. Defaults to 1.0e3.
power : float
Power of the reactor in [W]. 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.
"""
def __init__(self):
@ -40,9 +34,7 @@ class Settings(object):
self.chain_file = os.environ["OPENMC_DEPLETE_CHAIN"]
except KeyError:
self.chain_file = None
self.dt_vec = None
self.output_dir = '.'
self.power = None
self.dilute_initial = 1.0e3
@property

24
openmc/deplete/integrator/cecm.py
View file

@ -1,12 +1,13 @@
"""The CE/CM integrator."""
import copy
from collections.abc import Iterable
from .cram import deplete
from .save_results import save_results
def cecm(operator, print_out=True):
def cecm(operator, timesteps, power, print_out=True):
r"""The CE/CM integrator.
Implements the second order CE/CM Predictor-Corrector algorithm [ref]_.
@ -32,25 +33,36 @@ def cecm(operator, print_out=True):
----------
operator : openmc.deplete.Operator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]
power : float or iterable of float
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].
print_out : bool, optional
Whether or not to print out time.
"""
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
chain = operator.chain
t = 0.0
for i, dt in enumerate(operator.settings.dt_vec):
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep reaction rates
x = [copy.deepcopy(vec)]
results = [operator(x[0])]
results = [operator(x[0], p)]
# Deplete for first half of timestep
x_middle = deplete(chain, x[0], results[0], dt/2, print_out)
# Get middle-of-timestep reaction rates
x.append(x_middle)
results.append(operator(x_middle))
results.append(operator(x_middle, p))
# Deplete for full timestep using beginning-of-step materials
x_end = deplete(chain, x[0], results[1], dt, print_out)
@ -64,7 +76,7 @@ def cecm(operator, print_out=True):
# Perform one last simulation
x = [copy.deepcopy(vec)]
results = [operator(x[0])]
results = [operator(x[0], power[-1])]
# Create results, write to disk
save_results(operator, x, results, [t, t], len(operator.settings.dt_vec))
save_results(operator, x, results, [t, t], len(timesteps))

26
openmc/deplete/integrator/predictor.py
View file

@ -1,13 +1,14 @@
"""The Predictor algorithm."""
"""First-order predictor algorithm."""
import copy
from collections.abc import Iterable
from .cram import deplete
from .save_results import save_results
def predictor(operator, print_out=True):
r"""The basic predictor integrator.
def predictor(operator, timesteps, power, print_out=True):
r"""Deplete using a first-order predictor algorithm.
Implements the first-order predictor algorithm. This algorithm is
mathematically defined as:
@ -23,18 +24,29 @@ def predictor(operator, print_out=True):
----------
operator : openmc.deplete.Operator
The operator object to simulate on.
timesteps : iterable of float
Array of timesteps in units of [s]
power : float or iterable of float
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].
print_out : bool, optional
Whether or not to print out time.
"""
if not isinstance(power, Iterable):
power = [power]*len(timesteps)
# Generate initial conditions
with operator as vec:
chain = operator.chain
t = 0.0
for i, dt in enumerate(operator.settings.dt_vec):
for i, (dt, p) in enumerate(zip(timesteps, power)):
# Get beginning-of-timestep reaction rates
x = [copy.deepcopy(vec)]
results = [operator(x[0])]
results = [operator(x[0], p)]
# Create results, write to disk
save_results(operator, x, results, [t, t + dt], i)
@ -48,7 +60,7 @@ def predictor(operator, print_out=True):
# Perform one last simulation
x = [copy.deepcopy(vec)]
results = [operator(x[0])]
results = [operator(x[0], power[-1])]
# Create results, write to disk
save_results(operator, x, results, [t, t], len(operator.settings.dt_vec))
save_results(operator, x, results, [t, t], len(timesteps))

23
openmc/deplete/openmc_wrapper.py
View file

@ -45,8 +45,6 @@ class OpenMCSettings(Settings):
Attributes
----------
dt_vec : numpy.array
Array of time steps to in units of [s]
output_dir : pathlib.Path
Path to output directory to save results.
chain_file : str
@ -68,8 +66,8 @@ class OpenMCSettings(Settings):
"""
_depletion_attrs = {'dt_vec', '_output_dir', 'chain_file', 'dilute_initial',
'round_number', 'power'}
_depletion_attrs = {'_output_dir', 'chain_file', 'dilute_initial',
'round_number'}
def __init__(self):
super().__init__()
@ -152,13 +150,15 @@ class OpenMCOperator(Operator):
self.reaction_rates = ReactionRates(
self.local_mats, self._burnable_nucs, index_rx)
def __call__(self, vec, print_out=True):
def __call__(self, vec, power, print_out=True):
"""Runs a simulation.
Parameters
----------
vec : list of numpy.array
Total atoms to be used in function.
power : float
Power of the reactor in [W]
print_out : bool, optional
Whether or not to print out time.
@ -187,7 +187,7 @@ class OpenMCOperator(Operator):
time_openmc = time.time()
# Extract results
op_result = self._unpack_tallies_and_normalize()
op_result = self._unpack_tallies_and_normalize(power)
if comm.rank == 0:
time_unpack = time.time()
@ -424,7 +424,7 @@ class OpenMCOperator(Operator):
self._tally.scores = self.chain.reactions
self._tally.filters = [mat_filter]
def _unpack_tallies_and_normalize(self):
def _unpack_tallies_and_normalize(self, power):
"""Unpack tallies from OpenMC and return an operator result
This method uses OpenMC's C API bindings to determine the k-effective
@ -432,6 +432,11 @@ class OpenMCOperator(Operator):
normalized by the user-specified power, summing the product of the
fission reaction rate times the fission Q value for each material.
Parameters
----------
power : float
Power of the reactor in [W]
Returns
-------
openmc.deplete.OperatorResult
@ -509,10 +514,10 @@ class OpenMCOperator(Operator):
energy = comm.allreduce(energy)
# Determine power in eV/s
power = self.settings.power / JOULE_PER_EV
power /= JOULE_PER_EV
# Scale reaction rates to obtain units of reactions/sec
rates[:, :, :] *= power / energy
rates *= power / energy
return OperatorResult(k_combined, rates)

4
tests/dummy_geometry.py
View file

@ -21,13 +21,15 @@ class DummyGeometry(Operator):
def __init__(self, settings):
super().__init__(settings)
def __call__(self, vec, print_out=False):
def __call__(self, vec, power, print_out=False):
"""Evaluates F(y)
Parameters
----------
vec : list of numpy.array
Total atoms to be used in function.
power : float
Power in [W]
print_out : bool, optional, ignored
Whether or not to print out time.

17
tests/regression_tests/test_deplete_full.py
View file

@ -32,18 +32,10 @@ def test_full(run_in_tmpdir):
# Load geometry from example
geometry, lower_left, upper_right = generate_problem(n_rings, n_wedges)
# Create dt vector for 3 steps with 15 day timesteps
dt1 = 15.*24*60*60 # 15 days
dt2 = 1.5*30*24*60*60 # 1.5 months
N = floor(dt2/dt1)
dt = np.full(N, dt1)
# Depletion settings
settings = openmc.deplete.OpenMCSettings()
settings.chain_file = str(Path(__file__).parents[2] / 'chains' /
'chain_simple.xml')
settings.power = 2.337e15*4*JOULE_PER_EV*1e6 # MeV/second cm from CASMO
settings.dt_vec = dt
settings.round_number = True
# Add OpenMC-specific settings
@ -57,8 +49,15 @@ def test_full(run_in_tmpdir):
op = openmc.deplete.OpenMCOperator(geometry, settings)
# Power and timesteps
dt1 = 15.*24*60*60 # 15 days
dt2 = 1.5*30*24*60*60 # 1.5 months
N = floor(dt2/dt1)
dt = np.full(N, dt1)
power = 2.337e15*4*JOULE_PER_EV*1e6 # MeV/second cm from CASMO
# Perform simulation using the predictor algorithm
openmc.deplete.integrator.predictor(op)
openmc.deplete.integrator.predictor(op, dt, power)
# Get path to test and reference results
path_test = settings.output_dir / 'depletion_results.h5'

5
tests/unit_tests/test_deplete_cecm.py
View file

@ -15,13 +15,14 @@ def test_cecm(run_in_tmpdir):
"""Integral regression test of integrator algorithm using CE/CM."""
settings = openmc.deplete.Settings()
settings.dt_vec = [0.75, 0.75]
settings.output_dir = "test_integrator_regression"
op = dummy_geometry.DummyGeometry(settings)
# Perform simulation using the MCNPX/MCNP6 algorithm
openmc.deplete.cecm(op, print_out=False)
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.cecm(op, dt, power, print_out=False)
# Load the files
res = results.read_results(settings.output_dir / "depletion_results.h5")

5
tests/unit_tests/test_deplete_predictor.py
View file

@ -15,13 +15,14 @@ def test_predictor(run_in_tmpdir):
"""Integral regression test of integrator algorithm using predictor/corrector"""
settings = openmc.deplete.Settings()
settings.dt_vec = [0.75, 0.75]
settings.output_dir = "test_integrator_regression"
op = dummy_geometry.DummyGeometry(settings)
# Perform simulation using the predictor algorithm
openmc.deplete.predictor(op, print_out=False)
dt = [0.75, 0.75]
power = 1.0
openmc.deplete.predictor(op, dt, power, print_out=False)
# Load the files
res = results.read_results(settings.output_dir / "depletion_results.h5")