From ccfee55115a10144a04e63eac2e744b7411d1c82 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 19 Feb 2018 17:26:44 -0600 Subject: [PATCH] Make timesteps and power specified in integrator functions --- openmc/deplete/abc.py | 8 ------- openmc/deplete/integrator/cecm.py | 24 ++++++++++++++----- openmc/deplete/integrator/predictor.py | 26 +++++++++++++++------ openmc/deplete/openmc_wrapper.py | 23 +++++++++++------- tests/dummy_geometry.py | 4 +++- tests/regression_tests/test_deplete_full.py | 17 +++++++------- tests/unit_tests/test_deplete_cecm.py | 5 ++-- tests/unit_tests/test_deplete_predictor.py | 5 ++-- 8 files changed, 68 insertions(+), 44 deletions(-) diff --git a/openmc/deplete/abc.py b/openmc/deplete/abc.py index 6916066816..7fda93aa28 100644 --- a/openmc/deplete/abc.py +++ b/openmc/deplete/abc.py @@ -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 diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index 760e3c89d8..7d6c11190f 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -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)) diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index 0640783318..038c0778da 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -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)) diff --git a/openmc/deplete/openmc_wrapper.py b/openmc/deplete/openmc_wrapper.py index 56a248c30c..e5898997bc 100644 --- a/openmc/deplete/openmc_wrapper.py +++ b/openmc/deplete/openmc_wrapper.py @@ -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) diff --git a/tests/dummy_geometry.py b/tests/dummy_geometry.py index c013bb0054..d662616698 100644 --- a/tests/dummy_geometry.py +++ b/tests/dummy_geometry.py @@ -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. diff --git a/tests/regression_tests/test_deplete_full.py b/tests/regression_tests/test_deplete_full.py index 6033e3f75a..1ae2cee1fa 100644 --- a/tests/regression_tests/test_deplete_full.py +++ b/tests/regression_tests/test_deplete_full.py @@ -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' diff --git a/tests/unit_tests/test_deplete_cecm.py b/tests/unit_tests/test_deplete_cecm.py index 66c3ee156c..3daa7a0489 100644 --- a/tests/unit_tests/test_deplete_cecm.py +++ b/tests/unit_tests/test_deplete_cecm.py @@ -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") diff --git a/tests/unit_tests/test_deplete_predictor.py b/tests/unit_tests/test_deplete_predictor.py index f1133f87d2..be8497f5f7 100644 --- a/tests/unit_tests/test_deplete_predictor.py +++ b/tests/unit_tests/test_deplete_predictor.py @@ -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")