From 9c8eddddd75f827d11e63d513184946425ff4141 Mon Sep 17 00:00:00 2001 From: liangjg Date: Mon, 7 Jan 2019 20:42:10 -0500 Subject: [PATCH 01/24] pass reaction rates instead of op_results to depletion integrator --- openmc/deplete/integrator/cecm.py | 6 +++--- openmc/deplete/integrator/cram.py | 17 +++++------------ openmc/deplete/integrator/predictor.py | 4 ++-- 3 files changed, 10 insertions(+), 17 deletions(-) diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index ac519cf87..4554a3f22 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -95,10 +95,10 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True): # Scale reaction rates by ratio of powers power_res = operator.prev_res[-1].power ratio_power = p / power_res - op_results[0].rates[0] *= ratio_power[0] + op_results[0].rates *= ratio_power[0] # Deplete for first half of timestep - x_middle = deplete(chain, x[0], op_results[0], dt/2, print_out) + x_middle = deplete(chain, x[0], op_results[0].rates, dt/2, print_out) # Get middle-of-timestep reaction rates x.append(x_middle) @@ -106,7 +106,7 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True): # Deplete for full timestep using beginning-of-step materials # and middle-of-timestep reaction rates - x_end = deplete(chain, x[0], op_results[1], dt, print_out) + x_end = 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) diff --git a/openmc/deplete/integrator/cram.py b/openmc/deplete/integrator/cram.py index 85954603a..ea2d55b7d 100644 --- a/openmc/deplete/integrator/cram.py +++ b/openmc/deplete/integrator/cram.py @@ -14,7 +14,7 @@ import scipy.sparse.linalg as sla from .. import comm -def deplete(chain, x, op_result, dt, print_out): +def deplete(chain, x, rates, dt, print_out): """Deplete materials using given reaction rates for a specified time Parameters @@ -23,8 +23,8 @@ def deplete(chain, x, op_result, dt, print_out): Depletion chain x : list of numpy.ndarray Atom number vectors for each material - op_result : openmc.deplete.OperatorResult - Result of applying transport operator (contains reaction rates) + rates : openmc.deplete.ReactionRates + Reaction rates (from transport operator) dt : float Time in [s] to deplete for print_out : bool @@ -38,16 +38,9 @@ def deplete(chain, x, op_result, dt, print_out): """ t_start = time.time() - # Set up iterators - n_mats = len(x) - chains = repeat(chain, n_mats) - vecs = (x[i] for i in range(n_mats)) - rates = (op_result.rates[i, :, :] for i in range(n_mats)) - dts = repeat(dt, n_mats) - # Use multiprocessing pool to distribute work with Pool() as pool: - iters = zip(chains, vecs, rates, dts) + iters = zip(repeat(chain), x, rates, repeat(dt)) x_result = list(pool.starmap(_cram_wrapper, iters)) t_end = time.time() @@ -63,7 +56,7 @@ def _cram_wrapper(chain, n0, rates, dt): Parameters ---------- - chain : DepletionChain + chain : openmc.deplete.Chain Depletion chain used to construct the burnup matrix n0 : numpy.array Vector to operate a matrix exponent on. diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index 969144f68..f480f55fb 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -90,10 +90,10 @@ def predictor(operator, timesteps, power=None, power_density=None, # Scale reaction rates by ratio of powers power_res = operator.prev_res[-1].power ratio_power = p / power_res - op_results[0].rates[0] *= ratio_power[0] + op_results[0].rates *= ratio_power[0] # Deplete for full timestep - x_end = deplete(chain, x[0], op_results[0], dt, print_out) + x_end = deplete(chain, x[0], op_results[0].rates, dt, print_out) # Advance time, update vector t += dt From e36d1386587624aeafca828c74a2abfbe3fa6fb2 Mon Sep 17 00:00:00 2001 From: liangjg Date: Mon, 7 Jan 2019 20:52:34 -0500 Subject: [PATCH 02/24] allow to pass matrix function for high order methods --- openmc/deplete/integrator/cram.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/openmc/deplete/integrator/cram.py b/openmc/deplete/integrator/cram.py index ea2d55b7d..8a8940df8 100644 --- a/openmc/deplete/integrator/cram.py +++ b/openmc/deplete/integrator/cram.py @@ -14,7 +14,7 @@ import scipy.sparse.linalg as sla from .. import comm -def deplete(chain, x, rates, dt, print_out): +def deplete(chain, x, rates, dt, print_out=True, matrix_func=None): """Deplete materials using given reaction rates for a specified time Parameters @@ -27,8 +27,10 @@ def deplete(chain, x, rates, dt, print_out): Reaction rates (from transport operator) dt : float Time in [s] to deplete for - print_out : bool + print_out : bool, optional Whether to show elapsed time + maxtrix_func : function, optional + Function to form the depletion matrix Returns ------- @@ -40,7 +42,7 @@ def deplete(chain, x, rates, dt, print_out): # Use multiprocessing pool to distribute work with Pool() as pool: - iters = zip(repeat(chain), x, rates, repeat(dt)) + iters = zip(repeat(chain), x, rates, repeat(dt), repeat(matrix_func)) x_result = list(pool.starmap(_cram_wrapper, iters)) t_end = time.time() @@ -51,7 +53,7 @@ def deplete(chain, x, rates, dt, print_out): return x_result -def _cram_wrapper(chain, n0, rates, dt): +def _cram_wrapper(chain, n0, rates, dt, matrix_func=None): """Wraps depletion matrix creation / CRAM solve for multiprocess execution Parameters @@ -64,13 +66,19 @@ def _cram_wrapper(chain, n0, rates, dt): 2D array indexed by nuclide then by cell. dt : float Time to integrate to. + maxtrix_func : function, optional + Function to form the depletion matrix Returns ------- numpy.array Results of the matrix exponent. """ - A = chain.form_matrix(rates) + + if matrix_func is None: + A = chain.form_matrix(rates) + else: + A = matrix_func(chain, rates) return CRAM48(A, n0, dt) From 43881c94540aec323418d3ac12d4d978cae24917 Mon Sep 17 00:00:00 2001 From: liangjg Date: Mon, 7 Jan 2019 20:54:09 -0500 Subject: [PATCH 03/24] move cf4 implementation cf4.py from Colin's repo --- openmc/deplete/integrator/cf4.py | 169 +++++++++++++++++++++++++++++++ 1 file changed, 169 insertions(+) create mode 100644 openmc/deplete/integrator/cf4.py diff --git a/openmc/deplete/integrator/cf4.py b/openmc/deplete/integrator/cf4.py new file mode 100644 index 000000000..b099792fa --- /dev/null +++ b/openmc/deplete/integrator/cf4.py @@ -0,0 +1,169 @@ +""" The CF4 integrator. + +Implements the CF4 algorithm. +""" + +import copy +import os +import time + +from mpi4py import MPI + +from .cram import CRAM48 +from .save_results import save_results + +def cf4(operator, print_out=True): + """ Performs integration of an operator using the CF4 algorithm. + + Parameters + ---------- + operator : Operator + The operator object to simulate on. + print_out : bool, optional + Whether or not to print out time. + """ + + # Save current directory + dir_home = os.getcwd() + + # Move to folder + os.makedirs(operator.settings.output_dir, exist_ok=True) + os.chdir(operator.settings.output_dir) + + # Generate initial conditions + vec = operator.initial_condition() + + n_mats = len(vec) + + t = 0.0 + + for i, dt in enumerate(operator.settings.dt_vec): + # Create vectors + x = [copy.deepcopy(vec)] + seeds = [] + eigvls = [] + rates_array = [] + + eigvl, rates, seed = operator.eval(x[0]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrix + f0 = operator.form_matrix(rates_array[0], mat) + + x_new = CRAM48(1/2 * f0, x[0][mat], dt) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + x.append(x_result) + + eigvl, rates, seed = operator.eval(x[1]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrix + f1 = operator.form_matrix(rates_array[1], mat) + + x_new = CRAM48(1/2 * f1, x[0][mat], dt) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + x.append(x_result) + + eigvl, rates, seed = operator.eval(x[2]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrix + f0 = operator.form_matrix(rates_array[0], mat) + f2 = operator.form_matrix(rates_array[2], mat) + + x_new = CRAM48(-1/2 * f0 + f2, x[1][mat], dt) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + x.append(x_result) + + eigvl, rates, seed = operator.eval(x[3]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrix + f0 = operator.form_matrix(rates_array[0], mat) + f1 = operator.form_matrix(rates_array[1], mat) + f2 = operator.form_matrix(rates_array[2], mat) + f3 = operator.form_matrix(rates_array[3], mat) + + x_new = copy.deepcopy(x[0][mat]) + x_new = CRAM48(1/4*f0 + 1/6*f1 + 1/6*f2 - 1/12*f3, x_new, dt) + x_new = CRAM48(-1/12*f0 + 1/6*f1 + 1/6*f2 + 1/4*f3, x_new, dt) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i) + + t += dt + vec = copy.deepcopy(x_result) + + # Perform one last simulation + x = [copy.deepcopy(vec)] + seeds = [] + eigvls = [] + rates_array = [] + eigvl, rates, seed = operator.eval(x[0]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(rates) + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t], + len(operator.settings.dt_vec)) + + # Return to origin + os.chdir(dir_home) From a4d59f6150be2eec1264224394c5616136271ec4 Mon Sep 17 00:00:00 2001 From: liangjg Date: Mon, 7 Jan 2019 22:25:52 -0500 Subject: [PATCH 04/24] added cf4 integrator --- openmc/deplete/integrator/__init__.py | 1 + openmc/deplete/integrator/cf4.py | 291 +++++++++++++------------- 2 files changed, 144 insertions(+), 148 deletions(-) diff --git a/openmc/deplete/integrator/__init__.py b/openmc/deplete/integrator/__init__.py index cf8caffdf..18accdcc8 100644 --- a/openmc/deplete/integrator/__init__.py +++ b/openmc/deplete/integrator/__init__.py @@ -5,6 +5,7 @@ Integrator The integrator subcomponents. """ +from .cf4 import * from .cecm import * from .cram import * from .predictor import * diff --git a/openmc/deplete/integrator/cf4.py b/openmc/deplete/integrator/cf4.py index b099792fa..6fdd8bd42 100644 --- a/openmc/deplete/integrator/cf4.py +++ b/openmc/deplete/integrator/cf4.py @@ -1,169 +1,164 @@ -""" The CF4 integrator. - -Implements the CF4 algorithm. -""" +"""The CF4 integrator.""" import copy -import os -import time +from collections.abc import Iterable -from mpi4py import MPI +from .cram import deplete +from ..results import Results -from .cram import CRAM48 -from .save_results import save_results -def cf4(operator, print_out=True): - """ Performs integration of an operator using the CF4 algorithm. +# Functions to form the special matrix for depletion +def mf1(chain, rates): + return 1/2 * chain.form_matrix(rates) + +def mf3(chain, rates): + return -1/2 * chain.form_matrix(rates[0]) + chain.form_matrix(rates[1]) + +def mf4(chain, rates): + return 1/4 * chain.form_matrix(rates[0]) + \ + 1/6 * chain.form_matrix(rates[1]) + \ + 1/6 * chain.form_matrix(rates[2]) + \ + -1/12 * chain.form_matrix(rates[3]) + +def mf5(chain, rates): + return -1/12 * chain.form_matrix(rates[0]) + \ + 1/6 * chain.form_matrix(rates[1]) + \ + 1/6 * chain.form_matrix(rates[2]) + \ + 1/4 * chain.form_matrix(rates[3]) + +def cf4(operator, timesteps, power=None, power_density=None, print_out=True): + r"""Deplete using the CF4 algorithm. + + Implements the fourth order commutator-free Lie algorithm. + This algorithm is mathematically defined as: + + .. math:: + F_1 = h A(y_0) + + y_1 &= \text{expm}(1/2 F_1) y_0 + + F_2 &= h A(y_1) + + y_2 &= \text{expm}(1/2 F_2) y_0 + + F_3 &= h A(y_2) + + y_3 &= \text{expm}(-1/2 F_1 + F_3) y_1 + + F_4 &= h A(y_3) + + y_4 &= \text{expm}( 1/4 F_1 + 1/6 F_2 + 1/6 F_3 - 1/12 F_4) + \text{expm}(-1/12 F_1 + 1/6 F_2 + 1/6 F_3 + 1/4 F_4) y_0 Parameters ---------- - operator : Operator + 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] - # Save current directory - dir_home = os.getcwd() - - # Move to folder - os.makedirs(operator.settings.output_dir, exist_ok=True) - os.chdir(operator.settings.output_dir) + if not isinstance(power, Iterable): + power = [power]*len(timesteps) # Generate initial conditions - vec = operator.initial_condition() + with operator as vec: + chain = operator.chain - n_mats = len(vec) + # Initialize time + if operator.prev_res is None: + t = 0.0 + else: + t = operator.prev_res[-1].time[-1] - t = 0.0 + # Initialize starting index for saving results + if operator.prev_res is None: + i_res = 0 + else: + i_res = len(operator.prev_res) - for i, dt in enumerate(operator.settings.dt_vec): - # Create vectors + 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)] + + else: + # Get initial concentration + x = [operator.prev_res[-1].data[0]] + + # 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] + + # 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] + + # Step 1: deplete with matrix 1/2*A(y0) + x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out, + matrix_func=mf1) + x.append(x_new) + op_results.append(operator(x[1], p)) + + # Step 2: deplete with matrix 1/2*A(y1) + x_new = deplete(chain, x[0], op_results[1].rates, dt, print_out, + matrix_func=mf1) + x.append(x_new) + op_results.append(operator(x[2], p)) + + # Step 3: deplete with matrix -1/2*A(y0)+A(y2) + rates = list(zip(op_results[0].rates, op_results[2].rates)) + x_new = deplete(chain, x[1], rates, dt, print_out, + matrix_func=mf3) + x.append(x_new) + op_results.append(operator(x[3], p)) + + # Step 4: deplete with two matrix exponentials + rates = list(zip(op_results[0].rates, op_results[1].rates, + op_results[2].rates, op_results[3].rates)) + x_end = deplete(chain, x[0], rates, dt, print_out, + matrix_func=mf4) + x_end = deplete(chain, x_end, rates, dt, print_out, + matrix_func=mf5) + + # Create results, write to disk + Results.save(operator, x, op_results, [t, t + dt], p, i_res + i) + + # Advance time, update vector + t += dt + vec = copy.deepcopy(x_end) + + # Perform one last simulation x = [copy.deepcopy(vec)] - seeds = [] - eigvls = [] - rates_array = [] - - eigvl, rates, seed = operator.eval(x[0]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrix - f0 = operator.form_matrix(rates_array[0], mat) - - x_new = CRAM48(1/2 * f0, x[0][mat], dt) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - x.append(x_result) - - eigvl, rates, seed = operator.eval(x[1]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrix - f1 = operator.form_matrix(rates_array[1], mat) - - x_new = CRAM48(1/2 * f1, x[0][mat], dt) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - x.append(x_result) - - eigvl, rates, seed = operator.eval(x[2]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrix - f0 = operator.form_matrix(rates_array[0], mat) - f2 = operator.form_matrix(rates_array[2], mat) - - x_new = CRAM48(-1/2 * f0 + f2, x[1][mat], dt) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - x.append(x_result) - - eigvl, rates, seed = operator.eval(x[3]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrix - f0 = operator.form_matrix(rates_array[0], mat) - f1 = operator.form_matrix(rates_array[1], mat) - f2 = operator.form_matrix(rates_array[2], mat) - f3 = operator.form_matrix(rates_array[3], mat) - - x_new = copy.deepcopy(x[0][mat]) - x_new = CRAM48(1/4*f0 + 1/6*f1 + 1/6*f2 - 1/12*f3, x_new, dt) - x_new = CRAM48(-1/12*f0 + 1/6*f1 + 1/6*f2 + 1/4*f3, x_new, dt) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) + op_results = [operator(x[0], power[-1])] # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i) - - t += dt - vec = copy.deepcopy(x_result) - - # Perform one last simulation - x = [copy.deepcopy(vec)] - seeds = [] - eigvls = [] - rates_array = [] - eigvl, rates, seed = operator.eval(x[0]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(rates) - - # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t], - len(operator.settings.dt_vec)) - - # Return to origin - os.chdir(dir_home) + Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps)) From d3009d910146b62a5fba3adb34a6cef2046c1aff Mon Sep 17 00:00:00 2001 From: liangjg Date: Tue, 8 Jan 2019 11:41:45 -0500 Subject: [PATCH 05/24] moved more algorithms from Colin: epc-rk4, si-celi, si-leqi --- openmc/deplete/integrator/celi_cfq4_imp.py | 199 +++++++++++++++++++++ openmc/deplete/integrator/epc_rk4.py | 166 +++++++++++++++++ openmc/deplete/integrator/leqi_cfq4_imp.py | 191 ++++++++++++++++++++ 3 files changed, 556 insertions(+) create mode 100644 openmc/deplete/integrator/celi_cfq4_imp.py create mode 100644 openmc/deplete/integrator/epc_rk4.py create mode 100644 openmc/deplete/integrator/leqi_cfq4_imp.py diff --git a/openmc/deplete/integrator/celi_cfq4_imp.py b/openmc/deplete/integrator/celi_cfq4_imp.py new file mode 100644 index 000000000..a17f2fda9 --- /dev/null +++ b/openmc/deplete/integrator/celi_cfq4_imp.py @@ -0,0 +1,199 @@ +""" The CE/LI CFQ4 integrator. + +Implements the CE/LI Predictor-Corrector algorithm using commutator free +high order integrators. + +This algorithm is mathematically defined as: + +.. math: + y' = A(y, t) y(t) + A_p = A(y_n, t_n) + y_p = expm(A_p h) y_n + A_c = A(y_p, t_n) + A(t) = t/dt * A_c + (dt - t)/dt * A_p + +Here, A(t) is integrated using the fourth order algorithm described below. + +From +---- + Thalhammer, Mechthild. "A fourth-order commutator-free exponential + integrator for nonautonomous differential equations." SIAM journal on + numerical analysis 44.2 (2006): 851-864. +""" + +import copy +import os +import time + +from mpi4py import MPI + +from .cram import CRAM48 +from .save_results import save_results + +def celi_cfq4_imp(operator, print_out=True): + """ Performs integration of an operator using the CE/LI CFQ4 algorithm. + + Parameters + ---------- + operator : Operator + The operator object to simulate on. + print_out : bool, optional + Whether or not to print out time. + """ + + m = 10 + + # Save current directory + dir_home = os.getcwd() + + # Move to folder + os.makedirs(operator.settings.output_dir, exist_ok=True) + os.chdir(operator.settings.output_dir) + + # Generate initial conditions + vec = operator.initial_condition() + + # Compute initial rates + x = copy.deepcopy(vec) + + operator.settings.particles *= m + eigvl_bos, rates, seed = operator.eval(vec) + operator.settings.particles = int(operator.settings.particles / m) + + rates_bos = copy.deepcopy(rates) + + t = 0.0 + + for i, dt in enumerate(operator.settings.dt_vec): + vec, t, rates_eos, eigvl_bos = celi_cfq4_imp_inner(operator, vec, rates_bos, eigvl_bos, i, t, dt, print_out) + rates_bos = copy.deepcopy(rates_eos) + + # Perform one last simulation + x = [copy.deepcopy(vec)] + seeds = [0] + eigvls = [eigvl_bos] + rates_array = [copy.deepcopy(rates_bos)] + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t], + len(operator.settings.dt_vec)) + + # Return to origin + os.chdir(dir_home) + +def celi_cfq4_imp_inner(operator, vec, rates_bos, eigvl_bos, i, t, dt, print_out): + """ The inner loop of CE/LI CFQ4. + + Parameters + ---------- + operator : Operator + The operator object to simulate on. + vec : list of numpy.array + Nuclide vector, beginning of time. + i : Int + Current iteration number. + t : Float + Time at start of step. + dt : Float + Time step. + print_out : bool + Whether or not to print out time. + + Returns + ------- + x_result : list of numpy.array + Nuclide vector, end of time. + Float + Next time + ReactionRates + Reaction rates from beginning of step. + """ + + m = 10 + + n_mats = len(vec) + + # Create vectors + x = [copy.deepcopy(vec)] + seeds = [] + eigvls = [] + rates_array = [] + + eigvls.append(eigvl_bos) + seeds.append(0) + rates_array.append(copy.deepcopy(rates_bos)) + + # Deplete to end + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrix + f = operator.form_matrix(rates_array[0], mat) + + x_new = CRAM48(f, x[0][mat], dt) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + x.append(copy.deepcopy(x_result)) + + eigvl_bar = 0.0 + rates_bar = [] + + for j in range(0, m + 1): + eigvl, rates, seed = operator.eval(x_result) + + if j <= 1: + rates_bar = copy.deepcopy(rates) + eigvl_bar = eigvl + else: + rates_bar.rates = 1/j * rates.rates + (1 - 1/j) * rates_bar.rates + eigvl_bar = 1/j * eigvl + (1 - 1/j) * eigvl_bar + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrices + f1 = dt * operator.form_matrix(rates_array[0], mat) + f2 = dt * operator.form_matrix(rates_bar, mat) + + # Perform commutator-free integral + x_new = copy.deepcopy(x[0][mat]) + + # Compute linearly interpolated f at points + # A{1,2} = f(1/2 -/+ sqrt(3)/6) + # Then + # a{1,2} = 1/4 +/- sqrt(3)/6 + # m1 = a2 * A1 + a1 * A2 + # m2 = a1 * A1 + a2 * A2 + m1 = 1/12 * (f1 + 5 * f2) + m2 = 1/12 * (5 * f1 + f2) + + x_new = CRAM48(m2, x_new, 1.0) + x_new = CRAM48(m1, x_new, 1.0) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + eigvls.append(eigvl_bar) + seeds.append(0) + rates_array.append(copy.deepcopy(rates_bar)) + + # print(len(eigvls)) + # print(len(seeds)) + # print(len(rates_array)) + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i) + + return x_result, t + dt, copy.deepcopy(rates_bar), eigvl_bar \ No newline at end of file diff --git a/openmc/deplete/integrator/epc_rk4.py b/openmc/deplete/integrator/epc_rk4.py new file mode 100644 index 000000000..1347838a6 --- /dev/null +++ b/openmc/deplete/integrator/epc_rk4.py @@ -0,0 +1,166 @@ +""" The EPC-RK4 integrator. + +Implements the EPC-RK4 algorithm. +""" + +import copy +import os +import time + +from mpi4py import MPI + +from .cram import CRAM48 +from .save_results import save_results + +def epc_rk4(operator, print_out=True): + """ Performs integration of an operator using the EPC-RK4 algorithm. + + Parameters + ---------- + operator : Operator + The operator object to simulate on. + print_out : bool, optional + Whether or not to print out time. + """ + + # Save current directory + dir_home = os.getcwd() + + # Move to folder + os.makedirs(operator.settings.output_dir, exist_ok=True) + os.chdir(operator.settings.output_dir) + + # Generate initial conditions + vec = operator.initial_condition() + + n_mats = len(vec) + + t = 0.0 + + for i, dt in enumerate(operator.settings.dt_vec): + # Create vectors + x = [copy.deepcopy(vec)] + seeds = [] + eigvls = [] + rates_array = [] + + eigvl, rates, seed = operator.eval(x[0]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrix + f0 = operator.form_matrix(rates_array[0], mat) + + x_new = CRAM48(1/2 * f0, x[0][mat], dt) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + x.append(x_result) + + eigvl, rates, seed = operator.eval(x[1]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrix + f1 = operator.form_matrix(rates_array[1], mat) + + x_new = CRAM48(1/2 * f1, x[0][mat], dt) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + x.append(x_result) + + eigvl, rates, seed = operator.eval(x[2]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrix + f2 = operator.form_matrix(rates_array[2], mat) + + x_new = CRAM48(f2, x[0][mat], dt) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + x.append(x_result) + + eigvl, rates, seed = operator.eval(x[3]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrix + f0 = operator.form_matrix(rates_array[0], mat) + f1 = operator.form_matrix(rates_array[1], mat) + f2 = operator.form_matrix(rates_array[2], mat) + f3 = operator.form_matrix(rates_array[3], mat) + + x_new = CRAM48(1/6*f0 + 1/3*f1 + 1/3*f2 + 1/6*f3, x[0][mat], dt) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i) + + t += dt + vec = copy.deepcopy(x_result) + + # Perform one last simulation + x = [copy.deepcopy(vec)] + seeds = [] + eigvls = [] + rates_array = [] + eigvl, rates, seed = operator.eval(x[0]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(rates) + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t], + len(operator.settings.dt_vec)) + + # Return to origin + os.chdir(dir_home) diff --git a/openmc/deplete/integrator/leqi_cfq4_imp.py b/openmc/deplete/integrator/leqi_cfq4_imp.py new file mode 100644 index 000000000..f5b2642ff --- /dev/null +++ b/openmc/deplete/integrator/leqi_cfq4_imp.py @@ -0,0 +1,191 @@ +""" The LE/QI CFQ4 integrator. + +Implements the LE/QI Predictor-Corrector algorithm using commutator free +high order integrators. + +This algorithm is mathematically defined as: + +.. math: + y' = A(y, t) y(t) + A_m1 = A(y_n-1, t_n-1) + A_0 = A(y_n, t_n) + A_l(t) linear extrapolation of A_m1, A_0 + Integrate to t_n+1 to get y_p + A_c = A(y_p, y_n+1) + A_q(t) quadratic interpolation of A_m1, A_0, A_c + +Here, A(t) is integrated using the fourth order algorithm described below. + +From +---- + Thalhammer, Mechthild. "A fourth-order commutator-free exponential + integrator for nonautonomous differential equations." SIAM journal on + numerical analysis 44.2 (2006): 851-864. + +It is initialized using the CE/LI algorithm. +""" + +import copy +import os +import time + +from mpi4py import MPI + +from .celi_cfq4_imp import celi_cfq4_imp_inner +from .cram import CRAM48 +from .save_results import save_results + +def leqi_cfq4_imp(operator, print_out=True): + """ Performs integration of an operator using the LE/QI CFQ4 algorithm. + + Parameters + ---------- + operator : Operator + The operator object to simulate on. + print_out : bool, optional + Whether or not to print out time. + """ + + m = 10 + + # Save current directory + dir_home = os.getcwd() + + # Move to folder + os.makedirs(operator.settings.output_dir, exist_ok=True) + os.chdir(operator.settings.output_dir) + + # Generate initial conditions + vec = operator.initial_condition() + + n_mats = len(vec) + + t = 0.0 + + # Compute initial rates + operator.settings.particles *= 10 + eigvl_last, rates, seed = operator.eval(vec) + rates_last = copy.deepcopy(rates) + operator.settings.particles = int(operator.settings.particles / 10) + + # Perform single step of CE/LI CFQ4 Implicit + dt_l = operator.settings.dt_vec[0] + vec, t, rates_bos, eigvl_bos = celi_cfq4_imp_inner(operator, vec, rates_last, eigvl_last, 0, t, dt_l, print_out) + + rates_bar = [] + + # Perform remaining LE/QI + for i, dt in enumerate(operator.settings.dt_vec[1::]): + # Create vectors + x = [copy.deepcopy(vec)] + + seeds = [0] + eigvls = [eigvl_bos] + rates_array = [copy.deepcopy(rates_bos)] + + # Perform extrapolation + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrices + f1 = dt * operator.form_matrix(rates_last, mat) + f2 = dt * operator.form_matrix(rates_bos, mat) + + # Perform commutator-free integral + x_new = copy.deepcopy(x[0][mat]) + + # Compute linearly extrapolated f at points + # A{1,2} = f(1/2 -/+ sqrt(3)/6) + # Then + # a{1,2} = 1/4 +/- sqrt(3)/6 + # m1 = a2 * A1 + a1 * A2 + # m2 = a1 * A1 + a2 * A2 + m1 = -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 + m2 = -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 + + x_new = CRAM48(m2, x_new, 1.0) + x_new = CRAM48(m1, x_new, 1.0) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + x.append(copy.deepcopy(x_result)) + eigvl_bar = 0.0 + rates_bar = [] + + # Loop on inner + for j in range(0, m + 1): + eigvl, rates, seed = operator.eval(x_result) + + if j <= 1: + rates_bar = copy.deepcopy(rates) + eigvl_bar = eigvl + else: + rates_bar.rates = 1/j * rates.rates + (1 - 1/j) * rates_bar.rates + eigvl_bar = 1/j * eigvl + (1 - 1/j) * eigvl_bar + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrices + f1 = dt * operator.form_matrix(rates_last, mat) + f2 = dt * operator.form_matrix(rates_bos, mat) + f3 = dt * operator.form_matrix(rates_bar, mat) + + # Perform commutator-free integral + x_new = copy.deepcopy(x[0][mat]) + + # Compute quadratically interpolated f at points + # A{1,2} = f(1/2 -/+ sqrt(3)/6) + # Then + # a{1,2} = 1/4 +/- sqrt(3)/6 + # m1 = a2 * A1 + a1 * A2 + # m2 = a1 * A1 + a2 * A2 + m1 = (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + + (dt**2 + 2 * dt * dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + + (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3) + m2 = (-dt**2/(12 * dt_l * (dt + dt_l)) * f1 + + (dt**2 + 6 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + + dt_l / (12 * (dt + dt_l)) * f3) + + x_new = CRAM48(m2, x_new, 1.0) + x_new = CRAM48(m1, x_new, 1.0) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + eigvls.append(eigvl_bar) + seeds.append(0) + rates_array.append(copy.deepcopy(rates_bar)) + + save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i+1) + + rates_last = copy.deepcopy(rates_bos) + rates_bos = copy.deepcopy(rates_bar) + eigvl_bos = eigvl_bar + t += dt + dt_l = dt + vec = copy.deepcopy(x_result) + + # Perform one last simulation + x = [copy.deepcopy(vec)] + seeds = [0] + eigvls = [eigvl_bos] + rates_array = [copy.deepcopy(rates_bos)] + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t], + len(operator.settings.dt_vec)) + + # Return to origin + os.chdir(dir_home) From 2c1060ebaaede75153f3634aba48f230cbc1cac3 Mon Sep 17 00:00:00 2001 From: liangjg Date: Tue, 8 Jan 2019 17:17:00 -0500 Subject: [PATCH 06/24] updated epc_rk4 algorithm --- openmc/deplete/integrator/cecm.py | 12 +- openmc/deplete/integrator/cf4.py | 41 +-- openmc/deplete/integrator/epc_rk4.py | 277 +++++++++--------- openmc/deplete/integrator/predictor.py | 14 +- .../{celi_cfq4_imp.py => si_celi.py} | 0 .../{leqi_cfq4_imp.py => si_leqi.py} | 0 6 files changed, 163 insertions(+), 181 deletions(-) rename openmc/deplete/integrator/{celi_cfq4_imp.py => si_celi.py} (100%) rename openmc/deplete/integrator/{leqi_cfq4_imp.py => si_leqi.py} (100%) diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index 4554a3f22..926e27257 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -59,20 +59,16 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True): # Generate initial conditions with operator as vec: - chain = operator.chain - - # Initialize time + # Initialize time and starting index if operator.prev_res is None: t = 0.0 - else: - t = operator.prev_res[-1].time[-1] - - # Initialize starting index for saving results - if operator.prev_res is None: i_res = 0 else: + t = operator.prev_res[-1].time[-1] i_res = len(operator.prev_res) + chain = operator.chain + 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 diff --git a/openmc/deplete/integrator/cf4.py b/openmc/deplete/integrator/cf4.py index 6fdd8bd42..79cb3b3d9 100644 --- a/openmc/deplete/integrator/cf4.py +++ b/openmc/deplete/integrator/cf4.py @@ -8,19 +8,20 @@ from ..results import Results # Functions to form the special matrix for depletion -def mf1(chain, rates): +def _cf4_f1(chain, rates): return 1/2 * chain.form_matrix(rates) -def mf3(chain, rates): - return -1/2 * chain.form_matrix(rates[0]) + chain.form_matrix(rates[1]) +def _cf4_f2(chain, rates): + return -1/2 * chain.form_matrix(rates[0]) + \ + chain.form_matrix(rates[1]) -def mf4(chain, rates): +def _cf4_f3(chain, rates): return 1/4 * chain.form_matrix(rates[0]) + \ 1/6 * chain.form_matrix(rates[1]) + \ 1/6 * chain.form_matrix(rates[2]) + \ -1/12 * chain.form_matrix(rates[3]) -def mf5(chain, rates): +def _cf4_f4(chain, rates): return -1/12 * chain.form_matrix(rates[0]) + \ 1/6 * chain.form_matrix(rates[1]) + \ 1/6 * chain.form_matrix(rates[2]) + \ @@ -29,7 +30,7 @@ def mf5(chain, rates): def cf4(operator, timesteps, power=None, power_density=None, print_out=True): r"""Deplete using the CF4 algorithm. - Implements the fourth order commutator-free Lie algorithm. + Implements the fourth order [commutator-free Lie algorithm]_. This algorithm is mathematically defined as: .. math:: @@ -69,6 +70,12 @@ def cf4(operator, timesteps, power=None, power_density=None, print_out=True): print_out : bool, optional Whether or not to print out time. + References + ---------- + .. [commutator-free Lie algorithm] + Celledoni, E., Marthinsen, A. and Owren, B., 2003. Commutator-free Lie + group methods. Future Generation Computer Systems, 19(3), pp.341-352. + """ if power is None: if power_density is None: @@ -84,20 +91,16 @@ def cf4(operator, timesteps, power=None, power_density=None, print_out=True): # Generate initial conditions with operator as vec: - chain = operator.chain - - # Initialize time + # Initialize time and starting index if operator.prev_res is None: t = 0.0 - else: - t = operator.prev_res[-1].time[-1] - - # Initialize starting index for saving results - if operator.prev_res is None: i_res = 0 else: + t = operator.prev_res[-1].time[-1] i_res = len(operator.prev_res) + chain = operator.chain + 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 @@ -124,20 +127,20 @@ def cf4(operator, timesteps, power=None, power_density=None, print_out=True): # Step 1: deplete with matrix 1/2*A(y0) x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out, - matrix_func=mf1) + matrix_func=_cf4_f1) x.append(x_new) op_results.append(operator(x[1], p)) # Step 2: deplete with matrix 1/2*A(y1) x_new = deplete(chain, x[0], op_results[1].rates, dt, print_out, - matrix_func=mf1) + matrix_func=_cf4_f1) x.append(x_new) op_results.append(operator(x[2], p)) # Step 3: deplete with matrix -1/2*A(y0)+A(y2) rates = list(zip(op_results[0].rates, op_results[2].rates)) x_new = deplete(chain, x[1], rates, dt, print_out, - matrix_func=mf3) + matrix_func=_cf4_f2) x.append(x_new) op_results.append(operator(x[3], p)) @@ -145,9 +148,9 @@ def cf4(operator, timesteps, power=None, power_density=None, print_out=True): rates = list(zip(op_results[0].rates, op_results[1].rates, op_results[2].rates, op_results[3].rates)) x_end = deplete(chain, x[0], rates, dt, print_out, - matrix_func=mf4) + matrix_func=_cf4_f3) x_end = deplete(chain, x_end, rates, dt, print_out, - matrix_func=mf5) + matrix_func=_cf4_f4) # Create results, write to disk Results.save(operator, x, op_results, [t, t + dt], p, i_res + i) diff --git a/openmc/deplete/integrator/epc_rk4.py b/openmc/deplete/integrator/epc_rk4.py index 1347838a6..98e8ac58b 100644 --- a/openmc/deplete/integrator/epc_rk4.py +++ b/openmc/deplete/integrator/epc_rk4.py @@ -1,166 +1,153 @@ -""" The EPC-RK4 integrator. - -Implements the EPC-RK4 algorithm. -""" +"""The EPC-RK4 integrator.""" import copy -import os -import time +from collections.abc import Iterable -from mpi4py import MPI +from .cram import deplete +from ..results import Results -from .cram import CRAM48 -from .save_results import save_results -def epc_rk4(operator, print_out=True): - """ Performs integration of an operator using the EPC-RK4 algorithm. +# Functions to form the special matrix for depletion +def _rk4_f1(chain, rates): + return 1/2 * chain.form_matrix(rates) + +def _rk4_f4(chain, rates): + return 1/6 * chain.form_matrix(rates[0]) + \ + 1/3 * chain.form_matrix(rates[1]) + \ + 1/3 * chain.form_matrix(rates[2]) + \ + 1/6 * chain.form_matrix(rates[3]) + +def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True): + r"""Deplete using the EPC-RK4 algorithm. + + Implements an extended predictor-corrector algorithm with traditional + Runge-Kutta 4 method. + This algorithm is mathematically defined as: + + .. math:: + F_1 = h A(y_0) + + y_1 &= \text{expm}(1/2 F_1) y_0 + + F_2 &= h A(y_1) + + y_2 &= \text{expm}(1/2 F_2) y_0 + + F_3 &= h A(y_2) + + y_3 &= \text{expm}(F_3) y_0 + + F_4 &= h A(y_3) + + y_4 &= \text{expm}(1/6 F_1 + 1/3 F_2 + 1/3 F_3 + 1/6 F_4) y_0 Parameters ---------- - operator : Operator + 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] - # Save current directory - dir_home = os.getcwd() - - # Move to folder - os.makedirs(operator.settings.output_dir, exist_ok=True) - os.chdir(operator.settings.output_dir) + if not isinstance(power, Iterable): + power = [power]*len(timesteps) # Generate initial conditions - vec = operator.initial_condition() + 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) - n_mats = len(vec) + chain = operator.chain - t = 0.0 + # Initialize starting index for saving results + if operator.prev_res is None: + i_res = 0 + else: + i_res = len(operator.prev_res) - for i, dt in enumerate(operator.settings.dt_vec): - # Create vectors + 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)] + + else: + # Get initial concentration + x = [operator.prev_res[-1].data[0]] + + # 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] + + # 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] + + # Step 1: deplete with matrix 1/2*A(y0) + x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out, + matrix_func=_rk4_f1) + x.append(x_new) + op_results.append(operator(x[1], p)) + + # Step 2: deplete with matrix 1/2*A(y1) + x_new = deplete(chain, x[0], op_results[1].rates, dt, print_out, + matrix_func=_rk4_f1) + x.append(x_new) + op_results.append(operator(x[2], p)) + + # Step 3: deplete with matrix A(y2) + x_new = deplete(chain, x[0], op_results[2].rates, dt, print_out) + x.append(x_new) + op_results.append(operator(x[3], p)) + + # Step 4: deplete with two matrix exponentials + rates = list(zip(op_results[0].rates, op_results[1].rates, + op_results[2].rates, op_results[3].rates)) + x_end = deplete(chain, x[0], rates, dt, print_out, + matrix_func=_rk4_f4) + + # Create results, write to disk + Results.save(operator, x, op_results, [t, t + dt], p, i_res + i) + + # Advance time, update vector + t += dt + vec = copy.deepcopy(x_end) + + # Perform one last simulation x = [copy.deepcopy(vec)] - seeds = [] - eigvls = [] - rates_array = [] - - eigvl, rates, seed = operator.eval(x[0]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrix - f0 = operator.form_matrix(rates_array[0], mat) - - x_new = CRAM48(1/2 * f0, x[0][mat], dt) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - x.append(x_result) - - eigvl, rates, seed = operator.eval(x[1]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrix - f1 = operator.form_matrix(rates_array[1], mat) - - x_new = CRAM48(1/2 * f1, x[0][mat], dt) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - x.append(x_result) - - eigvl, rates, seed = operator.eval(x[2]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrix - f2 = operator.form_matrix(rates_array[2], mat) - - x_new = CRAM48(f2, x[0][mat], dt) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - x.append(x_result) - - eigvl, rates, seed = operator.eval(x[3]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrix - f0 = operator.form_matrix(rates_array[0], mat) - f1 = operator.form_matrix(rates_array[1], mat) - f2 = operator.form_matrix(rates_array[2], mat) - f3 = operator.form_matrix(rates_array[3], mat) - - x_new = CRAM48(1/6*f0 + 1/3*f1 + 1/3*f2 + 1/6*f3, x[0][mat], dt) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) + op_results = [operator(x[0], power[-1])] # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i) - - t += dt - vec = copy.deepcopy(x_result) - - # Perform one last simulation - x = [copy.deepcopy(vec)] - seeds = [] - eigvls = [] - rates_array = [] - eigvl, rates, seed = operator.eval(x[0]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(rates) - - # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t], - len(operator.settings.dt_vec)) - - # Return to origin - os.chdir(dir_home) + Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps)) diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index f480f55fb..71aa8fd58 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -55,19 +55,15 @@ def predictor(operator, timesteps, power=None, power_density=None, # Generate initial conditions with operator as vec: - chain = operator.chain - - # Initialize time + # Initialize time and starting index if operator.prev_res is None: t = 0.0 - else: - t = operator.prev_res[-1].time[-1] - - # Initialize starting index for saving results - if operator.prev_res is None: i_res = 0 else: - i_res = len(operator.prev_res) - 1 + t = operator.prev_res[-1].time[-1] + i_res = len(operator.prev_res) + + chain = operator.chain for i, (dt, p) in enumerate(zip(timesteps, power)): # Get beginning-of-timestep concentrations and reaction rates diff --git a/openmc/deplete/integrator/celi_cfq4_imp.py b/openmc/deplete/integrator/si_celi.py similarity index 100% rename from openmc/deplete/integrator/celi_cfq4_imp.py rename to openmc/deplete/integrator/si_celi.py diff --git a/openmc/deplete/integrator/leqi_cfq4_imp.py b/openmc/deplete/integrator/si_leqi.py similarity index 100% rename from openmc/deplete/integrator/leqi_cfq4_imp.py rename to openmc/deplete/integrator/si_leqi.py From ff158af3c873dc50d3989c9a327db09db89a92d0 Mon Sep 17 00:00:00 2001 From: liangjg Date: Tue, 8 Jan 2019 17:17:49 -0500 Subject: [PATCH 07/24] added si-ce/li algorithm --- openmc/deplete/integrator/__init__.py | 2 + openmc/deplete/integrator/si_celi.py | 280 ++++++++++++-------------- 2 files changed, 133 insertions(+), 149 deletions(-) diff --git a/openmc/deplete/integrator/__init__.py b/openmc/deplete/integrator/__init__.py index 18accdcc8..34bfb0e58 100644 --- a/openmc/deplete/integrator/__init__.py +++ b/openmc/deplete/integrator/__init__.py @@ -8,4 +8,6 @@ The integrator subcomponents. from .cf4 import * from .cecm import * from .cram import * +from .epc_rk4 import * from .predictor import * +from .si_celi import * diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py index a17f2fda9..499ae9be3 100644 --- a/openmc/deplete/integrator/si_celi.py +++ b/openmc/deplete/integrator/si_celi.py @@ -1,199 +1,181 @@ -""" The CE/LI CFQ4 integrator. - -Implements the CE/LI Predictor-Corrector algorithm using commutator free -high order integrators. - -This algorithm is mathematically defined as: - -.. math: - y' = A(y, t) y(t) - A_p = A(y_n, t_n) - y_p = expm(A_p h) y_n - A_c = A(y_p, t_n) - A(t) = t/dt * A_c + (dt - t)/dt * A_p - -Here, A(t) is integrated using the fourth order algorithm described below. - -From ----- - Thalhammer, Mechthild. "A fourth-order commutator-free exponential - integrator for nonautonomous differential equations." SIAM journal on - numerical analysis 44.2 (2006): 851-864. -""" +"""The SI-CE/LI CFQ4 integrator.""" import copy -import os -import time +from collections.abc import Iterable -from mpi4py import MPI +from .cram import deplete +from ..results import Results -from .cram import CRAM48 -from .save_results import save_results -def celi_cfq4_imp(operator, print_out=True): - """ Performs integration of an operator using the CE/LI CFQ4 algorithm. +# Functions to form the special matrix for depletion +def _celi_f1(chain, rates): + return 1/12 * chain.form_matrix(rates[0]) + \ + 5/12 * chain.form_matrix(rates[1]) + +def _celi_f2(chain, rates): + return 5/12 * chain.form_matrix(rates[0]) + \ + 1/12 * chain.form_matrix(rates[1]) + +def si_celi(operator, timesteps, power=None, power_density=None, print_out=True): + r"""Deplete using the SI-CE/LI CFQ4 algorithm. + + Implements the Stochastic Implicit CE/LI Predictor-Corrector algorithm using + the [fourth order commutator-free integrator]_. + + The CE/LI algorithm is mathematically defined as: + + .. math: + y' = A(y, t) y(t) + A_p = A(y_n, t_n) + y_p = expm(A_p h) y_n + A_c = A(y_p, t_n) + A(t) = t/dt * A_c + (dt - t)/dt * A_p + + Here, A(t) is integrated using the fourth order algorithm CFQ4. Parameters ---------- - operator : Operator + 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. + + References + ---------- + .. [fourth order commutator-free integrator] + Thalhammer, Mechthild. "A fourth-order commutator-free exponential + integrator for nonautonomous differential equations." SIAM journal on + numerical analysis 44.2 (2006): 851-864. """ + 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] - m = 10 - - # Save current directory - dir_home = os.getcwd() - - # Move to folder - os.makedirs(operator.settings.output_dir, exist_ok=True) - os.chdir(operator.settings.output_dir) + if not isinstance(power, Iterable): + power = [power]*len(timesteps) # Generate initial conditions - vec = operator.initial_condition() + 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) - # Compute initial rates - x = copy.deepcopy(vec) + op_results = None + x = [copy.deepcopy(vec)] + for i, (dt, p) in enumerate(zip(timesteps, power)): + # run the inner loop + x, t, op_results = si_celi_inner(operator, x, op_results, p + i, i_res, t, dt, print_out) - operator.settings.particles *= m - eigvl_bos, rates, seed = operator.eval(vec) - operator.settings.particles = int(operator.settings.particles / m) + # Create results for last point, write to disk + Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps)) - rates_bos = copy.deepcopy(rates) - - t = 0.0 - - for i, dt in enumerate(operator.settings.dt_vec): - vec, t, rates_eos, eigvl_bos = celi_cfq4_imp_inner(operator, vec, rates_bos, eigvl_bos, i, t, dt, print_out) - rates_bos = copy.deepcopy(rates_eos) - - # Perform one last simulation - x = [copy.deepcopy(vec)] - seeds = [0] - eigvls = [eigvl_bos] - rates_array = [copy.deepcopy(rates_bos)] - - # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t], - len(operator.settings.dt_vec)) - - # Return to origin - os.chdir(dir_home) - -def celi_cfq4_imp_inner(operator, vec, rates_bos, eigvl_bos, i, t, dt, print_out): - """ The inner loop of CE/LI CFQ4. +def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out): + """ The inner loop of SI-CE/LI CFQ4. Parameters ---------- operator : Operator The operator object to simulate on. - vec : list of numpy.array + x : list of nuclide vector Nuclide vector, beginning of time. - i : Int + op_results : list of OperatorResult + Operator result at BOS. + p : float + Power of the reactor in [W] + i : int Current iteration number. - t : Float + i_res : int + Starting index, for restart calculation. + t : float Time at start of step. - dt : Float + dt : float Time step. print_out : bool Whether or not to print out time. Returns ------- - x_result : list of numpy.array + list of nuclide vector (numpy.array) Nuclide vector, end of time. - Float + float Next time - ReactionRates - Reaction rates from beginning of step. + list of OperatorResult + Operator result at end of time. """ - m = 10 + m = 10 # stage number - n_mats = len(vec) + # Get the concentrations and reaction rates for the first + # beginning-of-timestep (BOS) + # Compute with s (stage number) times as many neutrons for statistics + # reasons if no previous calculation results loaded + if i == 0: + if operator.prev_res is None: + operator.settings.particles *= m + op_results = [operator(x[0], p)] + operator.settings.particles //= m + else: + # Get initial concentration + x = [operator.prev_res[-1].data[0]] - # Create vectors - x = [copy.deepcopy(vec)] - seeds = [] - eigvls = [] - rates_array = [] + # Get rates + op_results = [operator.prev_res[-1]] + op_results[0].rates = op_results[0].rates[0] - eigvls.append(eigvl_bos) - seeds.append(0) - rates_array.append(copy.deepcopy(rates_bos)) + # Set first stage value of keff + op_results[0].k = op_results[0].k[0] + + # 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] + + chain = operator.chain # Deplete to end - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrix - f = operator.form_matrix(rates_array[0], mat) - - x_new = CRAM48(f, x[0][mat], dt) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - x.append(copy.deepcopy(x_result)) - - eigvl_bar = 0.0 - rates_bar = [] + x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out) + x.append(x_new) for j in range(0, m + 1): - eigvl, rates, seed = operator.eval(x_result) + op_res = operator.eval(x_new, p) if j <= 1: - rates_bar = copy.deepcopy(rates) - eigvl_bar = eigvl + op_res_bar = copy.deepcopy(op_res) else: - rates_bar.rates = 1/j * rates.rates + (1 - 1/j) * rates_bar.rates - eigvl_bar = 1/j * eigvl + (1 - 1/j) * eigvl_bar + op_res_bar.rates = 1/j * op_res.rates + (1 - 1/j) * op_res_bar.rates + op_res_bar.k = 1/j * op_res.k + (1 - 1/j) * op_res_bar.k - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrices - f1 = dt * operator.form_matrix(rates_array[0], mat) - f2 = dt * operator.form_matrix(rates_bar, mat) - - # Perform commutator-free integral - x_new = copy.deepcopy(x[0][mat]) - - # Compute linearly interpolated f at points - # A{1,2} = f(1/2 -/+ sqrt(3)/6) - # Then - # a{1,2} = 1/4 +/- sqrt(3)/6 - # m1 = a2 * A1 + a1 * A2 - # m2 = a1 * A1 + a2 * A2 - m1 = 1/12 * (f1 + 5 * f2) - m2 = 1/12 * (5 * f1 + f2) - - x_new = CRAM48(m2, x_new, 1.0) - x_new = CRAM48(m1, x_new, 1.0) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - eigvls.append(eigvl_bar) - seeds.append(0) - rates_array.append(copy.deepcopy(rates_bar)) - - # print(len(eigvls)) - # print(len(seeds)) - # print(len(rates_array)) + rates = list(zip(op_results[0].rates, op_res_bar.rates)) + x_new = deplete(chain, x[0], rates, dt, print_out, + matrix_func=_celi_f1) + x_new = deplete(chain, x_new, rates, dt, print_out, + matrix_func=_celi_f2) # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i) + op_results.append(op_res_bar) + Results.save(operator, x, op_results, [t, t + dt], p, i+i_res) - return x_result, t + dt, copy.deepcopy(rates_bar), eigvl_bar \ No newline at end of file + # return updated time and vectors + return [copy.deepcopy(x_new)], t + dt, [copy.deepcopy(op_res_bar)] From a594ff414bd5d3e4b77baeaa0c362910de0e4c5c Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 9 Jan 2019 01:00:03 -0500 Subject: [PATCH 08/24] added SI-LEQI CFQ4 algorithm --- openmc/deplete/integrator/__init__.py | 1 + openmc/deplete/integrator/si_celi.py | 82 +++---- openmc/deplete/integrator/si_leqi.py | 327 +++++++++++++------------- 3 files changed, 205 insertions(+), 205 deletions(-) diff --git a/openmc/deplete/integrator/__init__.py b/openmc/deplete/integrator/__init__.py index 34bfb0e58..4dad78e1c 100644 --- a/openmc/deplete/integrator/__init__.py +++ b/openmc/deplete/integrator/__init__.py @@ -11,3 +11,4 @@ from .cram import * from .epc_rk4 import * from .predictor import * from .si_celi import * +from .si_leqi import * diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py index 499ae9be3..8d3cc0392 100644 --- a/openmc/deplete/integrator/si_celi.py +++ b/openmc/deplete/integrator/si_celi.py @@ -5,17 +5,21 @@ from collections.abc import Iterable from .cram import deplete from ..results import Results +from ..abc import OperatorResult +# Stage number for CE/LI +_CELI_STAGE_M = 10 + # Functions to form the special matrix for depletion def _celi_f1(chain, rates): - return 1/12 * chain.form_matrix(rates[0]) + \ - 5/12 * chain.form_matrix(rates[1]) - -def _celi_f2(chain, rates): return 5/12 * chain.form_matrix(rates[0]) + \ 1/12 * chain.form_matrix(rates[1]) +def _celi_f2(chain, rates): + return 1/12 * chain.form_matrix(rates[0]) + \ + 5/12 * chain.form_matrix(rates[1]) + def si_celi(operator, timesteps, power=None, power_density=None, print_out=True): r"""Deplete using the SI-CE/LI CFQ4 algorithm. @@ -81,11 +85,33 @@ def si_celi(operator, timesteps, power=None, power_density=None, print_out=True) t = operator.prev_res[-1].time[-1] i_res = len(operator.prev_res) - op_results = None - x = [copy.deepcopy(vec)] + # Get the concentrations and reaction rates for the first + # beginning-of-timestep (BOS) + # Compute with s (stage number) times as many neutrons for statistics + # reasons if no previous calculation results loaded + if operator.prev_res is None: + x = [copy.deepcopy(vec)] + operator.settings.particles *= _CELI_STAGE_M + op_results = [operator(x[0], power[0])] + operator.settings.particles //= _CELI_STAGE_M + else: + # Get initial concentration + x = [operator.prev_res[-1].data[0]] + + # 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] + + # 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] + for i, (dt, p) in enumerate(zip(timesteps, power)): - # run the inner loop - x, t, op_results = si_celi_inner(operator, x, op_results, p + x, t, op_results = si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out) # Create results for last point, write to disk @@ -125,47 +151,21 @@ def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out): Operator result at end of time. """ - m = 10 # stage number - - # Get the concentrations and reaction rates for the first - # beginning-of-timestep (BOS) - # Compute with s (stage number) times as many neutrons for statistics - # reasons if no previous calculation results loaded - if i == 0: - if operator.prev_res is None: - operator.settings.particles *= m - op_results = [operator(x[0], p)] - operator.settings.particles //= m - else: - # Get initial concentration - x = [operator.prev_res[-1].data[0]] - - # 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] - - # 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] - chain = operator.chain # Deplete to end x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out) x.append(x_new) - for j in range(0, m + 1): - op_res = operator.eval(x_new, p) + for j in range(_CELI_STAGE_M + 1): + op_res = operator(x_new, p) if j <= 1: op_res_bar = copy.deepcopy(op_res) else: - op_res_bar.rates = 1/j * op_res.rates + (1 - 1/j) * op_res_bar.rates - op_res_bar.k = 1/j * op_res.k + (1 - 1/j) * op_res_bar.k + 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)) x_new = deplete(chain, x[0], rates, dt, print_out, @@ -175,7 +175,7 @@ def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out): # Create results, write to disk op_results.append(op_res_bar) - Results.save(operator, x, op_results, [t, t + dt], p, i+i_res) + Results.save(operator, x, op_results, [t, t+dt], p, i_res+i) # return updated time and vectors - return [copy.deepcopy(x_new)], t + dt, [copy.deepcopy(op_res_bar)] + return [x_new], t + dt, [op_res_bar] diff --git a/openmc/deplete/integrator/si_leqi.py b/openmc/deplete/integrator/si_leqi.py index f5b2642ff..7b3cb34c4 100644 --- a/openmc/deplete/integrator/si_leqi.py +++ b/openmc/deplete/integrator/si_leqi.py @@ -1,191 +1,190 @@ -""" The LE/QI CFQ4 integrator. - -Implements the LE/QI Predictor-Corrector algorithm using commutator free -high order integrators. - -This algorithm is mathematically defined as: - -.. math: - y' = A(y, t) y(t) - A_m1 = A(y_n-1, t_n-1) - A_0 = A(y_n, t_n) - A_l(t) linear extrapolation of A_m1, A_0 - Integrate to t_n+1 to get y_p - A_c = A(y_p, y_n+1) - A_q(t) quadratic interpolation of A_m1, A_0, A_c - -Here, A(t) is integrated using the fourth order algorithm described below. - -From ----- - Thalhammer, Mechthild. "A fourth-order commutator-free exponential - integrator for nonautonomous differential equations." SIAM journal on - numerical analysis 44.2 (2006): 851-864. - -It is initialized using the CE/LI algorithm. -""" +"""The SI-LE/QI CFQ4 integrator.""" import copy -import os -import time +from collections.abc import Iterable +from itertools import repeat -from mpi4py import MPI +from .si_celi import si_celi_inner +from .cram import deplete +from ..results import Results +from ..abc import OperatorResult -from .celi_cfq4_imp import celi_cfq4_imp_inner -from .cram import CRAM48 -from .save_results import save_results -def leqi_cfq4_imp(operator, print_out=True): - """ Performs integration of an operator using the LE/QI CFQ4 algorithm. +# Stage number for CE/LI +_LEQI_STAGE_M = 10 + +# Functions to form the special matrix for depletion +def _leqi_f1(chain, inputs): + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + dt_l, dt = inputs[2], inputs[3] + return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 + +def _leqi_f2(chain, inputs): + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + dt_l, dt = inputs[2], inputs[3] + return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 + +def _leqi_f3(chain, inputs): + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + f3 = chain.form_matrix(inputs[2]) + dt_l, dt = inputs[3], inputs[4] + return (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + + (dt**2 + 6*dt*dt_l + 5*dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + + dt_l / (12 * (dt + dt_l)) * f3) + +def _leqi_f4(chain, inputs): + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + f3 = chain.form_matrix(inputs[2]) + dt_l, dt = inputs[3], inputs[4] + return (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + + (dt**2 + 2*dt*dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + + (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3) + +def si_leqi(operator, timesteps, power=None, power_density=None, print_out=True): + r"""Deplete using the SI-LE/QI CFQ4 algorithm. + + Implements the Stochastic Implicit LE/QI Predictor-Corrector algorithm using + the [fourth order commutator-free integrator]_. + + The LE/QI algorithm is mathematically defined as: + + .. math: + y' = A(y, t) y(t) + A_m1 = A(y_n-1, t_n-1) + A_0 = A(y_n, t_n) + A_l(t) linear extrapolation of A_m1, A_0 + Integrate to t_n+1 to get y_p + A_c = A(y_p, y_n+1) + A_q(t) quadratic interpolation of A_m1, A_0, A_c + + Here, A(t) is integrated using the fourth order algorithm CFQ4. + + It is initialized using the CE/LI algorithm. Parameters ---------- - operator : Operator + 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. + + References + ---------- + .. [fourth order commutator-free integrator] + Thalhammer, Mechthild. "A fourth-order commutator-free exponential + integrator for nonautonomous differential equations." SIAM journal on + numerical analysis 44.2 (2006): 851-864. """ + 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] - m = 10 - - # Save current directory - dir_home = os.getcwd() - - # Move to folder - os.makedirs(operator.settings.output_dir, exist_ok=True) - os.chdir(operator.settings.output_dir) + if not isinstance(power, Iterable): + power = [power]*len(timesteps) # Generate initial conditions - vec = operator.initial_condition() + 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) - n_mats = len(vec) + # Get the concentrations and reaction rates for the first + # beginning-of-timestep (BOS) + # Compute with s (stage number) times as many neutrons for statistics + # reasons if no previous calculation results loaded + if operator.prev_res is None: + x = [copy.deepcopy(vec)] + operator.settings.particles *= _LEQI_STAGE_M + op_results = [operator(x[0], power[0])] + operator.settings.particles //= _LEQI_STAGE_M + else: + # Get initial concentration + x = [operator.prev_res[-1].data[0]] - t = 0.0 + # Get rates + op_results = [operator.prev_res[-1]] + op_results[0].rates = op_results[0].rates[0] - # Compute initial rates - operator.settings.particles *= 10 - eigvl_last, rates, seed = operator.eval(vec) - rates_last = copy.deepcopy(rates) - operator.settings.particles = int(operator.settings.particles / 10) + # Set first stage value of keff + op_results[0].k = op_results[0].k[0] - # Perform single step of CE/LI CFQ4 Implicit - dt_l = operator.settings.dt_vec[0] - vec, t, rates_bos, eigvl_bos = celi_cfq4_imp_inner(operator, vec, rates_last, eigvl_last, 0, t, dt_l, print_out) + # 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] - rates_bar = [] + for i, (dt, p) in enumerate(zip(timesteps, power)): + # Perform SI-CE/LI CFQ4 for the first step + if i == 0: + # Save results for the last step + op_res_last = copy.deepcopy(op_results[0]) + dt_l = dt + x, t, op_results = si_celi_inner(operator, x, op_results, p, + i, i_res, t, dt, print_out) + continue - # Perform remaining LE/QI - for i, dt in enumerate(operator.settings.dt_vec[1::]): - # Create vectors - x = [copy.deepcopy(vec)] + # Perform remaining LE/QI + inputs = list(zip(op_res_last.rates, op_results[0].rates, + repeat(dt_l), repeat(dt))) + x_new = deplete(chain, x[0], inputs, dt, print_out, + matrix_func=_leqi_f1) + x_new = deplete(chain, x_new, inputs, dt, print_out, + matrix_func=_leqi_f2) + x.append(x_new) - seeds = [0] - eigvls = [eigvl_bos] - rates_array = [copy.deepcopy(rates_bos)] + # Loop on inner + for j in range(_LEQI_STAGE_M + 1): + op_res = operator(x_new, p) - # Perform extrapolation - x_result = [] + 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) - t_start = time.time() - for mat in range(n_mats): - # Form matrices - f1 = dt * operator.form_matrix(rates_last, mat) - f2 = dt * operator.form_matrix(rates_bos, mat) + inputs = list(zip(op_res_last.rates, op_results[0].rates, + op_res_bar.rate, repeat(dt_l), repeat(dt))) + x_new = deplete(chain, x[0], rates, dt, print_out, + matrix_func=_leqi_f3) + x_new = deplete(chain, x_new, rates, dt, print_out, + matrix_func=_leqi_f4) - # Perform commutator-free integral - x_new = copy.deepcopy(x[0][mat]) + # Create results, write to disk + op_results.append(op_res_bar) + Results.save(operator, x, op_results, [t, t+dt], p, i_res+i) - # Compute linearly extrapolated f at points - # A{1,2} = f(1/2 -/+ sqrt(3)/6) - # Then - # a{1,2} = 1/4 +/- sqrt(3)/6 - # m1 = a2 * A1 + a1 * A2 - # m2 = a1 * A1 + a2 * A2 - m1 = -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 - m2 = -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 + # update results + x = [x_new] + op_res_last = copy.deepcopy(op_results[0]) + op_results = [op_res_bar] + t += dt + dt_l = dt - x_new = CRAM48(m2, x_new, 1.0) - x_new = CRAM48(m1, x_new, 1.0) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - x.append(copy.deepcopy(x_result)) - eigvl_bar = 0.0 - rates_bar = [] - - # Loop on inner - for j in range(0, m + 1): - eigvl, rates, seed = operator.eval(x_result) - - if j <= 1: - rates_bar = copy.deepcopy(rates) - eigvl_bar = eigvl - else: - rates_bar.rates = 1/j * rates.rates + (1 - 1/j) * rates_bar.rates - eigvl_bar = 1/j * eigvl + (1 - 1/j) * eigvl_bar - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrices - f1 = dt * operator.form_matrix(rates_last, mat) - f2 = dt * operator.form_matrix(rates_bos, mat) - f3 = dt * operator.form_matrix(rates_bar, mat) - - # Perform commutator-free integral - x_new = copy.deepcopy(x[0][mat]) - - # Compute quadratically interpolated f at points - # A{1,2} = f(1/2 -/+ sqrt(3)/6) - # Then - # a{1,2} = 1/4 +/- sqrt(3)/6 - # m1 = a2 * A1 + a1 * A2 - # m2 = a1 * A1 + a2 * A2 - m1 = (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + - (dt**2 + 2 * dt * dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + - (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3) - m2 = (-dt**2/(12 * dt_l * (dt + dt_l)) * f1 + - (dt**2 + 6 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + - dt_l / (12 * (dt + dt_l)) * f3) - - x_new = CRAM48(m2, x_new, 1.0) - x_new = CRAM48(m1, x_new, 1.0) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - eigvls.append(eigvl_bar) - seeds.append(0) - rates_array.append(copy.deepcopy(rates_bar)) - - save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i+1) - - rates_last = copy.deepcopy(rates_bos) - rates_bos = copy.deepcopy(rates_bar) - eigvl_bos = eigvl_bar - t += dt - dt_l = dt - vec = copy.deepcopy(x_result) - - # Perform one last simulation - x = [copy.deepcopy(vec)] - seeds = [0] - eigvls = [eigvl_bos] - rates_array = [copy.deepcopy(rates_bos)] - - # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t], - len(operator.settings.dt_vec)) - - # Return to origin - os.chdir(dir_home) + # Create results for last point, write to disk + Results.save(operator, x, op_results, [t, t], p, i_res+len(timesteps)) From 87aef6dbaba49888b64bcfe01a2f4f204a2284e6 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 9 Jan 2019 01:22:07 -0500 Subject: [PATCH 09/24] allow to specify stage number for si_celi and si_leqi --- openmc/deplete/integrator/si_celi.py | 22 ++++++++++++---------- openmc/deplete/integrator/si_leqi.py | 24 +++++++++++++----------- 2 files changed, 25 insertions(+), 21 deletions(-) diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py index 8d3cc0392..6e2d34bac 100644 --- a/openmc/deplete/integrator/si_celi.py +++ b/openmc/deplete/integrator/si_celi.py @@ -8,9 +8,6 @@ from ..results import Results from ..abc import OperatorResult -# Stage number for CE/LI -_CELI_STAGE_M = 10 - # Functions to form the special matrix for depletion def _celi_f1(chain, rates): return 5/12 * chain.form_matrix(rates[0]) + \ @@ -20,7 +17,8 @@ def _celi_f2(chain, rates): return 1/12 * chain.form_matrix(rates[0]) + \ 5/12 * chain.form_matrix(rates[1]) -def si_celi(operator, timesteps, power=None, power_density=None, print_out=True): +def si_celi(operator, timesteps, power=None, power_density=None, + print_out=True, m=10): r"""Deplete using the SI-CE/LI CFQ4 algorithm. Implements the Stochastic Implicit CE/LI Predictor-Corrector algorithm using @@ -55,6 +53,8 @@ def si_celi(operator, timesteps, power=None, power_density=None, print_out=True) 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. References ---------- @@ -91,9 +91,9 @@ def si_celi(operator, timesteps, power=None, power_density=None, print_out=True) # reasons if no previous calculation results loaded if operator.prev_res is None: x = [copy.deepcopy(vec)] - operator.settings.particles *= _CELI_STAGE_M + operator.settings.particles *= m op_results = [operator(x[0], power[0])] - operator.settings.particles //= _CELI_STAGE_M + operator.settings.particles //= m else: # Get initial concentration x = [operator.prev_res[-1].data[0]] @@ -107,17 +107,17 @@ def si_celi(operator, timesteps, power=None, power_density=None, print_out=True) # Scale reaction rates by ratio of powers power_res = operator.prev_res[-1].power - ratio_power = p / power_res + 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) + 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)) -def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out): +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 @@ -140,6 +140,8 @@ def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out): Time step. print_out : bool Whether or not to print out time. + m : int, optional + Number of stages. Returns ------- @@ -157,7 +159,7 @@ def si_celi_inner(operator, x, op_results, p, i, i_res, t, dt, print_out): x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out) x.append(x_new) - for j in range(_CELI_STAGE_M + 1): + for j in range(m + 1): op_res = operator(x_new, p) if j <= 1: diff --git a/openmc/deplete/integrator/si_leqi.py b/openmc/deplete/integrator/si_leqi.py index 7b3cb34c4..d3a0ff6af 100644 --- a/openmc/deplete/integrator/si_leqi.py +++ b/openmc/deplete/integrator/si_leqi.py @@ -10,9 +10,6 @@ from ..results import Results from ..abc import OperatorResult -# Stage number for CE/LI -_LEQI_STAGE_M = 10 - # Functions to form the special matrix for depletion def _leqi_f1(chain, inputs): f1 = chain.form_matrix(inputs[0]) @@ -44,7 +41,8 @@ def _leqi_f4(chain, inputs): (dt**2 + 2*dt*dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3) -def si_leqi(operator, timesteps, power=None, power_density=None, print_out=True): +def si_leqi(operator, timesteps, power=None, power_density=None, + print_out=True, m=10): r"""Deplete using the SI-LE/QI CFQ4 algorithm. Implements the Stochastic Implicit LE/QI Predictor-Corrector algorithm using @@ -83,6 +81,8 @@ def si_leqi(operator, timesteps, power=None, power_density=None, print_out=True) 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. References ---------- @@ -119,9 +119,9 @@ def si_leqi(operator, timesteps, power=None, power_density=None, print_out=True) # reasons if no previous calculation results loaded if operator.prev_res is None: x = [copy.deepcopy(vec)] - operator.settings.particles *= _LEQI_STAGE_M + operator.settings.particles *= m op_results = [operator(x[0], power[0])] - operator.settings.particles //= _LEQI_STAGE_M + operator.settings.particles //= m else: # Get initial concentration x = [operator.prev_res[-1].data[0]] @@ -135,9 +135,11 @@ def si_leqi(operator, timesteps, power=None, power_density=None, print_out=True) # Scale reaction rates by ratio of powers power_res = operator.prev_res[-1].power - ratio_power = p / power_res + 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)): # Perform SI-CE/LI CFQ4 for the first step if i == 0: @@ -158,7 +160,7 @@ def si_leqi(operator, timesteps, power=None, power_density=None, print_out=True) x.append(x_new) # Loop on inner - for j in range(_LEQI_STAGE_M + 1): + for j in range(m + 1): op_res = operator(x_new, p) if j <= 1: @@ -169,10 +171,10 @@ def si_leqi(operator, timesteps, power=None, power_density=None, print_out=True) op_res_bar = OperatorResult(k, rates) inputs = list(zip(op_res_last.rates, op_results[0].rates, - op_res_bar.rate, repeat(dt_l), repeat(dt))) - x_new = deplete(chain, x[0], rates, dt, print_out, + op_res_bar.rates, repeat(dt_l), repeat(dt))) + x_new = deplete(chain, x[0], inputs, dt, print_out, matrix_func=_leqi_f3) - x_new = deplete(chain, x_new, rates, dt, print_out, + x_new = deplete(chain, x_new, inputs, dt, print_out, matrix_func=_leqi_f4) # Create results, write to disk From 952bddb61ec335b4e71c6ee0e16ac0bf080d988f Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 9 Jan 2019 09:40:45 -0500 Subject: [PATCH 10/24] update doc of openmc.deplete and openmc.Model.deplete --- docs/source/pythonapi/deplete.rst | 8 ++++++-- openmc/model/model.py | 21 ++++++--------------- 2 files changed, 12 insertions(+), 17 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index d4055f0fd..e38c00e96 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -6,8 +6,8 @@ .. module:: openmc.deplete -Two functions are provided that implement different time-integration algorithms -for depletion calculations. +Several functions are provided that implement different time-integration +algorithms for depletion calculations. .. autosummary:: :toctree: generated @@ -16,6 +16,10 @@ for depletion calculations. integrator.predictor integrator.cecm + integrator.cf4 + integrator.epc_rk4 + integrator.si_celi + integrator.si_leqi Each of these functions expects a "transport operator" to be passed. An operator specific to OpenMC is available using the following class: diff --git a/openmc/model/model.py b/openmc/model/model.py index 28d19713c..b98133208 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -122,7 +122,7 @@ class Model(object): for plot in plots: self._plots.append(plot) - def deplete(self, timesteps, power, chain_file=None, method='cecm', + def deplete(self, timesteps, chain_file=None, method='cecm', **kwargs): """Deplete model using specified timesteps/power @@ -131,17 +131,11 @@ class Model(object): timesteps : iterable of float Array of timesteps in units of [s]. Note that values are not cumulative. - 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]. chain_file : str, optional Path to the depletion chain XML file. Defaults to the :envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists. - method : {'cecm', 'predictor'} - Integration method used for depletion + method : string + Integration method used for depletion (e.g., 'cecm', 'predictor') **kwargs Keyword arguments passed to integration function (e.g., :func:`openmc.deplete.integrator.cecm`) @@ -156,12 +150,9 @@ class Model(object): op = dep.Operator(self.geometry, self.settings, chain_file) # Perform depletion - if method == 'predictor': - dep.integrator.predictor(op, timesteps, power, **kwargs) - elif method == 'cecm': - dep.integrator.cecm(op, timesteps, power, **kwargs) - else: - check_value('method', method, ('cecm', 'predictor')) + check_value('method', method, ('cecm', 'predictor', 'cf4', 'epc_rk4', + 'si_celi', 'si_leqi')) + getattr(dep.integrator, method)(op, timesteps, **kwargs) def export_to_xml(self): """Export model to XML files.""" From 57ea52a46ce9eb93d84d2bb5aff756b83f778574 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 9 Jan 2019 10:39:40 -0500 Subject: [PATCH 11/24] added unit tests for the new integrators --- openmc/deplete/integrator/si_celi.py | 12 ++++---- openmc/deplete/integrator/si_leqi.py | 12 ++++---- tests/unit_tests/test_deplete_cf4.py | 37 ++++++++++++++++++++++++ tests/unit_tests/test_deplete_epc_rk4.py | 37 ++++++++++++++++++++++++ tests/unit_tests/test_deplete_si_celi.py | 37 ++++++++++++++++++++++++ tests/unit_tests/test_deplete_si_leqi.py | 37 ++++++++++++++++++++++++ 6 files changed, 162 insertions(+), 10 deletions(-) create mode 100644 tests/unit_tests/test_deplete_cf4.py create mode 100644 tests/unit_tests/test_deplete_epc_rk4.py create mode 100644 tests/unit_tests/test_deplete_si_celi.py create mode 100644 tests/unit_tests/test_deplete_si_leqi.py diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py index 6e2d34bac..dc9bccb70 100644 --- a/openmc/deplete/integrator/si_celi.py +++ b/openmc/deplete/integrator/si_celi.py @@ -86,14 +86,16 @@ def si_celi(operator, timesteps, power=None, power_density=None, i_res = len(operator.prev_res) # Get the concentrations and reaction rates for the first - # beginning-of-timestep (BOS) - # Compute with s (stage number) times as many neutrons for statistics - # reasons if no previous calculation results loaded + # 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)] - operator.settings.particles *= m + if hasattr(operator, "settings"): + operator.settings.particles *= m op_results = [operator(x[0], power[0])] - operator.settings.particles //= m + if hasattr(operator, "settings"): + operator.settings.particles //= m else: # Get initial concentration x = [operator.prev_res[-1].data[0]] diff --git a/openmc/deplete/integrator/si_leqi.py b/openmc/deplete/integrator/si_leqi.py index d3a0ff6af..609b3f8c4 100644 --- a/openmc/deplete/integrator/si_leqi.py +++ b/openmc/deplete/integrator/si_leqi.py @@ -114,14 +114,16 @@ def si_leqi(operator, timesteps, power=None, power_density=None, i_res = len(operator.prev_res) # Get the concentrations and reaction rates for the first - # beginning-of-timestep (BOS) - # Compute with s (stage number) times as many neutrons for statistics - # reasons if no previous calculation results loaded + # 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)] - operator.settings.particles *= m + if hasattr(operator, "settings"): + operator.settings.particles *= m op_results = [operator(x[0], power[0])] - operator.settings.particles //= m + if hasattr(operator, "settings"): + operator.settings.particles //= m else: # Get initial concentration x = [operator.prev_res[-1].data[0]] diff --git a/tests/unit_tests/test_deplete_cf4.py b/tests/unit_tests/test_deplete_cf4.py new file mode 100644 index 000000000..784d228be --- /dev/null +++ b/tests/unit_tests/test_deplete_cf4.py @@ -0,0 +1,37 @@ +"""Regression tests for openmc.deplete.integrator.cf4 algorithm. + +These tests integrate a simple test problem described in dummy_geometry.py. +""" + +from pytest import approx +import openmc.deplete + +from tests import dummy_operator + + +def test_cf4(run_in_tmpdir): + """Integral regression test of integrator algorithm using CF4""" + + op = dummy_operator.DummyOperator() + op.output_dir = "test_integrator_regression" + + # Perform simulation using the cf4 algorithm + dt = [0.75, 0.75] + power = 1.0 + openmc.deplete.cf4(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [2.06101629, 1.37783588] + s2 = [2.57241318, 2.63731630] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) diff --git a/tests/unit_tests/test_deplete_epc_rk4.py b/tests/unit_tests/test_deplete_epc_rk4.py new file mode 100644 index 000000000..dd15bf8e1 --- /dev/null +++ b/tests/unit_tests/test_deplete_epc_rk4.py @@ -0,0 +1,37 @@ +"""Regression tests for openmc.deplete.integrator.epc_rk4 algorithm. + +These tests integrate a simple test problem described in dummy_geometry.py. +""" + +from pytest import approx +import openmc.deplete + +from tests import dummy_operator + + +def test_epc_rk4(run_in_tmpdir): + """Integral regression test of integrator algorithm using epc_rk4""" + + op = dummy_operator.DummyOperator() + op.output_dir = "test_integrator_regression" + + # Perform simulation using the epc_rk4 algorithm + dt = [0.75, 0.75] + power = 1.0 + openmc.deplete.epc_rk4(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [2.01978516, 1.42038037] + s2 = [2.05246421, 3.06177191] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) diff --git a/tests/unit_tests/test_deplete_si_celi.py b/tests/unit_tests/test_deplete_si_celi.py new file mode 100644 index 000000000..8e84bcba8 --- /dev/null +++ b/tests/unit_tests/test_deplete_si_celi.py @@ -0,0 +1,37 @@ +"""Regression tests for openmc.deplete.integrator.si_celi algorithm. + +These tests integrate a simple test problem described in dummy_geometry.py. +""" + +from pytest import approx +import openmc.deplete + +from tests import dummy_operator + + +def test_si_celi(run_in_tmpdir): + """Integral regression test of integrator algorithm using si_celi""" + + op = dummy_operator.DummyOperator() + op.output_dir = "test_integrator_regression" + + # 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) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [2.03325094, 1.16826254] + s2 = [2.69291933, 0.37907772] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) diff --git a/tests/unit_tests/test_deplete_si_leqi.py b/tests/unit_tests/test_deplete_si_leqi.py new file mode 100644 index 000000000..4814d8ab8 --- /dev/null +++ b/tests/unit_tests/test_deplete_si_leqi.py @@ -0,0 +1,37 @@ +"""Regression tests for openmc.deplete.integrator.si_leqi algorithm. + +These tests integrate a simple test problem described in dummy_geometry.py. +""" + +from pytest import approx +import openmc.deplete + +from tests import dummy_operator + + +def test_si_leqi(run_in_tmpdir): + """Integral regression test of integrator algorithm using si_leqi""" + + op = dummy_operator.DummyOperator() + op.output_dir = "test_integrator_regression" + + # 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) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [2.03325094, 1.16826254] + s2 = [2.92711288, 0.53753236] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) From cddfbf41c0d39687533130c99bea6b7e2174542b Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 9 Jan 2019 13:40:58 -0500 Subject: [PATCH 12/24] updated test results --- openmc/deplete/integrator/predictor.py | 2 +- tests/unit_tests/test_deplete_restart.py | 156 +++++++++++++++++++++++ 2 files changed, 157 insertions(+), 1 deletion(-) diff --git a/openmc/deplete/integrator/predictor.py b/openmc/deplete/integrator/predictor.py index 71aa8fd58..f670a125e 100644 --- a/openmc/deplete/integrator/predictor.py +++ b/openmc/deplete/integrator/predictor.py @@ -61,7 +61,7 @@ def predictor(operator, timesteps, power=None, power_density=None, i_res = 0 else: t = operator.prev_res[-1].time[-1] - i_res = len(operator.prev_res) + i_res = len(operator.prev_res) - 1 chain = operator.chain diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py index 745f2cce5..9c7e6c6c6 100644 --- a/tests/unit_tests/test_deplete_restart.py +++ b/tests/unit_tests/test_deplete_restart.py @@ -166,3 +166,159 @@ def test_restart_cecm_predictor(run_in_tmpdir): assert y1[2] == approx(s2[0]) assert y2[2] == approx(s2[1]) + + +def test_restart_cf4(run_in_tmpdir): + """Integral regression test of integrator algorithm using CF4.""" + + op = dummy_operator.DummyOperator() + output_dir = "test_restart_cf4" + op.output_dir = output_dir + + # Perform simulation + dt = [0.75] + power = 1.0 + openmc.deplete.cf4(op, dt, power, print_out=False) + + # Load the files + prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Re-create depletion operator and load previous results + op = dummy_operator.DummyOperator(prev_res) + op.output_dir = output_dir + + # Perform restarts simulation + openmc.deplete.cf4(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [2.06101629, 1.37783588] + s2 = [2.57241318, 2.63731630] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[3] == approx(s2[0]) + assert y2[3] == approx(s2[1]) + + +def test_restart_epc_rk4(run_in_tmpdir): + """Integral regression test of integrator algorithm using EPC-RK4.""" + + op = dummy_operator.DummyOperator() + output_dir = "test_restart_epc_rk4" + op.output_dir = output_dir + + # Perform simulation + dt = [0.75] + power = 1.0 + openmc.deplete.epc_rk4(op, dt, power, print_out=False) + + # Load the files + prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Re-create depletion operator and load previous results + op = dummy_operator.DummyOperator(prev_res) + op.output_dir = output_dir + + # Perform restarts simulation + openmc.deplete.epc_rk4(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [2.01978516, 1.42038037] + s2 = [2.05246421, 3.06177191] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[3] == approx(s2[0]) + assert y2[3] == approx(s2[1]) + + +def test_restart_si_celi(run_in_tmpdir): + """Integral regression test of integrator algorithm using SI-CELI.""" + + op = dummy_operator.DummyOperator() + output_dir = "test_restart_si_celi" + op.output_dir = output_dir + + # Perform simulation + dt = [0.75] + power = 1.0 + openmc.deplete.si_celi(op, dt, power, print_out=False) + + # Load the files + prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Re-create depletion operator and load previous results + op = dummy_operator.DummyOperator(prev_res) + op.output_dir = output_dir + + # Perform restarts simulation + openmc.deplete.si_celi(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [2.03325094, 1.16826254] + s2 = [2.69291933, 0.37907772] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[3] == approx(s2[0]) + assert y2[3] == approx(s2[1]) + + +def test_restart_si_leqi(run_in_tmpdir): + """Integral regression test of integrator algorithm using SI-LEQI.""" + + op = dummy_operator.DummyOperator() + output_dir = "test_restart_si_leqi" + op.output_dir = output_dir + + # Perform simulation + dt = [0.75] + power = 1.0 + openmc.deplete.si_leqi(op, dt, power, print_out=False) + + # Load the files + prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Re-create depletion operator and load previous results + op = dummy_operator.DummyOperator(prev_res) + op.output_dir = output_dir + + # Perform restarts simulation + openmc.deplete.si_leqi(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [2.03325094, 1.16826254] + s2 = [2.69291933, 0.37907772] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[3] == approx(s2[0]) + assert y2[3] == approx(s2[1]) From bfa9d75268e49c65546c7cc9ce84e47fd7e8e0de Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 18 Jan 2019 15:13:07 -0500 Subject: [PATCH 13/24] install mpi4py pakage in docker container --- Dockerfile | 1 + 1 file changed, 1 insertion(+) diff --git a/Dockerfile b/Dockerfile index 6466c7685..1bf82fb25 100644 --- a/Dockerfile +++ b/Dockerfile @@ -28,6 +28,7 @@ RUN git clone https://github.com/njoy/NJOY2016 /opt/NJOY2016 && \ cmake -Dstatic=on .. && make 2>/dev/null && make install # Clone and install OpenMC +RUN pip3 install --no-binary=mpi4py mpi4py RUN git clone https://github.com/openmc-dev/openmc.git /opt/openmc && \ cd /opt/openmc && mkdir -p build && cd build && \ cmake -Doptimize=on -DHDF5_PREFER_PARALLEL=on .. && \ From 56e6c4088de1892ed88c389d735d5fac8357c89c Mon Sep 17 00:00:00 2001 From: liangjg Date: Sat, 19 Jan 2019 09:40:04 -0500 Subject: [PATCH 14/24] fix: differentiate burnable mats before openmc.reset_auto_ids --- openmc/deplete/operator.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/deplete/operator.py b/openmc/deplete/operator.py index 754e22f7d..66cac392d 100644 --- a/openmc/deplete/operator.py +++ b/openmc/deplete/operator.py @@ -125,6 +125,10 @@ class Operator(TransportOperator): else: self.prev_res = None + # Differentiate burnable materials with multiple instances + if self.diff_burnable_mats: + self._differentiate_burnable_mats() + # Clear out OpenMC, create task lists, distribute openmc.reset_auto_ids() self.burnable_mats, volume, nuclides = self._get_burnable_mats() @@ -227,10 +231,6 @@ class Operator(TransportOperator): """ - if self.diff_burnable_mats: - # Automatically distribute burnable materials - self._differentiate_burnable_mats() - burnable_mats = set() model_nuclides = set() volume = OrderedDict() From 9b93c809eb47497b33fa1f2b7ff2b0c8de07328b Mon Sep 17 00:00:00 2001 From: liangjg Date: Sun, 20 Jan 2019 17:01:58 -0500 Subject: [PATCH 15/24] read WMP and write message in load_nuclide capi function --- src/nuclide_header.F90 | 43 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 75254cda1..290ceba70 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -11,6 +11,9 @@ module nuclide_header use endf_header, only: Function1D, Polynomial, Tabulated1D use error use hdf5_interface + use math, only: faddeeva, w_derivative, & + broaden_wmp_polynomials + use multipole_header use message_passing use random_lcg, only: prn, future_prn, prn_set_stream use reaction_header, only: Reaction @@ -728,6 +731,8 @@ contains integer(C_INT) :: err integer :: n + logical :: file_exists + character(7) :: readable integer(HID_T) :: file_id integer(HID_T) :: group_id character(:), allocatable :: name_ @@ -751,6 +756,9 @@ contains filename = library_path(LIBRARY_NEUTRON, to_lower(name_)) + call write_message('Reading ' // trim(name_) // ' from ' // & + trim(filename), 6) + ! Open file and make sure version is sufficient file_id = file_open(filename, 'r') call check_data_version(file_id) @@ -769,6 +777,41 @@ contains ! Initialize nuclide grid call nuclides(n) % init_grid() + + ! Read WMP data if necessary + nuclides(n) % mp_present = .false. + if (temperature_multipole .and. & + library_present(LIBRARY_WMP, to_lower(name_))) then + filename = library_path(LIBRARY_WMP, to_lower(name_)) + + ! Check if Multipole library exists and is readable + inquire(FILE=filename, EXIST=file_exists, READ=readable) + if (file_exists) then + if (readable(1:3) == 'NO') then + call fatal_error("Multipole library '" // trim(filename) // "' & + &is not readable! Change file permissions.") + end if + + ! Display message + call write_message("Reading " // trim(name_) // " WMP data from " & + // filename, 6) + + ! Open file and make sure version is sufficient + file_id = file_open(filename, 'r') + call check_wmp_version(file_id) + + ! Read nuclide data from HDF5 + group_id = open_group(file_id, name_) + allocate(nuclides(n) % multipole) + call nuclides(n) % multipole % from_hdf5(group_id) + + nuclides(n) % mp_present = .true. + + ! Close the group and file. + call close_group(group_id) + call file_close(file_id) + end if + end if else err = E_DATA call set_errmsg("Nuclide '" // trim(name_) // "' is not present & From 6f1da3077f3d2d2065862043fc42bf56160dd702 Mon Sep 17 00:00:00 2001 From: liangjg Date: Sun, 20 Jan 2019 21:56:34 -0500 Subject: [PATCH 16/24] defer write_summary calling in case the model is changed, for example materials are updated by capi in depletion after reading xml --- src/input_xml.F90 | 4 ---- src/simulation.F90 | 6 ++++++ 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index be93a164e..ef0835abd 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -29,7 +29,6 @@ module input_xml use string, only: to_lower, to_str, str_to_int, str_to_real, & starts_with, ends_with, split_string, & zero_padded, to_c_string - use summary, only: write_summary use tally use tally_header, only: openmc_extend_tallies use tally_derivative_header @@ -158,9 +157,6 @@ contains ! Normalize atom/weight percents call normalize_ao() - ! Write summary information - if (master .and. output_summary) call write_summary() - ! Warn if overlap checking is on if (master .and. check_overlaps) & call warning("Cell overlap checking is ON.") diff --git a/src/simulation.F90 b/src/simulation.F90 index f00b58ed1..c79b2188b 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -6,6 +6,9 @@ module simulation use nuclide_header, only: micro_xs, n_nuclides use photon_header, only: micro_photon_xs, n_elements use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer + use summary, only: write_summary + use settings, only: output_summary + use message_passing, only: master implicit none private @@ -20,6 +23,9 @@ contains integer :: i + ! Write summary information + if (master .and. output_summary) call write_summary() + ! Set up material nuclide index mapping do i = 1, n_materials call materials(i) % init_nuclide_index() From 02a977c5c12745c8ca1c87791eeb2260a000309b Mon Sep 17 00:00:00 2001 From: liangjg Date: Sun, 20 Jan 2019 23:46:23 -0500 Subject: [PATCH 17/24] rebase fix --- src/input_xml.F90 | 64 -------------------------- src/nuclide_header.F90 | 102 +++++++++++++++++++++++++---------------- 2 files changed, 63 insertions(+), 103 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ef0835abd..3c616edb6 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2237,68 +2237,4 @@ contains end subroutine read_ce_cross_sections -!=============================================================================== -! READ_MULTIPOLE_DATA checks for the existence of a multipole library in the -! directory and loads it using multipole_read -!=============================================================================== - - subroutine read_multipole_data(i_table) - - integer, intent(in) :: i_table ! index in nuclides/sab_tables - - logical :: file_exists ! Does multipole library exist? - character(7) :: readable ! Is multipole library readable? - character(MAX_FILE_LEN) :: filename ! Path to multipole xs library - integer(HID_T) :: file_id - integer(HID_T) :: group_id - - interface - subroutine nuclide_load_multipole(ptr, group) bind(C) - import C_PTR, HID_T - type(C_PTR), value :: ptr - integer(HID_T), value :: group - end subroutine - end interface - - associate (nuc => nuclides(i_table)) - - ! Look for WMP data in cross_sections.xml - if (library_present(LIBRARY_WMP, to_lower(nuc % name))) then - filename = library_path(LIBRARY_WMP, to_lower(nuc % name)) - else - nuc % mp_present = .false. - return - end if - - ! Check if Multipole library exists and is readable - inquire(FILE=filename, EXIST=file_exists, READ=readable) - if (.not. file_exists) then - nuc % mp_present = .false. - return - elseif (readable(1:3) == 'NO') then - call fatal_error("Multipole library '" // trim(filename) // "' is not & - &readable! Change file permissions with chmod command.") - end if - - ! Display message - call write_message("Reading " // trim(nuc % name) // " WMP data from " & - // filename, 6) - - ! Open file and make sure version is sufficient - file_id = file_open(filename, 'r') - call check_wmp_version(file_id) - - ! Read nuclide data from HDF5 - group_id = open_group(file_id, nuc % name) - nuc % mp_present = .true. - call nuclide_load_multipole(nuc % ptr, group_id) - call close_group(group_id) - - ! Close the file - call file_close(file_id) - - end associate - - end subroutine read_multipole_data - end module input_xml diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 290ceba70..04c2f02c0 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -11,10 +11,8 @@ module nuclide_header use endf_header, only: Function1D, Polynomial, Tabulated1D use error use hdf5_interface - use math, only: faddeeva, w_derivative, & - broaden_wmp_polynomials - use multipole_header use message_passing + use multipole_header use random_lcg, only: prn, future_prn, prn_set_stream use reaction_header, only: Reaction use settings @@ -731,8 +729,6 @@ contains integer(C_INT) :: err integer :: n - logical :: file_exists - character(7) :: readable integer(HID_T) :: file_id integer(HID_T) :: group_id character(:), allocatable :: name_ @@ -778,40 +774,8 @@ contains ! Initialize nuclide grid call nuclides(n) % init_grid() - ! Read WMP data if necessary - nuclides(n) % mp_present = .false. - if (temperature_multipole .and. & - library_present(LIBRARY_WMP, to_lower(name_))) then - filename = library_path(LIBRARY_WMP, to_lower(name_)) - - ! Check if Multipole library exists and is readable - inquire(FILE=filename, EXIST=file_exists, READ=readable) - if (file_exists) then - if (readable(1:3) == 'NO') then - call fatal_error("Multipole library '" // trim(filename) // "' & - &is not readable! Change file permissions.") - end if - - ! Display message - call write_message("Reading " // trim(name_) // " WMP data from " & - // filename, 6) - - ! Open file and make sure version is sufficient - file_id = file_open(filename, 'r') - call check_wmp_version(file_id) - - ! Read nuclide data from HDF5 - group_id = open_group(file_id, name_) - allocate(nuclides(n) % multipole) - call nuclides(n) % multipole % from_hdf5(group_id) - - nuclides(n) % mp_present = .true. - - ! Close the group and file. - call close_group(group_id) - call file_close(file_id) - end if - end if + ! Read multipole file into the appropriate entry on the nuclides array + if (temperature_multipole) call read_multipole_data(n) else err = E_DATA call set_errmsg("Nuclide '" // trim(name_) // "' is not present & @@ -846,4 +810,64 @@ contains end if end function openmc_nuclide_name + subroutine read_multipole_data(i_table) + + integer, intent(in) :: i_table ! index in nuclides/sab_tables + + logical :: file_exists ! Does multipole library exist? + character(7) :: readable ! Is multipole library readable? + character(MAX_FILE_LEN) :: filename ! Path to multipole xs library + integer(HID_T) :: file_id + integer(HID_T) :: group_id + + interface + subroutine nuclide_load_multipole(ptr, group) bind(C) + import C_PTR, HID_T + type(C_PTR), value :: ptr + integer(HID_T), value :: group + end subroutine + end interface + + associate (nuc => nuclides(i_table)) + + ! Look for WMP data in cross_sections.xml + if (library_present(LIBRARY_WMP, to_lower(nuc % name))) then + filename = library_path(LIBRARY_WMP, to_lower(nuc % name)) + else + nuc % mp_present = .false. + return + end if + + ! Check if Multipole library exists and is readable + inquire(FILE=filename, EXIST=file_exists, READ=readable) + if (.not. file_exists) then + nuc % mp_present = .false. + return + elseif (readable(1:3) == 'NO') then + call fatal_error("Multipole library '" // trim(filename) // "' is not & + &readable! Change file permissions with chmod command.") + end if + + ! Display message + call write_message("Reading " // trim(nuc % name) // " WMP data from " & + // filename, 6) + + ! Open file and make sure version is sufficient + file_id = file_open(filename, 'r') + call check_wmp_version(file_id) + + ! Read nuclide data from HDF5 + group_id = open_group(file_id, nuc % name) + nuc % mp_present = .true. + call nuclide_load_multipole(nuc % ptr, group_id) + call close_group(group_id) + + ! Close the file + call file_close(file_id) + + end associate + + end subroutine read_multipole_data + + end module nuclide_header From 9cc923faad0ed163232b78ed78244a8060aaa097 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 23 Jan 2019 14:59:43 -0500 Subject: [PATCH 18/24] moved CE/LI and LE/QI implementation from opendeplete --- openmc/deplete/integrator/celi.py | 176 +++++++++++++++++++++++++++++ openmc/deplete/integrator/leqi.py | 178 ++++++++++++++++++++++++++++++ 2 files changed, 354 insertions(+) create mode 100644 openmc/deplete/integrator/celi.py create mode 100644 openmc/deplete/integrator/leqi.py diff --git a/openmc/deplete/integrator/celi.py b/openmc/deplete/integrator/celi.py new file mode 100644 index 000000000..3a8b5286c --- /dev/null +++ b/openmc/deplete/integrator/celi.py @@ -0,0 +1,176 @@ +""" The CE/LI CFQ4 integrator. + +Implements the CE/LI Predictor-Corrector algorithm using commutator free +high order integrators. + +This algorithm is mathematically defined as: + +.. math: + y' = A(y, t) y(t) + A_p = A(y_n, t_n) + y_p = expm(A_p h) y_n + A_c = A(y_p, t_n) + A(t) = t/dt * A_c + (dt - t)/dt * A_p + +Here, A(t) is integrated using the fourth order algorithm described below. + +From +---- + Thalhammer, Mechthild. "A fourth-order commutator-free exponential + integrator for nonautonomous differential equations." SIAM journal on + numerical analysis 44.2 (2006): 851-864. +""" + +import copy +import os +import time + +from mpi4py import MPI + +from .cram import CRAM48 +from .save_results import save_results + +def celi_cfq4(operator, print_out=True): + """ Performs integration of an operator using the CE/LI CFQ4 algorithm. + + Parameters + ---------- + operator : Operator + The operator object to simulate on. + print_out : bool, optional + Whether or not to print out time. + """ + + # Save current directory + dir_home = os.getcwd() + + # Move to folder + os.makedirs(operator.settings.output_dir, exist_ok=True) + os.chdir(operator.settings.output_dir) + + # Generate initial conditions + vec = operator.initial_condition() + + t = 0.0 + + for i, dt in enumerate(operator.settings.dt_vec): + vec, t, _ = celi_cfq4_inner(operator, vec, i, t, dt, print_out) + + # Perform one last simulation + x = [copy.deepcopy(vec)] + seeds = [] + eigvls = [] + rates_array = [] + eigvl, rates, seed = operator.eval(x[0]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(rates) + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t], + len(operator.settings.dt_vec)) + + # Return to origin + os.chdir(dir_home) + +def celi_cfq4_inner(operator, vec, i, t, dt, print_out): + """ The inner loop of CE/LI CFQ4. + + Parameters + ---------- + operator : Operator + The operator object to simulate on. + vec : list of numpy.array + Nuclide vector, beginning of time. + i : Int + Current iteration number. + t : Float + Time at start of step. + dt : Float + Time step. + print_out : bool + Whether or not to print out time. + + Returns + ------- + x_result : list of numpy.array + Nuclide vector, end of time. + Float + Next time + ReactionRates + Reaction rates from beginning of step. + """ + + n_mats = len(vec) + + # Create vectors + x = [copy.deepcopy(vec)] + seeds = [] + eigvls = [] + rates_array = [] + + eigvl, rates, seed = operator.eval(x[0]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrix + f = operator.form_matrix(rates_array[0], mat) + + x_new = CRAM48(f, x[0][mat], dt) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + x.append(x_result) + + eigvl, rates, seed = operator.eval(x[1]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrices + f1 = dt * operator.form_matrix(rates_array[0], mat) + f2 = dt * operator.form_matrix(rates_array[1], mat) + + # Perform commutator-free integral + x_new = copy.deepcopy(x[0][mat]) + + # Compute linearly interpolated f at points + # A{1,2} = f(1/2 -/+ sqrt(3)/6) + # Then + # a{1,2} = 1/4 +/- sqrt(3)/6 + # m1 = a2 * A1 + a1 * A2 + # m2 = a1 * A1 + a2 * A2 + m1 = 1/12 * (f1 + 5 * f2) + m2 = 1/12 * (5 * f1 + f2) + + x_new = CRAM48(m2, x_new, 1.0) + x_new = CRAM48(m1, x_new, 1.0) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i) + + return x_result, t + dt, rates_array[0] \ No newline at end of file diff --git a/openmc/deplete/integrator/leqi.py b/openmc/deplete/integrator/leqi.py new file mode 100644 index 000000000..a79ca983c --- /dev/null +++ b/openmc/deplete/integrator/leqi.py @@ -0,0 +1,178 @@ +""" The LE/QI CFQ4 integrator. + +Implements the LE/QI Predictor-Corrector algorithm using commutator free +high order integrators. + +This algorithm is mathematically defined as: + +.. math: + y' = A(y, t) y(t) + A_m1 = A(y_n-1, t_n-1) + A_0 = A(y_n, t_n) + A_l(t) linear extrapolation of A_m1, A_0 + Integrate to t_n+1 to get y_p + A_c = A(y_p, y_n+1) + A_q(t) quadratic interpolation of A_m1, A_0, A_c + +Here, A(t) is integrated using the fourth order algorithm described below. + +From +---- + Thalhammer, Mechthild. "A fourth-order commutator-free exponential + integrator for nonautonomous differential equations." SIAM journal on + numerical analysis 44.2 (2006): 851-864. + +It is initialized using the CE/LI algorithm. +""" + +import copy +import os +import time + +from mpi4py import MPI + +from .celi_cfq4 import celi_cfq4_inner +from .cram import CRAM48 +from .save_results import save_results + +def leqi_cfq4(operator, print_out=True): + """ Performs integration of an operator using the LE/QI CFQ4 algorithm. + + Parameters + ---------- + operator : Operator + The operator object to simulate on. + print_out : bool, optional + Whether or not to print out time. + """ + + # Save current directory + dir_home = os.getcwd() + + # Move to folder + os.makedirs(operator.settings.output_dir, exist_ok=True) + os.chdir(operator.settings.output_dir) + + # Generate initial conditions + vec = operator.initial_condition() + + n_mats = len(vec) + + t = 0.0 + + # Perform single step of CE/LI CFQ4 + dt_l = operator.settings.dt_vec[0] + vec, t, rates_last = celi_cfq4_inner(operator, vec, 0, t, dt_l, print_out) + + # Perform remaining LE/QI + for i, dt in enumerate(operator.settings.dt_vec[1::]): + # Create vectors + x = [copy.deepcopy(vec)] + seeds = [] + eigvls = [] + rates_array = [] + + eigvl, rates, seed = operator.eval(x[0]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrices + f1 = dt * operator.form_matrix(rates_last, mat) + f2 = dt * operator.form_matrix(rates_array[0], mat) + + # Perform commutator-free integral + x_new = copy.deepcopy(x[0][mat]) + + # Compute linearly extrapolated f at points + # A{1,2} = f(1/2 -/+ sqrt(3)/6) + # Then + # a{1,2} = 1/4 +/- sqrt(3)/6 + # m1 = a2 * A1 + a1 * A2 + # m2 = a1 * A1 + a2 * A2 + m1 = -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 + m2 = -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 + + x_new = CRAM48(m2, x_new, 1.0) + x_new = CRAM48(m1, x_new, 1.0) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + x.append(x_result) + + eigvl, rates, seed = operator.eval(x[1]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(copy.deepcopy(rates)) + + x_result = [] + + t_start = time.time() + for mat in range(n_mats): + # Form matrices + f1 = dt * operator.form_matrix(rates_last, mat) + f2 = dt * operator.form_matrix(rates_array[0], mat) + f3 = dt * operator.form_matrix(rates_array[1], mat) + + # Perform commutator-free integral + x_new = copy.deepcopy(x[0][mat]) + + # Compute quadratically interpolated f at points + # A{1,2} = f(1/2 -/+ sqrt(3)/6) + # Then + # a{1,2} = 1/4 +/- sqrt(3)/6 + # m1 = a2 * A1 + a1 * A2 + # m2 = a1 * A1 + a2 * A2 + m1 = (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + + (dt**2 + 2 * dt * dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + + (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3) + m2 = (-dt**2/(12 * dt_l * (dt + dt_l)) * f1 + + (dt**2 + 6 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + + dt_l / (12 * (dt + dt_l)) * f3) + + x_new = CRAM48(m2, x_new, 1.0) + x_new = CRAM48(m1, x_new, 1.0) + + x_result.append(x_new) + + t_end = time.time() + if MPI.COMM_WORLD.rank == 0: + if print_out: + print("Time to matexp: ", t_end - t_start) + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i + 1) + + rates_last = copy.deepcopy(rates_array[0]) + t += dt + dt_l = dt + vec = copy.deepcopy(x_result) + + # Perform one last simulation + x = [copy.deepcopy(vec)] + seeds = [] + eigvls = [] + rates_array = [] + eigvl, rates, seed = operator.eval(x[0]) + + eigvls.append(eigvl) + seeds.append(seed) + rates_array.append(rates) + + # Create results, write to disk + save_results(operator, x, rates_array, eigvls, seeds, [t, t], + len(operator.settings.dt_vec)) + + # Return to origin + os.chdir(dir_home) From cbbac7df8d4c07c7a3a3e58da84c8ff813e9d999 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 23 Jan 2019 15:08:55 -0500 Subject: [PATCH 19/24] refactor CE/LI and LE/QI --- openmc/deplete/integrator/__init__.py | 2 + openmc/deplete/integrator/celi.py | 250 +++++++++++----------- openmc/deplete/integrator/leqi.py | 290 ++++++++++++-------------- openmc/deplete/integrator/si_celi.py | 10 +- openmc/deplete/integrator/si_leqi.py | 32 +-- openmc/model/model.py | 2 +- 6 files changed, 261 insertions(+), 325 deletions(-) diff --git a/openmc/deplete/integrator/__init__.py b/openmc/deplete/integrator/__init__.py index 4dad78e1c..db906b3c1 100644 --- a/openmc/deplete/integrator/__init__.py +++ b/openmc/deplete/integrator/__init__.py @@ -7,8 +7,10 @@ The integrator subcomponents. from .cf4 import * from .cecm import * +from .celi import * from .cram import * from .epc_rk4 import * +from .leqi import * from .predictor import * from .si_celi import * from .si_leqi import * diff --git a/openmc/deplete/integrator/celi.py b/openmc/deplete/integrator/celi.py index 3a8b5286c..61edaeb72 100644 --- a/openmc/deplete/integrator/celi.py +++ b/openmc/deplete/integrator/celi.py @@ -1,176 +1,170 @@ -""" The CE/LI CFQ4 integrator. - -Implements the CE/LI Predictor-Corrector algorithm using commutator free -high order integrators. - -This algorithm is mathematically defined as: - -.. math: - y' = A(y, t) y(t) - A_p = A(y_n, t_n) - y_p = expm(A_p h) y_n - A_c = A(y_p, t_n) - A(t) = t/dt * A_c + (dt - t)/dt * A_p - -Here, A(t) is integrated using the fourth order algorithm described below. - -From ----- - Thalhammer, Mechthild. "A fourth-order commutator-free exponential - integrator for nonautonomous differential equations." SIAM journal on - numerical analysis 44.2 (2006): 851-864. -""" +"""The CE/LI CFQ4 integrator.""" import copy -import os -import time +from collections.abc import Iterable -from mpi4py import MPI +from .cram import deplete +from ..results import Results +from ..abc import OperatorResult -from .cram import CRAM48 -from .save_results import save_results -def celi_cfq4(operator, print_out=True): - """ Performs integration of an operator using the CE/LI CFQ4 algorithm. +# Functions to form the special matrix for depletion +def _celi_f1(chain, rates): + return 5/12 * chain.form_matrix(rates[0]) + \ + 1/12 * chain.form_matrix(rates[1]) + +def _celi_f2(chain, rates): + return 1/12 * chain.form_matrix(rates[0]) + \ + 5/12 * chain.form_matrix(rates[1]) + +def celi(operator, timesteps, power=None, power_density=None, + print_out=True): + r"""Deplete using the CE/LI CFQ4 algorithm. + + Implements the CE/LI Predictor-Corrector algorithm using the [fourth order + commutator-free integrator]_. + + The CE/LI algorithm is mathematically defined as: + + .. math: + y' = A(y, t) y(t) + A_p = A(y_n, t_n) + y_p = expm(A_p h) y_n + A_c = A(y_p, t_n) + A(t) = t/dt * A_c + (dt - t)/dt * A_p + + Here, A(t) is integrated using the fourth order algorithm CFQ4. Parameters ---------- - operator : Operator + 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. + + References + ---------- + .. [fourth order commutator-free integrator] + Thalhammer, Mechthild. "A fourth-order commutator-free exponential + integrator for nonautonomous differential equations." SIAM journal on + numerical analysis 44.2 (2006): 851-864. """ + 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] - # Save current directory - dir_home = os.getcwd() - - # Move to folder - os.makedirs(operator.settings.output_dir, exist_ok=True) - os.chdir(operator.settings.output_dir) + if not isinstance(power, Iterable): + power = [power]*len(timesteps) # Generate initial conditions - vec = operator.initial_condition() + 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) - t = 0.0 + for i, (dt, p) in enumerate(zip(timesteps, power)): + vec, t, _ = celi_inner(operator, vec, p, i, i_res, t, dt, + print_out) - for i, dt in enumerate(operator.settings.dt_vec): - vec, t, _ = celi_cfq4_inner(operator, vec, i, t, dt, print_out) + # Perform one last simulation + x = [copy.deepcopy(vec)] + op_results = [operator(x[0], power[-1])] - # Perform one last simulation - x = [copy.deepcopy(vec)] - seeds = [] - eigvls = [] - rates_array = [] - eigvl, rates, seed = operator.eval(x[0]) + # Create results, write to disk + Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps)) - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(rates) - - # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t], - len(operator.settings.dt_vec)) - - # Return to origin - os.chdir(dir_home) - -def celi_cfq4_inner(operator, vec, i, t, dt, print_out): +def celi_inner(operator, vec, p, i, i_res, t, dt, print_out): """ The inner loop of CE/LI CFQ4. Parameters ---------- operator : Operator The operator object to simulate on. - vec : list of numpy.array + x : list of nuclide vector Nuclide vector, beginning of time. - i : Int + p : float + Power of the reactor in [W] + i : int Current iteration number. - t : Float + i_res : int + Starting index, for restart calculation. + t : float Time at start of step. - dt : Float + dt : float Time step. print_out : bool Whether or not to print out time. Returns ------- - x_result : list of numpy.array + list of numpy.array Nuclide vector, end of time. - Float + float Next time - ReactionRates - Reaction rates from beginning of step. + OperatorResult + Operator result from beginning of step. """ - n_mats = len(vec) + chain = operator.chain - # Create vectors - x = [copy.deepcopy(vec)] - seeds = [] - eigvls = [] - rates_array = [] + # 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)] - eigvl, rates, seed = operator.eval(x[0]) + else: + # Get initial concentration + x = [operator.prev_res[-1].data[0]] - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) + # Get rates + op_results = [operator.prev_res[-1]] + op_results[0].rates = op_results[0].rates[0] - x_result = [] + # Set first stage value of keff + op_results[0].k = op_results[0].k[0] - t_start = time.time() - for mat in range(n_mats): - # Form matrix - f = operator.form_matrix(rates_array[0], mat) + # 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] - x_new = CRAM48(f, x[0][mat], dt) + # Deplete to end + x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out) + x.append(x_new) + op_results.append(operator(x[1], p)) - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - x.append(x_result) - - eigvl, rates, seed = operator.eval(x[1]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrices - f1 = dt * operator.form_matrix(rates_array[0], mat) - f2 = dt * operator.form_matrix(rates_array[1], mat) - - # Perform commutator-free integral - x_new = copy.deepcopy(x[0][mat]) - - # Compute linearly interpolated f at points - # A{1,2} = f(1/2 -/+ sqrt(3)/6) - # Then - # a{1,2} = 1/4 +/- sqrt(3)/6 - # m1 = a2 * A1 + a1 * A2 - # m2 = a1 * A1 + a2 * A2 - m1 = 1/12 * (f1 + 5 * f2) - m2 = 1/12 * (5 * f1 + f2) - - x_new = CRAM48(m2, x_new, 1.0) - x_new = CRAM48(m1, x_new, 1.0) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) + # Deplete with two matrix exponentials + rates = list(zip(op_results[0].rates, op_results[1].rates)) + x_end = deplete(chain, x[0], rates, dt, print_out, + matrix_func=_celi_f1) + x_end = deplete(chain, x_end, rates, dt, print_out, + matrix_func=_celi_f2) # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i) + Results.save(operator, x, op_results, [t, t + dt], p, i_res + i) - return x_result, t + dt, rates_array[0] \ No newline at end of file + # return updated time and vectors + return x_end, t + dt, op_results[0] diff --git a/openmc/deplete/integrator/leqi.py b/openmc/deplete/integrator/leqi.py index a79ca983c..77875c40b 100644 --- a/openmc/deplete/integrator/leqi.py +++ b/openmc/deplete/integrator/leqi.py @@ -1,178 +1,156 @@ -""" The LE/QI CFQ4 integrator. - -Implements the LE/QI Predictor-Corrector algorithm using commutator free -high order integrators. - -This algorithm is mathematically defined as: - -.. math: - y' = A(y, t) y(t) - A_m1 = A(y_n-1, t_n-1) - A_0 = A(y_n, t_n) - A_l(t) linear extrapolation of A_m1, A_0 - Integrate to t_n+1 to get y_p - A_c = A(y_p, y_n+1) - A_q(t) quadratic interpolation of A_m1, A_0, A_c - -Here, A(t) is integrated using the fourth order algorithm described below. - -From ----- - Thalhammer, Mechthild. "A fourth-order commutator-free exponential - integrator for nonautonomous differential equations." SIAM journal on - numerical analysis 44.2 (2006): 851-864. - -It is initialized using the CE/LI algorithm. -""" +"""The LE/QI CFQ4 integrator.""" import copy -import os -import time +from collections.abc import Iterable +from itertools import repeat -from mpi4py import MPI +from .celi import celi_inner +from .cram import deplete +from ..results import Results +from ..abc import OperatorResult -from .celi_cfq4 import celi_cfq4_inner -from .cram import CRAM48 -from .save_results import save_results -def leqi_cfq4(operator, print_out=True): - """ Performs integration of an operator using the LE/QI CFQ4 algorithm. +# Functions to form the special matrix for depletion +def _leqi_f1(chain, inputs): + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + dt_l, dt = inputs[2], inputs[3] + return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 + +def _leqi_f2(chain, inputs): + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + dt_l, dt = inputs[2], inputs[3] + return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 + +def _leqi_f3(chain, inputs): + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + f3 = chain.form_matrix(inputs[2]) + dt_l, dt = inputs[3], inputs[4] + return (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + + (dt**2 + 6*dt*dt_l + 5*dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + + dt_l / (12 * (dt + dt_l)) * f3) + +def _leqi_f4(chain, inputs): + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + f3 = chain.form_matrix(inputs[2]) + dt_l, dt = inputs[3], inputs[4] + return (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + + (dt**2 + 2*dt*dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + + (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3) + +def leqi(operator, timesteps, power=None, power_density=None, print_out=True): + r"""Deplete using the LE/QI CFQ4 algorithm. + + Implements the LE/QI Predictor-Corrector algorithm using the [fourth order + commutator-free integrator]_. + + The LE/QI algorithm is mathematically defined as: + + .. math: + y' = A(y, t) y(t) + A_m1 = A(y_n-1, t_n-1) + A_0 = A(y_n, t_n) + A_l(t) linear extrapolation of A_m1, A_0 + Integrate to t_n+1 to get y_p + A_c = A(y_p, y_n+1) + A_q(t) quadratic interpolation of A_m1, A_0, A_c + + Here, A(t) is integrated using the fourth order algorithm CFQ4. + + It is initialized using the CE/LI algorithm. Parameters ---------- - operator : Operator + 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. + + References + ---------- + .. [fourth order commutator-free integrator] + Thalhammer, Mechthild. "A fourth-order commutator-free exponential + integrator for nonautonomous differential equations." SIAM journal on + numerical analysis 44.2 (2006): 851-864. """ + 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] - # Save current directory - dir_home = os.getcwd() - - # Move to folder - os.makedirs(operator.settings.output_dir, exist_ok=True) - os.chdir(operator.settings.output_dir) + if not isinstance(power, Iterable): + power = [power]*len(timesteps) # Generate initial conditions - vec = operator.initial_condition() + 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) - n_mats = len(vec) + for i, (dt, p) in enumerate(zip(timesteps, power)): + # Perform SI-CE/LI CFQ4 for the first step + if i == 0: + # Save results for the last step + dt_l = dt + vec, t, op_res_last = celi_inner(operator, vec, p, i, i_res, + t, dt, print_out) + continue - t = 0.0 + # Perform remaining LE/QI + x = [copy.deepcopy(vec)] + op_results = [operator(x[0], p)] - # Perform single step of CE/LI CFQ4 - dt_l = operator.settings.dt_vec[0] - vec, t, rates_last = celi_cfq4_inner(operator, vec, 0, t, dt_l, print_out) + inputs = list(zip(op_res_last.rates, op_results[0].rates, + repeat(dt_l), repeat(dt))) + x_new = deplete(chain, x[0], inputs, dt, print_out, + matrix_func=_leqi_f1) + x_new = deplete(chain, x_new, inputs, dt, print_out, + matrix_func=_leqi_f2) + x.append(x_new) + op_results.append(operator(x[1], p)) - # Perform remaining LE/QI - for i, dt in enumerate(operator.settings.dt_vec[1::]): - # Create vectors + inputs = list(zip(op_res_last.rates, op_results[0].rates, + op_results[1].rates, repeat(dt_l), repeat(dt))) + x_new = deplete(chain, x[0], inputs, dt, print_out, + matrix_func=_leqi_f3) + x_new = deplete(chain, x_new, inputs, dt, print_out, + matrix_func=_leqi_f4) + + # Create results, write to disk + Results.save(operator, x, op_results, [t, t+dt], p, i_res+i) + + # update results + x = [x_new] + op_res_last = copy.deepcopy(op_results[0]) + t += dt + dt_l = dt + + # Perform one last simulation x = [copy.deepcopy(vec)] - seeds = [] - eigvls = [] - rates_array = [] - - eigvl, rates, seed = operator.eval(x[0]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrices - f1 = dt * operator.form_matrix(rates_last, mat) - f2 = dt * operator.form_matrix(rates_array[0], mat) - - # Perform commutator-free integral - x_new = copy.deepcopy(x[0][mat]) - - # Compute linearly extrapolated f at points - # A{1,2} = f(1/2 -/+ sqrt(3)/6) - # Then - # a{1,2} = 1/4 +/- sqrt(3)/6 - # m1 = a2 * A1 + a1 * A2 - # m2 = a1 * A1 + a2 * A2 - m1 = -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 - m2 = -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 - - x_new = CRAM48(m2, x_new, 1.0) - x_new = CRAM48(m1, x_new, 1.0) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) - - x.append(x_result) - - eigvl, rates, seed = operator.eval(x[1]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(copy.deepcopy(rates)) - - x_result = [] - - t_start = time.time() - for mat in range(n_mats): - # Form matrices - f1 = dt * operator.form_matrix(rates_last, mat) - f2 = dt * operator.form_matrix(rates_array[0], mat) - f3 = dt * operator.form_matrix(rates_array[1], mat) - - # Perform commutator-free integral - x_new = copy.deepcopy(x[0][mat]) - - # Compute quadratically interpolated f at points - # A{1,2} = f(1/2 -/+ sqrt(3)/6) - # Then - # a{1,2} = 1/4 +/- sqrt(3)/6 - # m1 = a2 * A1 + a1 * A2 - # m2 = a1 * A1 + a2 * A2 - m1 = (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + - (dt**2 + 2 * dt * dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + - (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3) - m2 = (-dt**2/(12 * dt_l * (dt + dt_l)) * f1 + - (dt**2 + 6 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + - dt_l / (12 * (dt + dt_l)) * f3) - - x_new = CRAM48(m2, x_new, 1.0) - x_new = CRAM48(m1, x_new, 1.0) - - x_result.append(x_new) - - t_end = time.time() - if MPI.COMM_WORLD.rank == 0: - if print_out: - print("Time to matexp: ", t_end - t_start) + op_results = [operator(x[0], power[-1])] # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t + dt], i + 1) - - rates_last = copy.deepcopy(rates_array[0]) - t += dt - dt_l = dt - vec = copy.deepcopy(x_result) - - # Perform one last simulation - x = [copy.deepcopy(vec)] - seeds = [] - eigvls = [] - rates_array = [] - eigvl, rates, seed = operator.eval(x[0]) - - eigvls.append(eigvl) - seeds.append(seed) - rates_array.append(rates) - - # Create results, write to disk - save_results(operator, x, rates_array, eigvls, seeds, [t, t], - len(operator.settings.dt_vec)) - - # Return to origin - os.chdir(dir_home) + Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps)) diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py index dc9bccb70..96aa7e2ba 100644 --- a/openmc/deplete/integrator/si_celi.py +++ b/openmc/deplete/integrator/si_celi.py @@ -6,17 +6,9 @@ from collections.abc import Iterable from .cram import deplete from ..results import Results from ..abc import OperatorResult +from .celi import _celi_f1, _celi_f2 -# Functions to form the special matrix for depletion -def _celi_f1(chain, rates): - return 5/12 * chain.form_matrix(rates[0]) + \ - 1/12 * chain.form_matrix(rates[1]) - -def _celi_f2(chain, rates): - return 1/12 * chain.form_matrix(rates[0]) + \ - 5/12 * chain.form_matrix(rates[1]) - def si_celi(operator, timesteps, power=None, power_density=None, print_out=True, m=10): r"""Deplete using the SI-CE/LI CFQ4 algorithm. diff --git a/openmc/deplete/integrator/si_leqi.py b/openmc/deplete/integrator/si_leqi.py index 609b3f8c4..4433a2768 100644 --- a/openmc/deplete/integrator/si_leqi.py +++ b/openmc/deplete/integrator/si_leqi.py @@ -5,42 +5,12 @@ from collections.abc import Iterable from itertools import repeat from .si_celi import si_celi_inner +from .leqi import _leqi_f1, _leqi_f2, _leqi_f3, _leqi_f4 from .cram import deplete from ..results import Results from ..abc import OperatorResult -# Functions to form the special matrix for depletion -def _leqi_f1(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - dt_l, dt = inputs[2], inputs[3] - return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 - -def _leqi_f2(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - dt_l, dt = inputs[2], inputs[3] - return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 - -def _leqi_f3(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - f3 = chain.form_matrix(inputs[2]) - dt_l, dt = inputs[3], inputs[4] - return (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + - (dt**2 + 6*dt*dt_l + 5*dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + - dt_l / (12 * (dt + dt_l)) * f3) - -def _leqi_f4(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - f3 = chain.form_matrix(inputs[2]) - dt_l, dt = inputs[3], inputs[4] - return (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + - (dt**2 + 2*dt*dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + - (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3) - def si_leqi(operator, timesteps, power=None, power_density=None, print_out=True, m=10): r"""Deplete using the SI-LE/QI CFQ4 algorithm. diff --git a/openmc/model/model.py b/openmc/model/model.py index b98133208..dfa3decd3 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -151,7 +151,7 @@ class Model(object): # Perform depletion check_value('method', method, ('cecm', 'predictor', 'cf4', 'epc_rk4', - 'si_celi', 'si_leqi')) + 'si_celi', 'si_leqi', 'celi', 'leqi')) getattr(dep.integrator, method)(op, timesteps, **kwargs) def export_to_xml(self): From 0f639937b71a89d2dcf9a6c19cfdaa3b6132fb89 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 23 Jan 2019 15:50:01 -0500 Subject: [PATCH 20/24] added unit tests for CE/LI and LE/QI --- openmc/deplete/integrator/leqi.py | 4 +- tests/unit_tests/test_deplete_celi.py | 37 ++++++++++++ tests/unit_tests/test_deplete_leqi.py | 37 ++++++++++++ tests/unit_tests/test_deplete_restart.py | 77 ++++++++++++++++++++++++ 4 files changed, 154 insertions(+), 1 deletion(-) create mode 100644 tests/unit_tests/test_deplete_celi.py create mode 100644 tests/unit_tests/test_deplete_leqi.py diff --git a/openmc/deplete/integrator/leqi.py b/openmc/deplete/integrator/leqi.py index 77875c40b..3c6d26367 100644 --- a/openmc/deplete/integrator/leqi.py +++ b/openmc/deplete/integrator/leqi.py @@ -110,6 +110,8 @@ def leqi(operator, timesteps, power=None, power_density=None, print_out=True): t = operator.prev_res[-1].time[-1] i_res = len(operator.prev_res) + chain = operator.chain + for i, (dt, p) in enumerate(zip(timesteps, power)): # Perform SI-CE/LI CFQ4 for the first step if i == 0: @@ -149,7 +151,7 @@ def leqi(operator, timesteps, power=None, power_density=None, print_out=True): dt_l = dt # Perform one last simulation - x = [copy.deepcopy(vec)] + x = [copy.deepcopy(x_new)] op_results = [operator(x[0], power[-1])] # Create results, write to disk diff --git a/tests/unit_tests/test_deplete_celi.py b/tests/unit_tests/test_deplete_celi.py new file mode 100644 index 000000000..0b23b1a43 --- /dev/null +++ b/tests/unit_tests/test_deplete_celi.py @@ -0,0 +1,37 @@ +"""Regression tests for openmc.deplete.integrator.celi algorithm. + +These tests integrate a simple test problem described in dummy_geometry.py. +""" + +from pytest import approx +import openmc.deplete + +from tests import dummy_operator + + +def test_celi(run_in_tmpdir): + """Integral regression test of integrator algorithm using celi""" + + op = dummy_operator.DummyOperator() + op.output_dir = "test_integrator_regression" + + # Perform simulation using the celi algorithm + dt = [0.75, 0.75] + power = 1.0 + openmc.deplete.celi(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [1.82078767, 0.97122898] + s2 = [2.68441779, 0.05125966] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) diff --git a/tests/unit_tests/test_deplete_leqi.py b/tests/unit_tests/test_deplete_leqi.py new file mode 100644 index 000000000..a1a32e092 --- /dev/null +++ b/tests/unit_tests/test_deplete_leqi.py @@ -0,0 +1,37 @@ +"""Regression tests for openmc.deplete.integrator.leqi algorithm. + +These tests integrate a simple test problem described in dummy_geometry.py. +""" + +from pytest import approx +import openmc.deplete + +from tests import dummy_operator + + +def test_leqi(run_in_tmpdir): + """Integral regression test of integrator algorithm using leqi""" + + op = dummy_operator.DummyOperator() + op.output_dir = "test_integrator_regression" + + # Perform simulation using the leqi algorithm + dt = [0.75, 0.75] + power = 1.0 + openmc.deplete.leqi(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [1.82078767, 0.97122898] + s2 = [2.74526197, 0.23339915] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[2] == approx(s2[0]) + assert y2[2] == approx(s2[1]) diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py index 9c7e6c6c6..1e3a69980 100644 --- a/tests/unit_tests/test_deplete_restart.py +++ b/tests/unit_tests/test_deplete_restart.py @@ -246,6 +246,83 @@ def test_restart_epc_rk4(run_in_tmpdir): assert y2[3] == approx(s2[1]) +def test_restart_celi(run_in_tmpdir): + """Integral regression test of integrator algorithm using CELI.""" + + op = dummy_operator.DummyOperator() + output_dir = "test_restart_celi" + op.output_dir = output_dir + + # Perform simulation + dt = [0.75] + power = 1.0 + openmc.deplete.celi(op, dt, power, print_out=False) + + # Load the files + prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Re-create depletion operator and load previous results + op = dummy_operator.DummyOperator(prev_res) + op.output_dir = output_dir + + # Perform restarts simulation + openmc.deplete.celi(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [1.82078767, 0.97122898] + s2 = [2.68441779, 0.05125966] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[3] == approx(s2[0]) + assert y2[3] == approx(s2[1]) + + +def test_restart_leqi(run_in_tmpdir): + """Integral regression test of integrator algorithm using LEQI.""" + + op = dummy_operator.DummyOperator() + output_dir = "test_restart_leqi" + op.output_dir = output_dir + + # Perform simulation + dt = [0.75] + power = 1.0 + openmc.deplete.leqi(op, dt, power, print_out=False) + + # Load the files + prev_res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + # Re-create depletion operator and load previous results + op = dummy_operator.DummyOperator(prev_res) + op.output_dir = output_dir + + # Perform restarts simulation + openmc.deplete.leqi(op, dt, power, print_out=False) + + # Load the files + res = openmc.deplete.ResultsList(op.output_dir / "depletion_results.h5") + + _, y1 = res.get_atoms("1", "1") + _, y2 = res.get_atoms("1", "2") + + # Reference solution + s1 = [1.82078767, 0.97122898] + s2 = [2.74526197, 0.23339915] + + assert y1[1] == approx(s1[0]) + assert y2[1] == approx(s1[1]) + + assert y1[3] == approx(s2[0]) + assert y2[3] == approx(s2[1]) + def test_restart_si_celi(run_in_tmpdir): """Integral regression test of integrator algorithm using SI-CELI.""" From 6f3aeed89a76a6a128c5e5d8fd207ea00c2b0863 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 23 Jan 2019 16:53:31 -0500 Subject: [PATCH 21/24] fix test error --- docs/source/pythonapi/deplete.rst | 2 ++ openmc/deplete/integrator/leqi.py | 5 ++--- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index e38c00e96..7b24326d3 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -16,6 +16,8 @@ algorithms for depletion calculations. integrator.predictor integrator.cecm + integrator.celi + integrator.leqi integrator.cf4 integrator.epc_rk4 integrator.si_celi diff --git a/openmc/deplete/integrator/leqi.py b/openmc/deplete/integrator/leqi.py index 3c6d26367..04415dc21 100644 --- a/openmc/deplete/integrator/leqi.py +++ b/openmc/deplete/integrator/leqi.py @@ -117,12 +117,12 @@ def leqi(operator, timesteps, power=None, power_density=None, print_out=True): if i == 0: # Save results for the last step dt_l = dt - vec, t, op_res_last = celi_inner(operator, vec, p, i, i_res, + x_new, t, op_res_last = celi_inner(operator, vec, p, i, i_res, t, dt, print_out) continue # Perform remaining LE/QI - x = [copy.deepcopy(vec)] + x = [copy.deepcopy(x_new)] op_results = [operator(x[0], p)] inputs = list(zip(op_res_last.rates, op_results[0].rates, @@ -145,7 +145,6 @@ def leqi(operator, timesteps, power=None, power_density=None, print_out=True): Results.save(operator, x, op_results, [t, t+dt], p, i_res+i) # update results - x = [x_new] op_res_last = copy.deepcopy(op_results[0]) t += dt dt_l = dt From d65bdaac49a364cb3d99328b4fba540a718b4327 Mon Sep 17 00:00:00 2001 From: liangjg Date: Wed, 23 Jan 2019 21:28:53 -0500 Subject: [PATCH 22/24] fix restart bug for LE/QI methods --- Dockerfile | 1 - openmc/deplete/integrator/leqi.py | 18 ++++++++++++------ openmc/deplete/integrator/si_leqi.py | 18 ++++++++++++------ tests/unit_tests/test_deplete_restart.py | 2 +- 4 files changed, 25 insertions(+), 14 deletions(-) diff --git a/Dockerfile b/Dockerfile index 1bf82fb25..6466c7685 100644 --- a/Dockerfile +++ b/Dockerfile @@ -28,7 +28,6 @@ RUN git clone https://github.com/njoy/NJOY2016 /opt/NJOY2016 && \ cmake -Dstatic=on .. && make 2>/dev/null && make install # Clone and install OpenMC -RUN pip3 install --no-binary=mpi4py mpi4py RUN git clone https://github.com/openmc-dev/openmc.git /opt/openmc && \ cd /opt/openmc && mkdir -p build && cd build && \ cmake -Doptimize=on -DHDF5_PREFER_PARALLEL=on .. && \ diff --git a/openmc/deplete/integrator/leqi.py b/openmc/deplete/integrator/leqi.py index 04415dc21..c247bbe63 100644 --- a/openmc/deplete/integrator/leqi.py +++ b/openmc/deplete/integrator/leqi.py @@ -113,13 +113,19 @@ def leqi(operator, timesteps, power=None, power_density=None, print_out=True): chain = operator.chain for i, (dt, p) in enumerate(zip(timesteps, power)): - # Perform SI-CE/LI CFQ4 for the first step + # LE/QI needs the last step results to start + # Perform CE/LI CFQ4 or restore results for the first step if i == 0: - # Save results for the last step - dt_l = dt - x_new, t, op_res_last = celi_inner(operator, vec, p, i, i_res, - t, dt, print_out) - continue + if i_res <= 1: + dt_l = dt + x_new, t, op_res_last = celi_inner(operator, vec, 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_new = operator.prev_res[-1].data[0] # Perform remaining LE/QI x = [copy.deepcopy(x_new)] diff --git a/openmc/deplete/integrator/si_leqi.py b/openmc/deplete/integrator/si_leqi.py index 4433a2768..71c603ed3 100644 --- a/openmc/deplete/integrator/si_leqi.py +++ b/openmc/deplete/integrator/si_leqi.py @@ -113,14 +113,20 @@ def si_leqi(operator, timesteps, power=None, power_density=None, chain = operator.chain for i, (dt, p) in enumerate(zip(timesteps, power)): - # Perform SI-CE/LI CFQ4 for the first step + # LE/QI needs the last step results to start + # Perform SI-CE/LI CFQ4 or restore results for the first step if i == 0: - # Save results for the last step - op_res_last = copy.deepcopy(op_results[0]) dt_l = dt - x, t, op_results = si_celi_inner(operator, x, op_results, p, - i, i_res, t, dt, print_out) - continue + 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, diff --git a/tests/unit_tests/test_deplete_restart.py b/tests/unit_tests/test_deplete_restart.py index 1e3a69980..ec62064a4 100644 --- a/tests/unit_tests/test_deplete_restart.py +++ b/tests/unit_tests/test_deplete_restart.py @@ -392,7 +392,7 @@ def test_restart_si_leqi(run_in_tmpdir): # Reference solution s1 = [2.03325094, 1.16826254] - s2 = [2.69291933, 0.37907772] + s2 = [2.92711288, 0.53753236] assert y1[1] == approx(s1[0]) assert y2[1] == approx(s1[1]) From 2f3a1d5361c61ebd7eb487d854e60a2f7a04e541 Mon Sep 17 00:00:00 2001 From: liangjg Date: Fri, 25 Jan 2019 11:32:23 -0500 Subject: [PATCH 23/24] move write_summary back to input_xml --- src/input_xml.F90 | 4 ++++ src/simulation.F90 | 6 ------ 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a75c48e6f..72c14a645 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -29,6 +29,7 @@ module input_xml use string, only: to_lower, to_str, str_to_int, str_to_real, & starts_with, ends_with, split_string, & zero_padded, to_c_string + use summary, only: write_summary use tally use tally_header, only: openmc_extend_tallies use tally_derivative_header @@ -157,6 +158,9 @@ contains ! Normalize atom/weight percents call normalize_ao() + ! Write summary information + if (master .and. output_summary) call write_summary() + ! Warn if overlap checking is on if (master .and. check_overlaps) & call warning("Cell overlap checking is ON.") diff --git a/src/simulation.F90 b/src/simulation.F90 index c79b2188b..f00b58ed1 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -6,9 +6,6 @@ module simulation use nuclide_header, only: micro_xs, n_nuclides use photon_header, only: micro_photon_xs, n_elements use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer - use summary, only: write_summary - use settings, only: output_summary - use message_passing, only: master implicit none private @@ -23,9 +20,6 @@ contains integer :: i - ! Write summary information - if (master .and. output_summary) call write_summary() - ! Set up material nuclide index mapping do i = 1, n_materials call materials(i) % init_nuclide_index() From 760baa43d6840d8067032393ea09fdf45a5db6fe Mon Sep 17 00:00:00 2001 From: liangjg Date: Thu, 31 Jan 2019 12:01:55 -0500 Subject: [PATCH 24/24] address review comments from @paulromano --- docs/source/pythonapi/deplete.rst | 4 +- openmc/deplete/integrator/cecm.py | 6 +- openmc/deplete/integrator/celi.py | 40 ++++++------ openmc/deplete/integrator/cf4.py | 40 ++++++------ openmc/deplete/integrator/epc_rk4.py | 20 +++--- openmc/deplete/integrator/leqi.py | 91 +++++++++++++++------------- openmc/deplete/integrator/si_celi.py | 22 ++----- openmc/deplete/integrator/si_leqi.py | 25 +------- openmc/model/model.py | 2 +- 9 files changed, 106 insertions(+), 144 deletions(-) diff --git a/docs/source/pythonapi/deplete.rst b/docs/source/pythonapi/deplete.rst index 7b24326d3..532e06655 100644 --- a/docs/source/pythonapi/deplete.rst +++ b/docs/source/pythonapi/deplete.rst @@ -7,7 +7,9 @@ .. module:: openmc.deplete Several functions are provided that implement different time-integration -algorithms for depletion calculations. +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 `_. .. autosummary:: :toctree: generated diff --git a/openmc/deplete/integrator/cecm.py b/openmc/deplete/integrator/cecm.py index 926e27257..42024863a 100644 --- a/openmc/deplete/integrator/cecm.py +++ b/openmc/deplete/integrator/cecm.py @@ -11,8 +11,10 @@ def cecm(operator, timesteps, power=None, power_density=None, print_out=True): r"""Deplete using the CE/CM algorithm. Implements the second order `CE/CM predictor-corrector algorithm - `_. This algorithm is mathematically - defined as: + `_. + + "CE/CM" stands for constant extrapolation on predictor and constant + midpoint on corrector. This algorithm is mathematically defined as: .. math:: y' &= A(y, t) y(t) diff --git a/openmc/deplete/integrator/celi.py b/openmc/deplete/integrator/celi.py index 61edaeb72..0ba86d585 100644 --- a/openmc/deplete/integrator/celi.py +++ b/openmc/deplete/integrator/celi.py @@ -5,35 +5,38 @@ from collections.abc import Iterable from .cram import deplete from ..results import Results -from ..abc import OperatorResult # Functions to form the special matrix for depletion def _celi_f1(chain, rates): - return 5/12 * chain.form_matrix(rates[0]) + \ - 1/12 * chain.form_matrix(rates[1]) + return 5/12 * chain.form_matrix(rates[0]) + \ + 1/12 * chain.form_matrix(rates[1]) def _celi_f2(chain, rates): - return 1/12 * chain.form_matrix(rates[0]) + \ - 5/12 * chain.form_matrix(rates[1]) + return 1/12 * chain.form_matrix(rates[0]) + \ + 5/12 * chain.form_matrix(rates[1]) def celi(operator, timesteps, power=None, power_density=None, print_out=True): r"""Deplete using the CE/LI CFQ4 algorithm. - Implements the CE/LI Predictor-Corrector algorithm using the [fourth order - commutator-free integrator]_. + Implements the CE/LI Predictor-Corrector algorithm using the `fourth order + commutator-free integrator `_. - The CE/LI algorithm is mathematically defined as: + "CE/LI" stands for constant extrapolation on predictor and linear + interpolation on corrector. This algorithm is mathematically defined as: - .. math: - y' = A(y, t) y(t) - A_p = A(y_n, t_n) - y_p = expm(A_p h) y_n - A_c = A(y_p, t_n) - A(t) = t/dt * A_c + (dt - t)/dt * A_p + .. math:: + y' &= A(y, t) y(t) - Here, A(t) is integrated using the fourth order algorithm CFQ4. + A_0 &= A(y_n, t_n) + + y_p &= \text{expm}(h A_0) y_n + + A_1 &= A(y_p, t_n + h) + + y_{n+1} &= \text{expm}(\frac{h}{12} A_0 + \frac{5h}{12} A1) + \text{expm}(\frac{5h}{12} A_0 + \frac{h}{12} A1) y_n Parameters ---------- @@ -53,13 +56,6 @@ def celi(operator, timesteps, power=None, power_density=None, heavy metal inventory to get total power if `power` is not speficied. print_out : bool, optional Whether or not to print out time. - - References - ---------- - .. [fourth order commutator-free integrator] - Thalhammer, Mechthild. "A fourth-order commutator-free exponential - integrator for nonautonomous differential equations." SIAM journal on - numerical analysis 44.2 (2006): 851-864. """ if power is None: if power_density is None: diff --git a/openmc/deplete/integrator/cf4.py b/openmc/deplete/integrator/cf4.py index 79cb3b3d9..e93b22d5c 100644 --- a/openmc/deplete/integrator/cf4.py +++ b/openmc/deplete/integrator/cf4.py @@ -9,32 +9,33 @@ from ..results import Results # Functions to form the special matrix for depletion def _cf4_f1(chain, rates): - return 1/2 * chain.form_matrix(rates) + return 1/2 * chain.form_matrix(rates) def _cf4_f2(chain, rates): - return -1/2 * chain.form_matrix(rates[0]) + \ - chain.form_matrix(rates[1]) + return -1/2 * chain.form_matrix(rates[0]) + \ + chain.form_matrix(rates[1]) def _cf4_f3(chain, rates): - return 1/4 * chain.form_matrix(rates[0]) + \ - 1/6 * chain.form_matrix(rates[1]) + \ - 1/6 * chain.form_matrix(rates[2]) + \ - -1/12 * chain.form_matrix(rates[3]) + return 1/4 * chain.form_matrix(rates[0]) + \ + 1/6 * chain.form_matrix(rates[1]) + \ + 1/6 * chain.form_matrix(rates[2]) + \ + -1/12 * chain.form_matrix(rates[3]) def _cf4_f4(chain, rates): - return -1/12 * chain.form_matrix(rates[0]) + \ - 1/6 * chain.form_matrix(rates[1]) + \ - 1/6 * chain.form_matrix(rates[2]) + \ - 1/4 * chain.form_matrix(rates[3]) + return -1/12 * chain.form_matrix(rates[0]) + \ + 1/6 * chain.form_matrix(rates[1]) + \ + 1/6 * chain.form_matrix(rates[2]) + \ + 1/4 * chain.form_matrix(rates[3]) def cf4(operator, timesteps, power=None, power_density=None, print_out=True): r"""Deplete using the CF4 algorithm. - Implements the fourth order [commutator-free Lie algorithm]_. + Implements the fourth order `commutator-free Lie algorithm + `_. This algorithm is mathematically defined as: .. math:: - F_1 = h A(y_0) + F_1 &= h A(y_0) y_1 &= \text{expm}(1/2 F_1) y_0 @@ -69,13 +70,6 @@ def cf4(operator, timesteps, power=None, power_density=None, print_out=True): heavy metal inventory to get total power if `power` is not speficied. print_out : bool, optional Whether or not to print out time. - - References - ---------- - .. [commutator-free Lie algorithm] - Celledoni, E., Marthinsen, A. and Owren, B., 2003. Commutator-free Lie - group methods. Future Generation Computer Systems, 19(3), pp.341-352. - """ if power is None: if power_density is None: @@ -129,20 +123,20 @@ def cf4(operator, timesteps, power=None, power_density=None, print_out=True): x_new = deplete(chain, x[0], op_results[0].rates, dt, print_out, matrix_func=_cf4_f1) x.append(x_new) - op_results.append(operator(x[1], p)) + op_results.append(operator(x_new, p)) # Step 2: deplete with matrix 1/2*A(y1) x_new = deplete(chain, x[0], op_results[1].rates, dt, print_out, matrix_func=_cf4_f1) x.append(x_new) - op_results.append(operator(x[2], p)) + op_results.append(operator(x_new, p)) # Step 3: deplete with matrix -1/2*A(y0)+A(y2) rates = list(zip(op_results[0].rates, op_results[2].rates)) x_new = deplete(chain, x[1], rates, dt, print_out, matrix_func=_cf4_f2) x.append(x_new) - op_results.append(operator(x[3], p)) + op_results.append(operator(x_new, p)) # Step 4: deplete with two matrix exponentials rates = list(zip(op_results[0].rates, op_results[1].rates, diff --git a/openmc/deplete/integrator/epc_rk4.py b/openmc/deplete/integrator/epc_rk4.py index 98e8ac58b..3789f3698 100644 --- a/openmc/deplete/integrator/epc_rk4.py +++ b/openmc/deplete/integrator/epc_rk4.py @@ -9,13 +9,13 @@ from ..results import Results # Functions to form the special matrix for depletion def _rk4_f1(chain, rates): - return 1/2 * chain.form_matrix(rates) + return 1/2 * chain.form_matrix(rates) def _rk4_f4(chain, rates): - return 1/6 * chain.form_matrix(rates[0]) + \ - 1/3 * chain.form_matrix(rates[1]) + \ - 1/3 * chain.form_matrix(rates[2]) + \ - 1/6 * chain.form_matrix(rates[3]) + return 1/6 * chain.form_matrix(rates[0]) + \ + 1/3 * chain.form_matrix(rates[1]) + \ + 1/3 * chain.form_matrix(rates[2]) + \ + 1/6 * chain.form_matrix(rates[3]) def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True): r"""Deplete using the EPC-RK4 algorithm. @@ -25,7 +25,7 @@ def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True) This algorithm is mathematically defined as: .. math:: - F_1 = h A(y_0) + F_1 &= h A(y_0) y_1 &= \text{expm}(1/2 F_1) y_0 @@ -85,12 +85,6 @@ def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True) chain = operator.chain - # Initialize starting index for saving results - if operator.prev_res is None: - i_res = 0 - else: - i_res = len(operator.prev_res) - 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 @@ -132,7 +126,7 @@ def epc_rk4(operator, timesteps, power=None, power_density=None, print_out=True) x.append(x_new) op_results.append(operator(x[3], p)) - # Step 4: deplete with two matrix exponentials + # Step 4: deplete with matrix 1/6*A(y0)+1/3*A(y1)+1/3*A(y2)+1/6*A(y3) rates = list(zip(op_results[0].rates, op_results[1].rates, op_results[2].rates, op_results[3].rates)) x_end = deplete(chain, x[0], rates, dt, print_out, diff --git a/openmc/deplete/integrator/leqi.py b/openmc/deplete/integrator/leqi.py index c247bbe63..9d221e877 100644 --- a/openmc/deplete/integrator/leqi.py +++ b/openmc/deplete/integrator/leqi.py @@ -7,58 +7,72 @@ from itertools import repeat from .celi import celi_inner from .cram import deplete from ..results import Results -from ..abc import OperatorResult # Functions to form the special matrix for depletion def _leqi_f1(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - dt_l, dt = inputs[2], inputs[3] - return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + dt_l, dt = inputs[2], inputs[3] + return -dt / (12 * dt_l) * f1 + (dt + 6 * dt_l) / (12 * dt_l) * f2 def _leqi_f2(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - dt_l, dt = inputs[2], inputs[3] - return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + dt_l, dt = inputs[2], inputs[3] + return -5 * dt / (12 * dt_l) * f1 + (5 * dt + 6 * dt_l) / (12 * dt_l) * f2 def _leqi_f3(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - f3 = chain.form_matrix(inputs[2]) - dt_l, dt = inputs[3], inputs[4] - return (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + - (dt**2 + 6*dt*dt_l + 5*dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + - dt_l / (12 * (dt + dt_l)) * f3) + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + f3 = chain.form_matrix(inputs[2]) + dt_l, dt = inputs[3], inputs[4] + return -dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + \ + (dt**2 + 6*dt*dt_l + 5*dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + \ + dt_l / (12 * (dt + dt_l)) * f3 def _leqi_f4(chain, inputs): - f1 = chain.form_matrix(inputs[0]) - f2 = chain.form_matrix(inputs[1]) - f3 = chain.form_matrix(inputs[2]) - dt_l, dt = inputs[3], inputs[4] - return (-dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + - (dt**2 + 2*dt*dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + - (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3) + f1 = chain.form_matrix(inputs[0]) + f2 = chain.form_matrix(inputs[1]) + f3 = chain.form_matrix(inputs[2]) + dt_l, dt = inputs[3], inputs[4] + return -dt**2 / (12 * dt_l * (dt + dt_l)) * f1 + \ + (dt**2 + 2*dt*dt_l + dt_l**2) / (12 * dt_l * (dt + dt_l)) * f2 + \ + (4 * dt * dt_l + 5 * dt_l**2) / (12 * dt_l * (dt + dt_l)) * f3 def leqi(operator, timesteps, power=None, power_density=None, print_out=True): r"""Deplete using the LE/QI CFQ4 algorithm. - Implements the LE/QI Predictor-Corrector algorithm using the [fourth order - commutator-free integrator]_. + Implements the LE/QI Predictor-Corrector algorithm using the `fourth order + commutator-free integrator `_. - The LE/QI algorithm is mathematically defined as: + "LE/QI" stands for linear extrapolation on predictor and quadratic + interpolation on corrector. This algorithm is mathematically defined as: - .. math: - y' = A(y, t) y(t) - A_m1 = A(y_n-1, t_n-1) - A_0 = A(y_n, t_n) - A_l(t) linear extrapolation of A_m1, A_0 - Integrate to t_n+1 to get y_p - A_c = A(y_p, y_n+1) - A_q(t) quadratic interpolation of A_m1, A_0, A_c + .. math:: + y' &= A(y, t) y(t) - Here, A(t) is integrated using the fourth order algorithm CFQ4. + A_{last} &= A(y_{n-1}, t_n - h_1) + + A_0 &= A(y_n, t_n) + + F_1 &= \frac{-h_2^2}{12h_1} A_{last} + \frac{h_2(6h_1+h_2)}{12h_1} A_0 + + F_2 &= \frac{-5h_2^2}{12h_1} A_{last} + \frac{h_2(6h_1+5h_2)}{12h_1} A_0 + + y_p &= \text{expm}(F_2) \text{expm}(F_1) y_n + + A_1 &= A(y_p, t_n + h_2) + + F_3 &= \frac{-h_2^3}{12 h_1 (h_1 + h_2)} A_{last} + + \frac{h_2 (5 h_1^2 + 6 h_2 h_1 + h_2^2)}{12 h_1 (h_1 + h_2)} A_0 + + \frac{h_2 h_1)}{12 (h_1 + h_2)} A_1 + + F_4 &= \frac{-h_2^3}{12 h_1 (h_1 + h_2)} A_{last} + + \frac{h_2 (h_1^2 + 2 h_2 h_1 + h_2^2)}{12 h_1 (h_1 + h_2)} A_0 + + \frac{h_2 (5 h_1^2 + 4 h_2 h_1)}{12 h_1 (h_1 + h_2)} A_1 + + y_{n+1} &= \text{expm}(F_4) \text{expm}(F_3) y_n It is initialized using the CE/LI algorithm. @@ -80,13 +94,6 @@ def leqi(operator, timesteps, power=None, power_density=None, print_out=True): heavy metal inventory to get total power if `power` is not speficied. print_out : bool, optional Whether or not to print out time. - - References - ---------- - .. [fourth order commutator-free integrator] - Thalhammer, Mechthild. "A fourth-order commutator-free exponential - integrator for nonautonomous differential equations." SIAM journal on - numerical analysis 44.2 (2006): 851-864. """ if power is None: if power_density is None: diff --git a/openmc/deplete/integrator/si_celi.py b/openmc/deplete/integrator/si_celi.py index 96aa7e2ba..4f34bace6 100644 --- a/openmc/deplete/integrator/si_celi.py +++ b/openmc/deplete/integrator/si_celi.py @@ -14,18 +14,10 @@ def si_celi(operator, timesteps, power=None, power_density=None, r"""Deplete using the SI-CE/LI CFQ4 algorithm. Implements the Stochastic Implicit CE/LI Predictor-Corrector algorithm using - the [fourth order commutator-free integrator]_. + the `fourth order commutator-free integrator `_. - The CE/LI algorithm is mathematically defined as: - - .. math: - y' = A(y, t) y(t) - A_p = A(y_n, t_n) - y_p = expm(A_p h) y_n - A_c = A(y_p, t_n) - A(t) = t/dt * A_c + (dt - t)/dt * A_p - - Here, A(t) is integrated using the fourth order algorithm CFQ4. + Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis + `_. Parameters ---------- @@ -47,13 +39,6 @@ def si_celi(operator, timesteps, power=None, power_density=None, Whether or not to print out time. m : int, optional Number of stages. - - References - ---------- - .. [fourth order commutator-free integrator] - Thalhammer, Mechthild. "A fourth-order commutator-free exponential - integrator for nonautonomous differential equations." SIAM journal on - numerical analysis 44.2 (2006): 851-864. """ if power is None: if power_density is None: @@ -111,6 +96,7 @@ def si_celi(operator, timesteps, power=None, power_density=None, # Create results for last point, write to disk Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps)) + 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. diff --git a/openmc/deplete/integrator/si_leqi.py b/openmc/deplete/integrator/si_leqi.py index 71c603ed3..05002880a 100644 --- a/openmc/deplete/integrator/si_leqi.py +++ b/openmc/deplete/integrator/si_leqi.py @@ -16,22 +16,10 @@ def si_leqi(operator, timesteps, power=None, power_density=None, r"""Deplete using the SI-LE/QI CFQ4 algorithm. Implements the Stochastic Implicit LE/QI Predictor-Corrector algorithm using - the [fourth order commutator-free integrator]_. + the `fourth order commutator-free integrator `_. - The LE/QI algorithm is mathematically defined as: - - .. math: - y' = A(y, t) y(t) - A_m1 = A(y_n-1, t_n-1) - A_0 = A(y_n, t_n) - A_l(t) linear extrapolation of A_m1, A_0 - Integrate to t_n+1 to get y_p - A_c = A(y_p, y_n+1) - A_q(t) quadratic interpolation of A_m1, A_0, A_c - - Here, A(t) is integrated using the fourth order algorithm CFQ4. - - It is initialized using the CE/LI algorithm. + Detailed algorithm can be found in Section 3.2 in `Colin Josey's thesis + `_. Parameters ---------- @@ -53,13 +41,6 @@ def si_leqi(operator, timesteps, power=None, power_density=None, Whether or not to print out time. m : int, optional Number of stages. - - References - ---------- - .. [fourth order commutator-free integrator] - Thalhammer, Mechthild. "A fourth-order commutator-free exponential - integrator for nonautonomous differential equations." SIAM journal on - numerical analysis 44.2 (2006): 851-864. """ if power is None: if power_density is None: diff --git a/openmc/model/model.py b/openmc/model/model.py index dfa3decd3..b6a89792d 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -134,7 +134,7 @@ class Model(object): chain_file : str, optional Path to the depletion chain XML file. Defaults to the :envvar:`OPENMC_DEPLETE_CHAIN` environment variable if it exists. - method : string + method : str Integration method used for depletion (e.g., 'cecm', 'predictor') **kwargs Keyword arguments passed to integration function (e.g.,