From 9c8eddddd75f827d11e63d513184946425ff4141 Mon Sep 17 00:00:00 2001 From: liangjg Date: Mon, 7 Jan 2019 20:42:10 -0500 Subject: [PATCH 01/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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/88] 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 2758770153ed3461651f98dcd824cc328c782424 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sun, 2 Dec 2018 19:38:49 -0600 Subject: [PATCH 24/88] Adding some documentation about how to configure and install with DAGMC. --- docs/source/io_formats/settings.rst | 8 ++++++++ docs/source/usersguide/geometry.rst | 19 +++++++++++++++++++ docs/source/usersguide/install.rst | 11 +++++++++++ 3 files changed, 38 insertions(+) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 68848014d..e3bd83431 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -88,6 +88,14 @@ you care. This element has the following attributes/sub-elements: *Default*: 0.0 +-------------------------------- +```` Element +-------------------------------- + +When the DAGMC mode is enabled, the OpenMC geometry will be read from the file +``dagmc.h5m``. If a `geometry.xml <../io_formats/geometry.html>`_ file if +present, it will be ignored. + -------------------------------- ```` Element -------------------------------- diff --git a/docs/source/usersguide/geometry.rst b/docs/source/usersguide/geometry.rst index bb800699b..3315131fe 100644 --- a/docs/source/usersguide/geometry.rst +++ b/docs/source/usersguide/geometry.rst @@ -398,3 +398,22 @@ if needed, lattices, the last step is to create an instance of .. _constructive solid geometry: http://en.wikipedia.org/wiki/Constructive_solid_geometry .. _quadratic surfaces: http://en.wikipedia.org/wiki/Quadric + +-------------------------- +Using a DAGMC CAD Model +-------------------------- + +A DAGMC run can be enabled in OpenMC by setting the ``dagmc`` property to +``True`` in the model settings either via the Python ``openmc.settings`` +module:: + + settings = openmc.Settings() + settings.dagmc = True + +or in the `settings.xml <../io_formats/settings.html>`_ file:: + + true + +With ``dagmc`` set to true, OpenMC will load the DAGMC model (named +``dagmc.h5m``) when initializing a simulation. If a `geometry.xml +<../io_formats/geometry.html>`_ is present, it will be ignored. diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 48dcd8c78..4ab66ff0b 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -159,14 +159,25 @@ Prerequisites sudo apt install mpich libmpich-dev sudo apt install openmpi-bin libopenmpi-dev + * DAGMC_ toolkit for simulation using CAD-based geometries + + OpenMC supports particle tracking in CAD-based geometries via the Direct + Accelerated Geometry Monte Carlo (DAGMC) toolkit (`installation + instructions + `_). For use in + OpenMC, only the ``MOAB_DIR`` and ``BUILD_TALLY`` variables need to be + specified in the CMake configuration step. + | * git_ version control software for obtaining source code + .. _gfortran: http://gcc.gnu.org/wiki/GFortran .. _gcc: https://gcc.gnu.org/ .. _CMake: http://www.cmake.org .. _OpenMPI: http://www.open-mpi.org .. _MPICH: http://www.mpich.org .. _HDF5: https://www.hdfgroup.org/solutions/hdf5/ +.. _DAGMC: https://svalinn.github.io/DAGMC/index.html Obtaining the Source -------------------- From 5c7e8a9cf83f027cc9c5f5e15aedbabd7fc5f79c Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sun, 2 Dec 2018 19:44:59 -0600 Subject: [PATCH 25/88] Adding a dagmc note in the CMake options list. --- docs/source/usersguide/install.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/source/usersguide/install.rst b/docs/source/usersguide/install.rst index 4ab66ff0b..116938f65 100644 --- a/docs/source/usersguide/install.rst +++ b/docs/source/usersguide/install.rst @@ -247,6 +247,10 @@ openmp Enables shared-memory parallelism using the OpenMP API. The Fortran compiler being used must support OpenMP. (Default: on) +dagmc + Enables use of CAD-based DAGMC_ geometries. Please see the note about DAGMC in + the optional dependencies list for more information on this feature.(Default: off) + coverage Compile and link code instrumented for coverage analysis. This is typically used in conjunction with gcov_. From f1bc9f76b009b34e81a1a8d7ac2b6d876618eb84 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 3 Dec 2018 21:44:25 -0600 Subject: [PATCH 26/88] Adding UWUW workflow support and test. --- src/dagmc.cpp | 106 ++++++++++++++---- src/input_xml.F90 | 2 +- tests/regression_tests/dagmc/dagmc.h5m | Bin 1233316 -> 1233452 bytes tests/regression_tests/dagmc/inputs_true.dat | 1 - tests/regression_tests/dagmc/results_true.dat | 6 +- tests/regression_tests/dagmc/test.py | 1 - tests/regression_tests/uwuw/__init__.py | 0 tests/regression_tests/uwuw/dagmc.h5m | Bin 0 -> 1244191 bytes tests/regression_tests/uwuw/inputs_true.dat | 23 ++++ tests/regression_tests/uwuw/results_true.dat | 1 + tests/regression_tests/uwuw/test.py | 41 +++++++ 11 files changed, 152 insertions(+), 29 deletions(-) create mode 100644 tests/regression_tests/uwuw/__init__.py create mode 100644 tests/regression_tests/uwuw/dagmc.h5m create mode 100644 tests/regression_tests/uwuw/inputs_true.dat create mode 120000 tests/regression_tests/uwuw/results_true.dat create mode 100644 tests/regression_tests/uwuw/test.py diff --git a/src/dagmc.cpp b/src/dagmc.cpp index e6ccee24d..bc3362892 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -2,6 +2,7 @@ #include "openmc/cell.h" #include "openmc/error.h" +#include "openmc/file_utils.h" #include "openmc/string_utils.h" #include "openmc/settings.h" #include "openmc/geometry.h" @@ -9,9 +10,15 @@ #include #include #include +#include #ifdef DAGMC +#include "uwuw.hpp" +#include "dagmcmetadata.hpp" + +#define DAGMC_FILENAME "dagmc.h5m" + namespace openmc { namespace model { @@ -20,28 +27,65 @@ moab::DagMC* DAG; } // namespace model +bool write_materials_xml(UWUW uwuw) { + // if there is a material library in the file, + // write a material.xml file + std::ofstream mats_xml("materials.xml"); + // write header + mats_xml << "\n"; + mats_xml << "\n"; + std::map ml = uwuw.material_library; + std::map::iterator it; + // write materials + for (it = ml.begin(); it != ml.end(); it++) { mats_xml << it->second.openmc("atom"); } + // write footer + mats_xml << ""; + mats_xml.close(); +} + void load_dagmc_geometry() { if (!model::DAG) { model::DAG = new moab::DagMC(); } - int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC + // create uwuw instance + UWUW uwuw(DAGMC_FILENAME); - moab::ErrorCode rval = model::DAG->load_file("dagmc.h5m"); + // check for uwuw material definitions + bool using_uwuw = (uwuw.material_library.size() == 0) ? false : true; + + if (using_uwuw) { + std::cout << "Found UWUW Materials in the DAGMC geometry file." << std::endl; + // don't overwrite an existing materials.xml file + if (file_exists("materials.xml")) { + std::stringstream err_msg; + err_msg << "A materials.xml file is present along with UWUW material definitions"; + fatal_error(err_msg.str()); + } + write_materials_xml(uwuw); + } + + int32_t dagmc_univ_id = 0; // universe is always 0 for DAGMC runs + + moab::ErrorCode rval = model::DAG->load_file(DAGMC_FILENAME); MB_CHK_ERR_CONT(rval); rval = model::DAG->init_OBBTree(); MB_CHK_ERR_CONT(rval); - std::vector prop_keywords; - prop_keywords.push_back("mat"); - prop_keywords.push_back("boundary"); - - std::map ph; - model::DAG->parse_properties(prop_keywords, ph, ":"); - MB_CHK_ERR_CONT(rval); - + dagmcMetaData DMD(model::DAG); + if (using_uwuw) { + DMD.load_property_data(); + } + + std::vector keywords; + keywords.push_back("mat"); + keywords.push_back("density"); + keywords.push_back("boundary"); + std::map dum; + rval = model::DAG->parse_properties(keywords, dum, ":/"); + // initialize cell objects model::n_cells = model::DAG->num_entities(3); @@ -75,22 +119,38 @@ void load_dagmc_geometry() continue; } - if (model::DAG->has_prop(vol_handle, "mat")){ - std::string mat_value; + std::string mat_value; + if (model::DAG->has_prop(vol_handle, "mat")) { rval = model::DAG->prop_value(vol_handle, "mat", mat_value); MB_CHK_ERR_CONT(rval); - to_lower(mat_value); - - if (mat_value == "void" || mat_value == "vacuum") { - c->material_.push_back(MATERIAL_VOID); - } else { - c->material_.push_back(std::stoi(mat_value)); - } } else { std::stringstream err_msg; err_msg << "Volume " << c->id_ << " has no material assignment."; fatal_error(err_msg.str()); } + + std::string cmp_str = mat_value; + to_lower(cmp_str); + if (cmp_str.find("void") != std::string::npos || + cmp_str.find("vacuum") != std::string::npos || + cmp_str.find("graveyard") != std::string::npos) { + c->material_.push_back(MATERIAL_VOID); + } else { + if (using_uwuw) { + std::string uwuw_mat = DMD.volume_material_property_data_eh[vol_handle]; + if (uwuw.material_library.count(uwuw_mat) != 0) { + int matnumber = uwuw.material_library[uwuw_mat].metadata["mat_number"].asInt(); + c->material_.push_back(matnumber); + } else { + std::stringstream err_msg; + err_msg << "Material with value " << mat_value << " not found "; + err_msg << "in the material library"; + fatal_error(err_msg.str()); + } + } else { + c->material_.push_back(std::stoi(mat_value)); + } + } } // Allocate the cell overlap count if necessary. @@ -110,12 +170,12 @@ void load_dagmc_geometry() s->id_ = model::DAG->id_by_index(2, i+1); s->dagmc_ptr_ = model::DAG; + std::string bc_value; if (model::DAG->has_prop(surf_handle, "boundary")) { - std::string bc_value; rval = model::DAG->prop_value(surf_handle, "boundary", bc_value); MB_CHK_ERR_CONT(rval); - to_lower(bc_value); - + to_lower(bc_value); + if (bc_value == "transmit" || bc_value == "transmission") { s->bc_ = BC_TRANSMIT; } else if (bc_value == "vacuum") { @@ -130,7 +190,7 @@ void load_dagmc_geometry() << "\" specified on surface " << s->id_; fatal_error(err_msg); } - } else { // if no BC property is found, set to transmit + } else { // if no condition is found, set to transmit s->bc_ = BC_TRANSMIT; } diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 088377624..c497c79b9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -124,9 +124,9 @@ contains type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b) call read_settings_xml() + call read_geometry_xml() call read_cross_sections_xml() call read_materials_xml() - call read_geometry_xml() ! Convert user IDs -> indices, assign temperatures call finalize_geometry(nuc_temps, sab_temps) diff --git a/tests/regression_tests/dagmc/dagmc.h5m b/tests/regression_tests/dagmc/dagmc.h5m index 5f84286de26d73a7538ff2ff0558c95e19e6d230..7fb43ad4d9bec08107c19f033ba4ff1b6d3abf53 100644 GIT binary patch delta 1400 zcmY*YUr19?9NxKi=AD1?dMi>>x*4K`OxKJJU3m@Z51OnL5)-jnG$O%)P^1ebYX};% z-=YPUlTER(aPDEVkW2{wWHfuKxBDQ`gAWot2>Q^u_nhnI;A20|`M&f0zVkab6s(&H z)p<*p(jVayZRuk%2I9FKcr#;4*;~nhomWflz&zay#t$*LjKp)d5E_N$Bil#^U-FCf zzdxn4!e+iKnogFLo6AgAMn`@e{Wn=ED=I3@*2rc-1&>8VMk=x!m#AG(G^vM+ryLAs zJho8Av?Wd$N*597w7NT9O?Ttf{_I9-qNo3j@Odf+W=D?m$^wi0Mzzqjr$J9k4rx}F z{03em0v`Hm%-MvSW6V+)dvyrhEr%(&W?cz-mU_xH>uN3*T_}Rpv!Y6RRx5{N925C- z&aPJr=@Q4Ph5ZW0(u|Szt8vEAKx}mlaM@y3?O3#BB|UaRcT2A3I+;Yhc+E`XII9TN zigx2Xs|~<~R!Ozw1Ou$B%(wU?%qQf6)zu`lvB(MzIH;<`?5z$e~P-X9|TIlS_ z*Ptl*+s(+JU3eq7oLX1$SwW3epXe$Qca|t4RUXq!#}Wfh=Rs(-l+eieABn=dCo(4G z!Qak818Hw+PJlV-D7g{IB`a@0>>6Qc_Z!gdCf!58y#R0TX}X<70lwu^-SopTczM5C z8-o16$It=VmY~N?28a(>r>6Bwd&H`#w{L@`rh||^L`}3-HL+Sgz5pw)3h6BRSII1N$+Pgnl~fHkDh4h56^Saz^qVh@*^gFgU0GN z*$FBZuCTD!I2X%65Fz}iXx{Cjixs=*qKjP=1a{Hud4H#Kyz39x;Mx10_dL)0`+C@B z5^ag8wuqfitfmf}DwDz`zc~nRU1A%7u+78mY?bDB5?S|KwR1g_CqImO1tGj zTGEqIG9z%e@vWLm;E2q8-i*vlvsm7RQ4P3u!8O5W;^>o}i#Q1E3GHPm=LFGYwlS6& z#{Bi=j8x$;%ZFEW9E|qtLb5-+LeLog7$WM~9 zqe9PiLYJ}&p2+d04i;P20iB+DGrKP1CGu%M4EIzjp{F&}^xK(tq{2X!H+5pkotl#K@0jfF#*$U)*FmM^fjLIC zmHsOkgb;RJAVjiXj`BE|@kJOJWrKnE=mNaDQjIQnyl8<>R&+t~`ZIWRwVAanGn7!| zUpC0|DE+Zfs);O)<@MOfIw+f9odoB^0$jXZjW68l3n3JT?mIZC$i2&own?QpC-SIk zRS$*xnwa(n?TW4)hNBoW9#y>ga|y;Lu?5LV6c*w{5=FNj+tj)Of2NRQODIbB2ldqX zd;wOTRHH#_KC%($w9|^f^wS_~VIgi%a;kx~Da=W|dBVcepT4sHFt*~m471prUd!Bv6dVYhX?vNBm!Bgld3>2FvHdB;P z7%5B?r4(fpTPU_tY)g}r`Sd8KkA)m$Z(4K^p)~#5OZxfPk7bK7J?JBK9Q?Xem74#~ z6;bR#(E~*f6g}|&9(XW+AvW(LhIDq;B5?TAdUMWVOV7}+(lq@K7EQ77RZI9U9zwo9 diff --git a/tests/regression_tests/dagmc/inputs_true.dat b/tests/regression_tests/dagmc/inputs_true.dat index df8e3267e..d88072223 100644 --- a/tests/regression_tests/dagmc/inputs_true.dat +++ b/tests/regression_tests/dagmc/inputs_true.dat @@ -8,7 +8,6 @@ - diff --git a/tests/regression_tests/dagmc/results_true.dat b/tests/regression_tests/dagmc/results_true.dat index fc13a9cda..134cab089 100644 --- a/tests/regression_tests/dagmc/results_true.dat +++ b/tests/regression_tests/dagmc/results_true.dat @@ -1,5 +1,5 @@ k-combined: -1.115067E+00 5.423808E-02 +8.847913E-01 1.700903E-02 tally 1: -8.543144E+00 -1.530584E+01 +8.016189E+00 +1.358910E+01 diff --git a/tests/regression_tests/dagmc/test.py b/tests/regression_tests/dagmc/test.py index 9766c9732..34fdbbdd7 100644 --- a/tests/regression_tests/dagmc/test.py +++ b/tests/regression_tests/dagmc/test.py @@ -38,7 +38,6 @@ def test_dagmc(): water = openmc.Material() water.add_nuclide('H1', 2.0, 'ao') water.add_nuclide('O16', 1.0, 'ao') - water.add_s_alpha_beta('c_H_in_H2O') water.set_density('g/cc', 1.0) water.id = 41 diff --git a/tests/regression_tests/uwuw/__init__.py b/tests/regression_tests/uwuw/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/uwuw/dagmc.h5m b/tests/regression_tests/uwuw/dagmc.h5m new file mode 100644 index 0000000000000000000000000000000000000000..6aa6f058a14d29516932d2b825054a61b2608eaf GIT binary patch literal 1244191 zcmeF430%zE|No~|Qb~*UnyH3_Bq2n{5+RjDWoxHMp%7V;rIelQ`sA-A7M`ZZb`(E)=in`%gS$gq%zXxeOf5s|?mlAWj!%p{*jj~!mgB;O4U zi;9jMIkvK`O7>|l-pkFslbpC$r|sG?8~5XqGT(by5qR|akIT{HH@f!YaMj~e*uhLh(T_X7%5+A)_dxO^hHsj|uVO*+|6+=Abfkp1e- z?^Oe|AMhs+LK(fiB>aZj7tIK~Pwk63d7?YE_j2yy?&(XvCdZQ(|HsBx3Op5_yN3 zf8anm6x#o!k$qb37TGT&-3{$mtjP06%gH{E(pT}H?^V?g$kyM_#shIp@wJfK`p++{ z`UwB18le4v54mb;pO=v5gWimD@ZbPHz_XLei6AZeQoI*r0q9cJin5B5)>H`{q4LgC^RfQcw`6}t?_un1H?Px zDH?6xv|p!3t*ZSlA?Hn(4BW)8pS53)KwPf6{o;D~xBZoEzqFm*BJGFPM;BN34((n2 zx^^Ofwo60{{NfhRBjjNok=FRf2cYe9IBDPX{Ic!+$4B}7KCP-BSfhvUl{3KGpY;Pi zVYnQv=ijfUs_*DWHBePQkVo1N{qFy!A5h|1krRY=plAC2miA9R?Z5vyrF8h{f++d^ zmW~@^I)49iO6hYEL!K+z|HwXPwEX||md3+m{h?|++saFk7bzmKW(2kDC=e3E;Yd>_|=3Qsb(loS)sEhQUX&rt!A zkO|KtGJNdd$OxWSNac`+=P`CCt6Mp14aB5v4Kdr7f$AElu7TKdr7 zf$AElu7TN5CCy&8gcYGCeo-eyk{VR3&PQ*8Wc=ok;_Fu^t zn12Nykoi|a_u?bRjUp*GmAkV36*V6A%*pX;R{E`?q60=oS9-4Sw>8BM3G9d$|#8j$s)^U3X5kvuJWw=^Cs zjsr)81Vx994vQWu`HZ$->X)YFf3~IZG}`Vs8~hD=yr*PeMJ8X&ikyg+^C0`fTmS5g zzvt&w#nVb8$y4H&CV9Sp)(@lxiR9| znZCcJ?K6Y4Z`u#Yw)gLU`}d_Zo+jIGs*!Ppyqr2uRXhSA`U*dZhC&=nqG`^w8pYA6 zACI$ANj!wcS!DMW!~q(~!DQogRgbTpTH)U`?&G(=&VOrsHIIx7c_gm*?K%F;`>MuS zCGGIHUC2Rx+P}8N1!?>JGv}@@^ZV3*Z2!t9w>DR1-16TZXH_*m&}@S*z?QVjpS54! zU|hbs{gNN0>HAySKA)23iW)ZA_Wu2E|Gt!t4`k;{ zzc@Y+J4b#zKF~q94yhkTcK_GM2cZ^NVwfkV{LRnv-x?pZ@}mc}rEh*|`?LD!MUF=E zOJw(OKc=#D<02ZoMpsb|E>DyLGr2o`grNbc~6>re9d$V_b{@DMZ(d=HW%%HBWA1jFb1|Mgdn_v8AQMJ|kPf4x4^8sPc}rA@5zi2vUDa1s8q z>I3(u)Sglu>|7VzFhWWD_}Tc3-xU|c&k*@i_2yqw1G3{YY$JHQlFrgqP&w}xBRVc3 z%}lz`FYIAIt`e+IvStYWGCh6<*?$_2xfMBf37MDEa{3nya;aoLV;ts*c*Ap;Bwfy*?8i+e1^z#6=yM<14POAcd>+-m_YmvfJ-dn( zExnw*Bqdm$v~p#yAqJJb53dD67Mk=PqaF0SnN@b4-|Y%-5m z=2CgM>Oa3&19aS-)(L+xi#(5P+^tgSmkSOV6fi10TI`aeHBEmp$%;^t-`lG40P=D* zD$5TJj|>b5ud=xKcmKX!SnZQ@~hkx$({@qGd)gSthc0}*@f7Kr*R*wFdnv0GPJ9YK!+R??; z*^Bgu|JUO~ed^ev6;QqXQ#C;QL)rOH)%_u@=|9;Ye!qB?{UPU#FZ<7q4?TQvc{-o? zQ#V~*=(ni>+5V7E?ktqFBYLB>KMV`@t6WlKbU^qh$&xs<6+6&^e{xI5hfBzcaI)jW zF3#>f{5o~*;q2C{tGjrt_Mh~JlKC`k-#?x&OFqRf&Avar5juXb^7(N&dR~6NAD8=G zKBfJkRd0MpE~Fh*-5=5#q%G+u^X2as7q`P7`|0U$pKgEMPhUos`WsGDXDV;-pX#Tl z^+30h`jEu@X!@Yc=q*aTv`vrrTd{S9#;=R9?~=lEI6w@1N4H zmIVD;eO8yoHSn``Mf23?_Sf50+#%eqJYY^MIrBfIU6oX>h*jNMezgXE)~={qlx}~$ zU1c58s=O$?{1v^ve@eS@p%Wu|;nmweRReP3tPIJxhmND^_Fej|+^cG@amjBWsdrqbsIqpdo@W{#!S;17bzxMN|GV1P zzfvFAFDu=$s;-ai$8de1ALE-?fv8BtMi6KYnH>^+#3y>|k&4UiII9OASCG-pS4?Aq&$7Tq;9@?=eMtcs_fQNqI;ng@V~NKvO1eT+O4LQm+*gW zxAZ5_JOBQc+O1Ub{CvL@^S@Y4tg6KC#sBwq%Q55MVYj4PR@HXPmDC54YDoV$yEO%0 z{-3p5w@CSFyCwCr|7&)tvRi8=OIN{v+iqdMoFZ=tISnHIN$R(g^;bI3%OUxQw47{y zVk#+@F&n?3PNjZl_(&hZgejL3pQtt1%q5e`?{VcWro!fVG_Uh`^#m~## zRdNr1-Tq5vu7BRYE9-aVPy2BV!Z?8RM^%;9Cm+~Oud2lF#eX&qpzWJ(e|`Qslh{4_ z3*h?u``W!9=Lgezpj$S7O_gd{RgVJ}{=@3y`}0?S5?A`O;{g9;{JYcV>-Tz!s_*7! zHBi+!AeF=m>G%KJ;{bYG*>x3Ub;Yu}I$2$b?0RP%#p^V5arY#u6dq63f6Mmqyja#w z%Em=2&(}+=VsZAflU(yR`)+>DZeBF6M2c6I$Dc?%$CB^4cjzG5RK`jqTl#!T-v2LN zvUWIi=8w-ptw4rMTw3;9w4X?yCHu+mx&Obn!z=z_cIU6!Vd<7twH@B@537&w>rxWQ ze1kscKZ}1_&BxCVeZGFrr>Xj`s;Ysi?656qNA!FDZ97c!+J-N{XGkR{qVvda>{_M$ z`%mYQR^<2!Jbq=KkhOSU5knc7Zi=edN1(zBi}DK%hzJg+>BzF{QqlPp-Kz0iDl1H)B;t$zr|T&F?02Q# zm#yue*;`%m1-kw9dGB^&Uua9m2mGhbduct;Et~gBKLJ(STd(st9^fJMll*<_ri>?a^hc3Vy8z3j0_71 z3=i=O4~ZBY9r}m$BWY*Owm+V?nUad7+h4b<8!liwQX=&e{eA7oU#t&t!^dqORo-XS zc6DLq?^GZ7Zh1VqW#bI#9-!*_n0gV{2MqHo>+kR0K4`y0`)&WVm6y#MjyJrYjWfzf zd3qf3wd&2kqz0-QXQ+^NL95|^HO@fm@fm!`iM~6|_{LxTUms_PJIim^LzW$9EFsrN z>p^yl`&u3^OQybL%UKbSzESoQ>3#Yfrf!ht|Fk_i&d6AWul(s6yk+Bzzw0<*>gC_b zF4FN5-NLv@+;DN*MU_8H)pjxZ%I{Pk_-?;9PLO_rtFDjWEHX}z!z~*B{rip+z8{Cr z+KR8+pWMXH#tA;#a6#HXk*`&6{v|aaJ5Jz}`?lR8T>*a_C;VsYIAVLO%*!IDSwrjz zefVVUiwnIjvfdFbC%gU;?)!MWM6&KaEhoDU|5I|jj7|6rEmuPJHOc+c=OVA0bUE2| ztm*8N+bLb{DcM)qO-}e<-cWlid)?~93eY@&zdpW_ubep1dQcD_SK`N!_30y zWX<+(KK#yG`a9W@{^vVIDiJ+z+K#G9{9gQLaR^D>OI~yfc3;eJog&G{;?GrZY)fAD z-?b0EitU_LWjXri?`!9N>@T79K)1|4Al*M=&KaTV-{^Txx zX7~Afa6$SUlCM>7{v|aaYxgzDeOt-cDQfqreRTN`>xbGaS-VK>70p-sReMEa`1Ies zwI$c7ApZArBz=PGx^5sA;R09(S1 z$rgXgYlKEADw1*#FY_M;? z6R+}RkBpqJZL|^d;q^7g-`RTx_?C&Y%dNf$kj*KZPYXL|0aq??+j_V{fbQdJ-}waix%vBQM7}=J}ul0#s^x^ibxya^RC6YRQk@DL|q8 zxZN$vJ_0|jyJ)arlK?%Kt8`7VM+I=5U^VMOD+TCw+kOLV)_w*ac*b#jr-cHfquXQY zw$d-a^aRM|_>tnqW7a4DQFzeO0(E^kd-m%%B zQF4g3tFzAax1AU&Kr7Uq9NcIk4_x!?v%2g31<2*ru>9CFj2{j=8D-X6fZqCF&3_!O z0ObpGH!HU6DnRA6Cl=q+Rs^1#dsCyxQGlX5MK$$5%s3jJd)vWQfD{A0pAYd@g7Ug! z5_VZw2vCQtS$)ctlz~qvD!RK>SAZOgTzAyi#`u6f{yQ$$5THr9du`{t*MRbgTJh0y z%K4~MOrAx`TgDwCtf#&x;v>rf#Y<|7RiM1t!Z%^344RE;5cb5_+9J1GDY{79IntG!g!;h^IW`6nxaz?Cwh1x0+axMhMGo#jf31v+q?R4?@@XtN0#5qWjRmN3Tiygx2Z z8)bqvt!%cSB#rSkYa;uO>S%&2w*)sAePF!*zT}s4G)<7uJ=5<2;R8pE0zZF;czu zdHNJj#y!;<8#t>Nqq1c$FAQA7c!PxA!Le73(7KfEdk@@ZT;rN!X5<^pC4bhag={Z{i8K0kI)NB?qL`TMC)NQ|x@uz`PJ#Q-- zqGeqabzL4ZJ}3CSV0MN9@-ul9c-+hYj(<8*hqoiy0EONXrmP#l_~?LVZY|vm(54y> zI=)@Tcwd#gC>u=!bf9jm>aLrNhbbSqx!|fknm;%r;kvdV9RJ3unFXQq_0a;KuLAMm*im_?Nw#N5qfRMH7llAC^92e8Z!xUY&AJ&NlV=%97k8ye_%Fn(oJUv=x%I%wq9OFJy%8Sm_p(&3D^4sx5TwRpx} z#+y&Ku`a5ugA&7=oo)Y^@r6B<2Q53NjrOm6i(yC97HH^o{^jiC3troiYrpA(#ON>7lUG5v{ zt%bs(6AoH`WPH@V(V^4Tw2*q8z}>a#nZfnUGmD6`KT!uMo)|K{WlzQ(Yd`Uh9aRVQ zd_J$uj#$R?+*@t#-=GdkH!Z$BbR*-z=6Vh8+|v~IdwU{AUS>S2apbahb2L$yYOh+{ zd&aq-@WW{xG*O7l{`>~D`EWg1yRsJ=7S%@M+JyF7WX||e-K$21Zne=#<4a#QG-v#p zhiXdDZVhzMS94%*7slnI?#}tBt%2ef#H(%Y%lMXsJ6GqAQb+Q7%U3TAXI!~vVL!8* zYADd<8mc89e%)^lb7xlEu>udY@6L0#zootj=y`aingY% z%-)d3c&o)<1Gn^6MN6C>74PR5FWT4t;iX;q%=D`sYcJG5Uk5fzu=>n+a2+-N_^pF?zIF^&(55v*mv_(|(-~jlHmY;;4e7A2vu>yK{(R+?)1(^<94|d$ zT=kgJ`T-^{oW`8nxGnf8<7-;4U#&Bw$jKq5=lccE8FzAA__gVQLMM-#NU>iD<7dCb zjyogw#A%7*c+sJ^jCbBD>5XdgA#y*yN`K?HA>)xZuIgG0ezT72gWvM9Qxi3uMCe3z{VLM^ z8m9Yv9^z(Vc~053&J=|(#@D@=Io1D_JonD!&4$bYjBA@TS+`@60=KU;GQ-A`@sl}5 zHxk<_a%CA>O&fP++}!%3=gUWmoO!(+mgibA9=qMd>&^@%&VH?LaYZA>qoP-6cWt7~ zP4FLk<+2&$qm@Fgw79O!)gHRXq^%a?ORc>!R>#!fJ{G;ObWvnHgzMHJ&RB&D5A4&T z|T_0&0ID?8K3QH)Qof3}U{8Fg;5f;!qB#CWi;puw3y4Q^H) z1FfWc>5WH7*G0IGN!1Y4p$b?%_F@H)hJ>d>P}5UCu7pbyA0O+Wzrk(=f)9 zhQ;Lh)zszelsukhwqX2Y$6Yxo?z&uR+Dp@;3XESqZFJu~MVE`|aDJ&m7Ms7+Uhfnz z^NudZpTBw9`sIuZ;zwWVYF3w9^~6kjvx@Eyp5tlQiH~NT(;Wir$#dwKV9hL=@qHRiKc0M zEbPU&UGcoD@6+@+*N%Rx=9_|xax3lhxjn=Dt}fop z_(1Pt+{%&q+?tO%UtWbXzNjqUQ*)0#H|s*&gq=+p_g>kt%i9ut?*0bH)a&oqJbF<# zD_iF#2ApzJpF2T^81FV+xL$FX0q3pJGitzi#x)jLw6fi5z}38^nsc%p<71Z$ntisw zfNP_BE;NP5__c$JI-E5(vIY5jQg#K_pVrC$gOJ} zIbv%E#)GUcsMuy3a;EVWib$UE9JzXj)4T<7NHT~9FoM*R9zz6T-DR} zNw|eEmlJzBc5*c1uO4JBZ|`f&=`~h14YB=(*LXZ?#B5`(^qR?mbH&V$6xp&@yH6*L zxr}(N?5&#^=hnRFp`2NVPaWrmFR>FL`uwZ93s)nT>D z;OoraHU9pKnzLt{aCa`Gsq~-4_~G@#ADk1Ja5uKaPPy&O_?KEM_Ij0=aC`6H(|`Dx z`Q_qP>W4Sun{sQE{DwyFXZ(23lsAFyrd-|wTYgG76(LHp?IUxgO)u&QUwxxbnHqnT5p$NzDIT z%k%x}Mt*!QCFtvosn3}o+vQ4jzZxUN*VA#q;r4Efk2~99TztIvIOpZ;ymm2v_1KVdl!DSrMrMT;(WhDB*)@7_YTx~^rEkT%b$F)Q#%93U7CDRb{Zt$65hl%Uo}D= z<`ESZfuF2Lh_8R(@kjH|%K?9N;ZWS}SONFA=M6qz!1#8@E56sJ3Amhseq7fW9+dBX z>qz%ma|B%5IsvxzSABKjRle-eyQR-e@|OrW|1QUP3f)=!OD>>hD_9*Fl^*n`F67F`r+*dG>f5nE$!=YC*u zy_*ZpbVxWM;P`J3Sua}q0ghw3ILtmOL%^klXKI}3#^Qt5>kU2f;j(~>uGiCjJCeTjz#=d#V30edF@CaGZ-Kol4K&5g+Hm4a315S-etXRsR6JJOS783BSwnA{O_Y z72fAg#bW{2eN-FEK^CvzIQDYs9mYKqa9_ur+pui-bKsK|LgMDU5OA3}O;0p9_6+#I z{8h1zuLayz??fAmnuWk^2P&COE){UfE>0U4c|8W+L}{`2fDZyLPA(^J!Rm*=pRV1l zbG2N+^?mTZUGM9;z&|Y4a=G_Oz#UZx>3Adb7VyvAZa$g#S-=ewyy$hV;u`P+Z6-WU z|03WP-ke&&n5k8F78)>(|-6M?uUWXetR_T zw}GpZetsqH=YdnZa2DGI;M9%`#C8O@6R|sa*zN$|LG08dY^Q)zyS5SAHQ=Sh4qnD~ z5ID7)cGzwLr*`%#wzI&gUG9zTGH@?q$NOMA4!n%m{pr~51E=GJ4tSgZoQ^B>@wfsw z9f!=t;}GC<+|nM8TY%GXjvgN80G~s~MSJnM2sj-_?Zx9L;MQc^RR@o|fKMgkv?M%E z15U?vX?R=*Jdun8MR*(toQ@mA@wgFqXEM&5jmMe5>9}+V9+v{AV#Pb5+ zbRLn4=Mli^ykjq(cK|<0<|#|?JO%g&GOsy{=QY4zk$KP?JP%@=%$p4Hya_m+XJzAg z7I2w)8SuJf9(NDVttS>is!|^W2NWOz%|LdTN%&0fxD1-dS5(G2fl>N>lN|59=I>@1Kh)Y z0N_iA-yj+L4S@R-Kf@~QX8>MG{1R!{F9BSI_%R}|9|O2O@q2W`eh=U$rTrwp9i;s# zz^6(3VSrP=%{uJ20q!X6=K(I`7Xn^~_>oMp9|?FC@jH#iekb6ciJxi{_EQ0$Decz+ zzJmC{+G0N#@N#Lt8SofsKO6Aw(tbJMX~d7Gg#CEHHHqKP8~go$XAnPOAM7UtUMlTZ z1m1@DAv3Ta5_nzWw`_?0mcZu`Kj(Yw=LG(c_(i*6zbJ5}aoCTlf&Hk!Q;6U7CHA`l zZz=7k1wLBZuM7MZ@dHoCeqi7Q#Bc0?{l>uarTxsn{iXfVz-9c{z-9d2z-9d8z-9dE zz?G%_@W49~zx{3Ow+C(~?dJzRgv142;eT#z*mwuMLCXB0G~wS z8X7pR0en1(gN(;<5a3g#<0inZrQS>fpE)@R`zaF5p%qF18HE#emDi(SXau-GIx)>3~N|$Mt}7(s4lG z$AV!+79@`)0_U*+??&=oI^et) z;Dyq8GQcO2yc!*xR|C8+$-{Y$^KgK7CV4w?IBy5|2enlIOM!=eYrom(Gg=K2N8k^o z^CW@GOl?@cAUK^&rk` z1#U<3V3TkjEbvYwZ`K3n%>ut9oo5TYOgb+Y_*#<3>w)ulfme{c-$^*{7kE>WCoI5u z!oWo&ulNhjD+bq>xA9ZCe&p#V1^x)p1zTLG>ot#bjsTUr+bd;!tX zjK?|};4-=!;4(TL;4-=%;4(TO;8ZshhIK>0ml2)O9;`D0PIXD!ur3Ms2%=+piFHiC zsqU#3);$5II;qZBCk349s$8+I3OLnad0-tD@JOQD^1!++;8f@J1nazjQ(ag&)`bC| zM|5OSSVsn2M098NSa$||6Va&+#X2?M&55pU6V|l>UqEzlU$71i_-vw^E62Jy;8bT< zfOU4jTM=EJ8P??ir#ik4SjPvP>i))K-5+qO6KsWbg1}{Tg}`NWh`_0CaWvK~0;f91 z3aoPkeu(HI-(g)OaH^y1hIN#{&k)_^8mzknPIa1du}%{>)pd@>x=!F!2O5lZpunka zRDgA(5a-e+I@4!ZX9}F^QeR_TDsZY}9fx(Sz^U$4*{o9c3Y_X>W3WyZIMvnm#kyMH z`9z1i2aH>1LiFL=osZN=Xb;`i0 zu6Z`rH3O$Q=muB^4O~Vy4V>z%Be2dIIMrou#ky?ZRL5N(>$rha-FIiK`vy*R;@hxJ z95~gLFT=WW;6Bq2FU$NYM1AhxoOM+$-HBKE66alPdoJ1MMy(-hDVu=*JYRb!9>{ zUS;{K(6217AbaV_{^jq4i2p=XW6T$phw$)>@O|MMA#ztz=LZ&tn#h4eCO z6j;u5Ec|!K$iD%%0H=_)uB?-~0KE4jU9x`6& za=Q22I3YUOC39)>mGX$Tp9T+R)N4Lkh-_Ae_Z0cCJR;*4hR33Y3el_6Zy!GYMkoGm z$sSi%e<7-+efi4a%Pddngu+Nwk(Us43%6Y>oWb&%TwVPdUvv?oxC`~~H|WUnqcmm& zzl?GeqME0-Xr8WMT>Fq);(;bY)JbLI>5}~{|0-|W#g+@rg(&%2oy6^8eHFd_h`eEK z)OCd@sh525+j=bT%evp%Y(Whnnx(mS(y9B5t48hXaIMS|O?Fs*M{g0!A3L0K^`dE! zCF=EZ{QMWLED!C2{MUmu?pmUN#8Ifd0?S)_H_^m(R;DF7dRjfM_94bKDhjUcI%0`( zUUz-v6~Xf2lG?mj9J$RB4Xt(9M#F;T*X;;d^7P^|OVl=V#>1938CO{<)XA7(iH7f( zyW?^a%hOXHv2kcXj3v@2OULiy)4o01_SEx zt}^cACK5eqZ;8;kEWH%5UXa%3$XdL*pKL5qm*AycZ%t-=>!p1zX-1al+G$t4jN6R2 z?r?Iil8PmIn6G|6uqn%v?C*YK->0`0DEZLBicN`(m-Y_WHvfSITBXt=d;SB)-ws(- zcHx2rnxa{J%CZH^+iaJZ@O8s}3pBX1uU)^Hj5ptZ%Tln$0(tW8OtH^rys)176{jQ% zWN#2}wyPz}FD>*QFi0-i0&Sn!AmaQ?#vke(^BOe30)?nX-I@H5@#phzD@As;K&5fn zSvk#FzUx*Gy+WfV7D#^0QqkENjN3PyJ9wC$1=`f}U0sh{#tYllDGT~)jxw#kx{UZH zFW2;dXI7DHc$zjIBcIaIlsArBmHmBb8sLArKedHB(J{0S| zI_k8#uI9tI|CiYhHSz?=Tw{lh$0o*~*mW!^5bJ`h-nFqZC}sRS=jcCVvjAN(4q7p& zJ+?f*8_mncr+wZ@i z_kfoE z+Wh>lj2qt)^|E_uhLT^%ncQ(^I*YNlbmOmYGeZ{P^9tK1GrnPY(e($RW+>ungXa2K zjO({(nA}@vhSXja9_y#WbRy3O&MIi}&=i&KT<=)Rhw=GQGa61?VT!u(-$w|=>q*dl zrbUF;+V{2+vYymu?1gW1$JfX1@2V=+dpw)Cx*$rt9tJ(m zxG9bG4W^o)xh{n*vOE~K{m{^TxuXdxPp_GEe>meg$qpUoml`AIJ6;B>=Q3{UHnRKI zZN|vvNOG)Z8sm=bW}J8wXpBOfp84BdVZ77BgL3T+jM2NN=FJtq(fP)mJvv1Hig;ax zWkW-zsR`gZ(|Oki44rC(EPIXi-*3(Mo&{;`+O#r4P7h-iEpTD{`HfLgGhP^?w7E88 z%?2^P`D|LS_gX`A{mPwKpDB#b_vtxdlb0d-e5t#3yETlf?p|*lQ^OEV)ymy}faW9zEl$uhhx_8^ zyqxjwl0MQm7$!X3kn!W?&vU{i>LWF8xq4eWGhX2`Iws6UABopPHE0#UxX+^H#qV7ZlZMAK{<_cnU7EA>&`XcpT=x{l?I-VWKGH@HmFliKn6s1dDFJ<&et1?Foobn~ zKrw^yE5S--XBO8*ALnM5-^yb=_I7JQ1LwLZ&QK@1S2^QzHjQav{zezwH$6Tps*VNX zRle*|NO|*9dsgY9Eg^;34hEo-KWoub%^oKy4S|7d9Q=y>rOS_?alaT zTXdGUP6u5{S=A#bg7JvC?JsDbB^9MlWqjmTkET7m6Zo;qrhjo=Z~4jc>Nt4VjC*fMg=#rA1z$Oc>Z7;yLukt zbqN}mSa@t^{PYXgX>*sWqdr%>6t)~>JW!|M)nNtVah{!}{pxd!&l#~G`h=|-T7Ij~ z5WAaPh zg|s*Q11B24^{R;mH~x~bybj~HZ!`)&zd{AIpW3^ijTz$(cB3P#AJ;(L-yI)3rV-;$ zXUT6cXjlV1QD~;&){61HhrA{(8=;KaEKOP=a%OyI4l1F9`f)cidon(K$&A+n zKP#far_-j59KiUnkmxa6Iw&HC)+;0ThBBVxpfi3!f&vQH+~zksit&)(gZw~+8whz^5Zhb zhj)z-l7ohJBBS89HW@z2dm_Llc3b829e*HiTZiiF+Q-_mhNr5-Z{N=8WcGphw&^qzxrIxTc^95ot2F4F^&#oPATwt<1{^R zL9E9^#?@Y&jop8z#Ocb-Nry}d7=KWEo_A>aOQ(c;4O8Qb7(cgSR#5igVkbYdEq0S$ zGM=N+pzq5y3SF0DI zvMUWg&3evwBg2v*6OBql8>93c2Rvom$#3bnTT5SyhV-v(arzPC>l)qNdtdLZXnyYx zXHMQ{yda_Grx_`wqSDKo9er;zzAUvti)^j;qG3ayu1UydJR?)P$KE*~L^Y=m*KoVU z_@*xl8(657iR zmo)LSNTZAQoW5y{pE2`$^rq~Ks6AS*VU)_aSEo*wgIe*p<8##qMXzSu@ul-4o>)qX^#~wMZ{=82e-^^rud7W?@=LmVu+91_jZwlk-ZF>ekX`sLzR6I4~ zMhxSLCY?R4pD1vFTkCCY9Km>f&Skvj8bz*I^1=gVL5vSyVi>xxuM+pbf8V2R-i)tQ z@40f1jxtwat{r-{;$^5d4@qV|I zUcDbso156EPmI=V#v49sxpTIPCO5EKy|$;vGd^w0+i5ofG`T@Rg+q0QGL9OL>3#XI zCTE>+%SXYN@oN?zqQ_|0;riB4={v78w9*uI_Of}Nt78R~?D_3KDe%ZK&_eW}R3rD*PYLv6_Sl2u6Tj!nE;#4+_xVrx_ z<9l55G86gQ+?)7|BBdvBPt96FWTGj1AXG-^zSF1Ki<|NH0E>`fu9Qyyh%^vcfb9FhVtG6r0UT6GH+0?QEEj=#v)ZL3QM;M>n zy6MxXU_I`!Utyr>I>rO+qleDhug9HH@%gV74`{bZpF^`od~M%~amUe{8qIm9&%J-3m}+avc&Nz{{jg33+@;Z% z6wb>tzId_I)(i6uxDPR2DYXjOybR^u8M-IWfUB|lNM6}D>#o!{izza-G31u$9T_@L zyuK-~@@0=c?&z^!Bi4}X+h^6W8{e#dav|r@fVgvp-1d$g{LV(R^4Iqn#f;Q8;#Tz0 zx*O!p_`VO#mpvM2#BJABKDP6lbx0Pjy50ZU79&o5>s?tu9*tS;a&d)4ikn0=fZ!zqn^)@`oOuSA- zvZrVV0f=Q3x=t;tVey!Kw#*!)mFS8||MZh{-*?w`*OI*=&VO>aImd95VBXh;|bT&(w-{{ zI1hz@n!6q{|J!w8igHs60XIH$Z<+aG#@~%_I%3gY!1YlZA*X?)AxzOC<;r+KWzB#&9x3+K1Ii2bDK6iW>fA(y*a&#>V zE_`F;r5DPKKi(gv=3!yM#hf$sUb2_@$rnGF`7Eoo1((p;F8;0`<13HUw|U>wf?K?O zdTx3x#-q=C94ZR6;0B3u%bOl#e*3jg$G2!S$%0E-kQ3$-$oS^Lx#o^37F>F_f)k(B z7{75}DdXi<3(jV;$C6$LSsdYp!;0e_ISX#VakFXN0~jw_rnt8C4GV7geG~ISRmQ)% z1l7uZCO*#ai>g-NtXqpFuUfA5#e!3)J>XvL0W97z*xmniK^;pjXZM!oYm^vY?5)2u zKxoM|dE4mVogFOBa&gJ$=fn?xJ7qME_m>Pqh6m`d}rrf!-I1N zSaQiGb7m~w!1zku{>d{UExCQ!%`}s{Gv5DY$9Me6mfVsW%U=gPX7Q$%8!C?VnQzJ6 zKm6jhZd1lj^nEhJdc7s*yFau^^9d|&WvXyfWyM}guJ(#s-A-k&_}BZ)iO;v5vgER} zk0qSaX8hxPzo&h&Ejg22)VfhF`xQtG2VK>gnNP(LN486W6w}u7MBbkGU34X4ni(}uHtaTH7t&~ zV6~|GH8=5b&Tc$!U!VlfQLKAt;QYQq?nc^@?23krH)%KN@db_nlPwBo91km}cVQX+rLT-r&yePK*m@*H!eL zBjl{KXFX^?MGndbEw{S7cZraj7xzg0;aMJV%k>8jZC)#`&zbu?ZB-dRn`XMgeyfnX zuQ5Ad&@#5pLx&~nUHk44a_xMV``j^OaqWibV>aj?7IH67Z%;Y3fvxk9v?qOHz`4_NNY!#6R+}RkKR}=wJo?TSZY1 z5ONNYxv@!m+4}IqqPT!#cZA%yW>;mxEI+~-Ql&6 z^D6z^!gV}b&!K(mS4Sq43b~Jt)-ksy6u^DOXWz}~`$5Q6H1b*eY2qW`qAo){u9OS8 z#5%6Gxk+rD|GJ9LjqiREa+eKj){Nw_e;1&?#%jBXpM_liW-m_M<*|Q9U}CfM8R=id z$NwyRkRf*&j=!hTK;xvZLhh_>ch~N6CxLf6vbErim~XFbdHaD}>K}M$tx7JUzY9Ii zz9yB&0Zz|bTv2&m;Pg7L+f`mCaC%>z@O=TN_0SR518`bzhj6_Cr}exa*E4YXyJ+Iy z1vvd3m*C$KIQ`uh;NKlMZ6||qI{_X_+SM`Ku7F39cIb`UA@Jp--Rk0Y3!JudU);`t z)8}FfelCC~ljlesKS#hDk>~Cbe(r$N=hOs0r@-lREyB+=aM}-y!TkV7e|Oq%#NmDe zIPGWjaX$k*neui<_hcpK8sZ^r#R@Vmq=#9_Mt+>zLkJZwjRKOlBz9JV{a1Bjg(gY6V>V`A6L zuw4UQNbKMPYzKjVB6iad+fCqQ#Lk|?b{4o9vCDGUE(1SF?08pf$AOm+yMGzmec%Jg zI3WR#6M(NG9;X3+PR4a%cw7fOhl~T0@Hi0oEi!JbkH?L`w~%q> zZamHezK)Dbr{HlZ@X=%(I}nd!fv+ax-nn?(3w$XVC&%G&GVlU2u8zdxYT!j=93F(n z;lQ_&ak~Z{w*wzT#`%GGoDX~+nHMz1^8(-;nMb_F^9bOz$h@O3o_7ErLgpz?@jM0i zUNWz_hUYcF7n6BV0iFi|e@^C2eek>qxHXw)HNf*M;7!Q9Y%-pg0goi}xaN2s2Rwqz z`)=TQAMgWYp7;*W6M?&udF5n0uLQn~%tH<3^H$(D$vih2&vSt%k$G`9 zJTC^GN#@ay@jM#18ku+N;(0glnPi?`is$LTdy;wmIy|oju0{L+Mc5AjT%Gt0&SJj- z@F&F2a2)#?fG;P0iA?O50N#%HF-~GX25>Xt_h^Lu9>AlCpQJDLlK?+T{3<)KUj_I| z;)f~2ei-0I#BZ|@`)z>FBz~S%*w4c_@eAc*zYuU$;zzQiQlOJ`<;L%5kHkR z_EP~rPW)QM*slfrG4X?)#C|Z~4T<0E8upt3&nJGiP1w%{d^PdQ)x>@|;0K5wZyWaG z0nZ_RzvbBP2mBK86K=(RLf{p|ujq^Yiok7&A2J2|A%S-$eoIB{w*b%zKdNBrK7*zXP8ocPI8u%8@w8S$$JV!t|Y zUE+tIhW+rs9}~ZQE9|!iewg_A9k8Dtcvlh^IDz8=z@11OVK|N>0PjWO4tH_f0r+7O zrzpX33gCT6TtgAZHGr=oagcd94gx%$#7!3BxCwAq5@(r)<1E0>lDN!M9G3xpgT!%s zaU2JDCldE*fa5;ElceKBz->ufsXvY@0Ut!-P~&kN3iy2zx7vc^R={_YI9C>qa{^O4w!pWOIIk6s^8&Xeap8J6E)3jLI*tr{h;-Z;_-pAnHSkmt*H*)EZQ##H9NZtr z!GRAaaq~GiZVr4miL)!>I6Ls`BrdOl|acomT+dO*#(& zxXl2Zx9|byEdci=c@Cd&o&)d=(s>cUJ@cjVD1Zk@=Uo6lBb}!K+^dUJUI*}mCsKJJ zz&}gpjR3bkEtO{ieB(o@ycFPp%cSyHfVY#*djWo+Kq^lLc=0Z&yc*!GrSour7w?eD z+W~$~ALscb;yfSV6Q%QlfJc%%B7dAm1biyVJDP#>j)3Pq$9YOwI8O<9P3gQQ;7g_R zpn&(1&YJ@6B%NmkypZH&S>wDc;CH0+xPS+eysrqH_XYgEbeaNZqoPm-tC8t3T& zA42l_!f;+6@Ga7LfWUW3=M4hSmCiE+-bgwx5%^WQt zywj&R?-Y0#$x{u%d8)u&NnYz!oYxBcl5`#{a8r^uYmM_}fj=gBwl#2`E$|qUmz#(4 za)G}gdAwV39xw1)B=6S^=luc?BzeNOah@=6JCax20Ou70_ab@7RyYrtagw*Z7UwMk z4<&idyKtT}@L-Y`eG}(J0}m&8)K_sHHSpFX?^+4xT?4;B^0aeso;L76lGnWr=XC?0 zN%Fw+a2`1D(IjttGR_+ZexBr+f5Lg@z$cQt^s_iG9e5Lx$NmoIu>*ff^4?2u-aBv) zk|!^~dGf&DlDzs%oL3M00LjDG#CiC@^+?{n5zgBO?nd(b7vMa9;3h;D;D&Vpzxi!52G%tIFCaRIJy-_;JdNlk+GE`W@V!K5 zQ4i}ZfR7`(j15?q0lXd2aeTl!4&W__?qd+veE_dZbRtc$P6YTRqAPL1x)R{^i4LU# z>rjA?Ai9->ShoT^ljvMtV4Vx_mP8k0gLN^$t%;6C9qVX-ClK9DZ>+lk{*35!K4P5? za9yJ7(ZISM;8X_`jdei4qlj*(0oDxxHzYcv*;r=;{5{bn&BVGS;Gsmv^aSgefS)G1 zrvq5`1bidWN%>%%6!0*jtJ;ZmRlqr-!}^4ESipUVZYv(^wt&YHomW2Ac>!NXbYXq5 zE(~}v(UI}7jtqDt(Vcz9x-;NjM5p!~>(qd+BD%IkSl0$TgXrM0unrFRW1^dDh;?(o z2NIoKeXO$s-ka$1x?)`(@KB=T>y34Mz)OklZzb0K0l!3ag6pwP5O_zTD{PH*g}^5e z9b#RqLj-O_bc^;_w+LL9=p0{Tog?tmL>Fm?b&ps)=>hdy35O0cL}^1(P^e) zohI-tMAxZ|b)CSAi4L?3>p+1&AiB{otQ!UXgy>9nW1T7R0-{T8jdiKO?TL=n7wcGo zw;{UMNUVDW{)Xseld(<~_-mr8-Gp_uznR7kDJm z1$$y$FmOkrBX-3)V&GJFybtS+fm5AwIo2ryr@Ce>tZN2Nb%@WAAiDC!SXU0*is;bS zV;ws1VMMpS2J6;=-zPfvwOHp4T#M-92Vh-1aH^x%z&d*1RChlO>+XSbM5o^i>-2$_ z5?#Lz*7XB7Bty5+F`}ovf?S2Cu6C{2*C}VI5wF0XSW^VJ=>`0k|bu=OGlY^8lQ# z3t>0Ba$N}EbRCJVcpVAgblr*Cc-;x$be)Qqc%2I1bX|*Bysiasx(-GuUIzm>T{oj2 zUN-|cU1#GSUS|V1U6*4XUY7&-HL{LJ7+%K%I9>N61F!o5oURkn|8nIzA;7PZbww`Y zbwz;Fbx4Msu3U!%I9<17nnUHfCBW%ACtDv@u5$uB-gkxL;;%;NaqFEW6~etvyvmn7 z%J-f8uKyP!)Z)?h?)!zOfL}dxzgf=;BlJ;U|K`1e7lEt3RHzYMZiKiyuZ}IU%m#kN z(QoRRcjABFeDL|KUC(a&h|XU@^dtM zTxgj6#0cf>8((b9XZbwF+t;njyKjW#Pmi$niGKpe>9=rBlF3aYWEEBt{pJx{FJ<=X zO(sc~jnKt)9m;QaV(X<`9-y)XRdXnSFkDP0v=o>Iy} z^>HP;jL;B=3q0@LY#o(hEqle7ZZJabd+qzE<*|IHM(zAoPY~zVo-W$>w2j7FxXy>( zyL)IP86jt%8~OLzuzae9PS*;vVvW#)$s3~zLRtRR)AJs;&V(DGpwf`b?kQ}2m7D*M zuJew|`Ro69TUyfIbyb(GZ;~C?A+q;Kw(OOW5gEydGPAQcSs`RzCy|jI5kjGeGD9T4 z?#_9?kL&y%u7B>w{k+e8uB+?%e9k$q_xttwNV~ASpr0dI;ACCXJX!d-n%p=vzfGBc z_g=~NRG0U{_Z7SBK;4Ip9m&lLwO>xM5k9fyA`&F$d%UQ3@GV65Xl<>2;FDk8l_F-B6 zI-$wB;wY4FSYu@2-? z+rxE-?iRe_!dnkI&2S*4p{32|*ViF%J@?^Tb&Yv;dB z*A@Ol&wGA>dS(C4eXHI3s?S>~_!Rr8dp*i>;umH-&|T>w{Ed}b&7Bb$XHP~xThO~W zUU2Q0*b|CydqRATO_S;ipX8~B)BJV%+7q9;2aT^K3!d|FrIoU=J!vyF#lKny;kP^& za?>To&YpbE%i4eUf#7FiHukUg&5i^aWHv}1ApD%A*YodsrrQx6E$xEOxq_dV;u_iK zh#fhx*Xveri13B3{a`%fVYD5|Z}$4NNs-`5Id86S8fHgsSx;zLP?q0?-^JR$W6S}S z9Vz|uzn9=Qs#kfaXKzP*H(&K_K417(!xtUCz5kOfc{tfw|4;?tcP&gDzHwvO^ElqV zT(&_KEqJ$^#>1a%vn7?E8Yqg4g`c)X(Xv6Kr`eJgp|9^;TO#<=F(%1vJK7Qxf9Kby zO@uGEO8CW#Gu&;-=Cz%-CC3Qv(A;oO={Fk^R`=`fLOtR8o!hNtLdJC)@~g7*ruqJY z=Qns+Q*Vn6sj|Dnvu$$)ce>nhD>f3o<+<*|3~IizCN>iqr`9U_ckJPNeiSh2qHUrzIXgXj<;?`abH=4u z2hX!6CR6Q8-o6vOeXB=fjk{Ts3w4&Zi>M-e)UnFf4Xxd*NlpFq)MbMN_icUkLqvfU zG1xu*f@PfGMh{AjyPUBiZ5~+7Gffx#$CC{v^A}nX_tLhbXBZ2=_jtRyQ(U@Pk>NX@ z-t6r!c+vLkn{Lilq{gW;#(QT7?l#41=aM{2qIzn$aL7@?oAdt;|2upS@c)PTfW!O2 z`vdpjpASASaQJ)R?*k5hKYSm+;rm*_zdzvc{lfPR+>HMoFZk~R{5=2t;P(V>$A545 z{ec_v>jSPA;I;Yn1lJdExc=aJ1P<5hTYmil|IM%OIsAGD4)+7xAHcWq`^S^tPr%`R zd(7`Y-~s&pg!>h^6~CY1{s#V@-~R=CJplfguNSa>0EhJj)*ImLKTQa%Pr&c+_3I>G z&wv;3^=>s^|A52#2;rt5cS>Rpyc^A&Vz(?`(F`SoyC-Cz$oUehy`5VsTz|;A89nSB-;d~G0ec%Me_yWO)bjOkjOjPEeEo3!rk1a_ zVeI-;%hzW|c0H@*>p7X}H?@5IbN!~4@0Z8y{!+{LS1`LD)$;uq!R}vmCf9TDyWe5= zyIQ{Aud?+)Ew2xDY&}u0{YK^e%b)2t zD(`n(zfpOA9Ln}nD(|Pwn0}-3{_D;5YYKV}zBjJlsJy>({YK^UKo6$h=vl7kV4l)H zTL05;Ww~jl=YSt)`i-9AdJZ_(Zxr+#aDM)w^7*HP={IV_^&E`%X8MiF=PRz?sC*vd z=Qk>!-^|!~kM`ht4(3;}^C9);dJecD({B{?9B{7RDCjxhQ<;9FPF&9cFJ|Xk3VIGW zKMzyTbHMrenabzq8tlAHYjHgXo2a?sI1pqn7*T+ z=V1OXOb^mpk{;ymM?ufQc&^{5tTzprKBes>eaioS3VIIa;rfk&o&(PH8wEWFoa;9V zdJcF7)8iEM9B{7RDCjxhT)$CS?{ob||C8zi4y>L)LC?WFT)$D!bHMraP33w8*Kd^H z9|Xp8{YF900q6RSf}R7;^&6GzEnL4*xjw`1fBI3X=kWD{f}Vr<`T9aZ&jIK9jmq^U zzJAfWQay_6HwtYr4nt!Dj`@2?c}9L&S_ zV+wi>IN!f%fK-pI!0NYDd$MNzmY)wO=sB2&pC>5jIpF;KLFIZfKd;a;QvI2qZ|G*J zUd{Czjgacw{QN}adN|i_R4LWZ`T30Am+I~OJVzf$^?9z}s9ev_W%Yj=BGv!-`I6>J z{Q`a-rR}8t0zbdfRH+}q&%3mP)W6{8V>(IdckuHxJudY}`1zZjllmzuSpS7ymHID1 ztY1T$Nc|eF2dJ6U-{Ja!=1Tn_t~cmgsei=v32iR*o4B5#tEB!E*FW@@)X(C2iOT&i zuCHjc)Gyj@Fark-Q_pQnG9E^wSA&FfN zI3BK_!(7k7c(~s7vFi=T!}ZDgame9%Zprl=%mdfIn(H~>aKE%<_Y2Mg_ZRQSA?Nyy z_ThRC=HdE{&gOa!INa~IxSj(J>jUq{;rv{`QFE^6U_7ipk6Aql$8-HgLC?W>-jAc8 z=YYd{n8@`Ua9BURSUrpahxHcp9B^2lLC*o_`i*L_^&IDc^}hvM|B=If;m7qH%me!i z=sDnAztK{z=YYfhHH+0Fael7fXmhUTU_9)Ppyz;d{YH5|4(H+ejW*$W4#vZN?Zx#R zaM<4iSiKhKf%5?9IT+9N8wEWFTq~_yzm=V5n4Sa9^&15}2b}9S3VIGWoPR*i0q6RS zf}R5o=PS^2z~MXwdJZ_8-$2g+=lYF;o&(PH8wEWFoa;9VdJZ_(Zxr+#aIW9ryh?$? z`4;pXaIW9rJdESHeuMKfa<1Pf=sB2&>o+)`<9M#$DCjvD&-EJxJqMiYHwto`dmRzfsV0z`1?{y@vC&VEs4>dJe{O z{YF900q6RSj${2moQL=0DCjvD&-EJxJqO&7^&_bd>rdi5wORj?f}Vr%T)$D!bHL5H zzJAX1HO>RSAN-y$p6fRXdJecL>&MY$Oz-17yuV99&%yXQtRGB4&jGK+`o}bn_2Y0J z-jAcaABTJfuW#71`UY~if8c(C`781I3FtZC?RdR~uzCy5!}S{lJqP2teuH`rj)(PP zChO-TZ@}wC*IEA`d4=)`O{gWSFCpiC0?Pdb$Um@t90fgB{(slo}`-;drj! zpgx9t0I#Qko`ZP`Sw9ZyZ#cdKU*BQ9hw+X0e!%;2I9{p;A|J!+hoI+Ro@=ZhM?;x^ z1jh&P{S5RRj9zSL(Z;0IVNx2^p`dP>gHTz%4<23ta$k%B0*N|UJGV-vf@$j5h`O69U(?%NoJLF9@ z`+dj{aX%&Kx$^1C)~gAc{X`r;j_Xy>b1*)qHuF=0o&)|_v%iV++~fKg^c;+T&VDbI z{Ci3LR-ESq*Ylv~V4hU2|3S|I-=W$6#d-Sj`ocQquS5=hO3-sK&pBSd@M8W-9Ph*H z9WR-`61mjxMjoiyA4k4ev!9MUjn`{H&uNvve93uU-&xK4l{g;!l=6C#!TkWB=U|?( zyxs(Q4tOlDPl28TzJu4ZJea={-&gW8Aa~+^O3-sK&pGaw=*;|;I9~EsAm{pxfS!YS zT5>-n=sDnnH2oep&v)*p+`{XNFkbSL;CQaz243?K+kEFzkJC-USCd<&L?lUzXo)>=4!GnG$9cd{33?93bNz;T4mkK5xW5wTN#cG8&~q>z{1YF! zo&yek3(#}G!B5#(s_%iHgZnG-eYt)kpyyy7uHOjgIpAEs5zuqMxqc&yrFtg#SGd0t z-#4D;7eLR!JY2sK&~w1SPXl@mIM;6k^c-;T>!ff!2fPN)Yk-~u4t~nUQau^`L)>49 zb!_k(ZDZ#h>LM{iOZ^_{pwI{Rr@%^_Th=;8)}RN_=1Nw_T9>BjAV2k@_j%pSvsd zU$}lFEu?-8*8}92)ZYO=-+QSa1pYrWuIJ!;0KXv5f8pon`i*Rt`cvRXwBmXW<^lhr z1J`rF!S5K)^&D{UN1l-SYruI94n7Zsd3X*EIlM3L10sjd$#ZbX!8gfsaLE7hJy7`G zcn%K7!}rX6tjOVa0Uss*dnx!Pxo*Sp@VkSKhw*Toa33y?hwBP_lrSEyL#gjZ;kxCz z4d;RDocoB8!+ilhN|+z+BdPC3;l6`93XF&Ql=sbXez>oB4h}i215)3O!n(nA8;*x{ z27Hw8zTlhWzHb~a<=~M2W8FjFB={&{e(+6lp92LB>#)>!L*FFVZ8-iP>pc1Pq z*!QKr8-?5j*KIftoGYMC4D&#)gXcVOJoqN1z8lo3AqR)!A!ovUO31-C$$d@8!8ZwY zbglB2FO_m|IQ}2!ItsZO?#sgQ;G5(=F693>XHv-daNUOE!8gf$Wyrxd2|h}#@|Q3A zCc#Gu{2%9P3b`imQNnmQw@ZCD^i4|Vd@AMOD9rN@9f7_{@KM5e@J)h`68JxK4fL@f z2Z!$qzDekV!94%aS?HVOx((+6-=w7Dpl?T@j|cOEZxZ@`!2h8up^r$QuL$G8Hwk@6 z;NY8-`fli(l=^Pyn*<*v%mcnjsqcoqNvZEf|Do%tl!K!%Klmo4x&icw3HT^syp)5( zdHzwyK;IX7swbxSJc;P8FHHwiu$m>+zTQXQ2_IXIl>A9Y&vO@dDc<_F)TR0pR2s2fu$2Z!$~ z<=~J@IXL7}4i34LgF`Om;E+o>IOI|e4!M+rL;jDx0F`oZI9|%ZA(wJ+$fX<{aw!Lg zT*|>AmvV5(|IwGBQVtHsOF1~?QVtHel!HSq<=~J@IXL7}4i34LgF`Om;E+o>IOI|e z4!M+rLoVgukV`o@m9Q2#iccW4c4#$Ij<9ZIcl!HSK`c3Nd zQYi0U44#!J5IOI|e4!M+rLoVgukV`o@U2`yjY>H<91r?U-e3OV0|DPZ z=sBE6%E2KA{U+~^QVtHs|HJn}r5qfN2mL1R@1Wo0^8olJ!AGf8{_;iNB;+=LgKrXi zl)yp1Nqsjs|3D6o!g%maf{zk7_$KA^73eqlJO;i=@KM4%;G2~CZs?l?A0>@J)h`5;*8L@KFK>{U)EE!8Zv$N*E9N4SbZqLBC0TH}p-)=YQ}` zf{zmB0pFz5cSGMK_$XmK=r>uvKwnDMJ5mmg!aU%cl=T$!$z=T{<=}80@J+)14)cR= z5`2`v;r(R&2);@9yf7Yoli;HS4u3y0@r4h}i!H>vLi zdSCY0qHmJt;BX%BO~UmE?+f}3u3zAAearOx*13l5%i39(>?D2ZtQ=8`pElr5qe`SWo5p82Bc64i3k|dd+ii z$YFiwIXL9ttC#D0QVtHsgKv`O;E==q#B*@S!8a+_Cm}~5*E6LY9L@v#E7x<#LBGlM zRoK6I4i3kIZ<6QWki+?a=irdTc|xuaLykkPCqwRo=ism|4d)y1QEHXHd@1-Qc@7T8 z!}$r$Q!pNUPCN&P<00q5b8yJPN5yk+$iX)$*Yl+u9FB*a4$r|Mhw~`U!6Ann5YNFO z2cH?w!6Ap75zoOPhx0Vg!666V8_&TZ2cI0z!666VB+tPi2R*=ZaLB>8$8&JVA-BbI zaLBjqb8yJv^Y9!Ta`^l_2ZtQ~UY>(P4&M*Y!6Aq5ljq=&!}rf~aLD2J;yF0v@cZ%{ z9CElGcn%IZTt7Sqha9dqo`XXU*QfMcO6bQ?o`b{jaQ*Wf9CEl{cn%IZ++WgjDWM++ zIXE0I<=~LR{myf6$YFhuo=Zuk92|~^^@r!+kb{2XIXL9to0OhQNnt(YIXD~->nG2_ z5x5Rvz2$liIjqlI&mos`aL8f(2OlNO5Bmkr!Qpt=U$~w_F6H2m!~Vr{aLA!v2|h}# z@|Q0>$Aoe{hvTIj9C9fKhaC26o`XXU`#bn3waQ<o~>2>@`dM^Kn@Pa!}$t)lrSFj8`pC<9`qa3<6%7LH?HS!Jm@!` zgF_Db4SbX^59l}U`@->{-=ya@z;jF}_$Xl>&~MUnA)p@zIXIjj^c(mnVII(L(sL@{ zIVO;U!+Aiz@qQd~&~IGNAqV{iK1z6B&~MOR0xsznoCow9*K^3F92|1cQ{bb7`9Z&N zJ%{6=-v<3R7!STl@KFK>{l@hi&I9@le3UR=%E94y&~IGNAqU?i_$Xl>DF=t+LBB!& z6UKvY5`2`vp`QvqO5mX1pdSkye3Rg#1P;DQ==TB#-z4}bfrD=ne3Zb!HwivU;NY7C zA0=?`O+tSgIQS;PM+qE!li;HS4*CuH?Z835fsYb6)PKN72^@Tr;G+Z%zDe*=0+(`d z1UUF6!AA*P%E94y&~M;(fbrm)1b+l@DF=u1fNzqo_sFFj9CE21h#Y*A;G@(kfB91I zO@faSIQS;PM+scY!4cq6{S!I(Ci#Ag9DI}Dql9^+dMu6y-z4}bVZ4-s!|~voAmvV5(r5qe`DF=sK z%E2L*`mM;N92{~f2ZvnB!6BD&aLA<`9C9fKhg{0RA(#5y$fX<{aw!LgT*|>AmvV5( zr5qe`DF=sK%E2L*a&X9{92|1V&wyOY!6BD&aL6Tp1#-}D;IGswfB8}=2Z!UO92{~f z2ZvnplOPBE2L4K|@|Q1_a&S0a^0y$Da&X8c{|s`_Z#)NwT*|>Am;5})r5qe`DF=sK z%E2L*a&XAC%3tVTLJsIOLLF z5jp5Lo`XX!<=~J@IXL8?-?*MbF8MQ&gMI^lrB?aNmr6M}954Ark%NBYdJegigF_Db zjq5q&px?NjLk{{)*5jbxxSqrDpx?NjLk{{)`gbefIVO;U!||ZsxSm4}`VIV*TIDZa z3i^%fIUEoAO|FN4e&c!$$Af<3dJZ}0H?HT9gMQA2mJ>AO0DvjFXmZz4i3jd{)Ok@kV9UE=irb-zJ}-Ekb`~$e{;Ipn!`4i34L|3WV1;E+SUjOXBxLmrLi;E+RpjpyKyL*9+&;E;Pv zueLw!kAhfEIJ3?^@q`+_PsdjGc6*ia302K!K0baK`Tp6h*yf&uhju}kQ?N7 zia%4LAhx=Rm7g7ZfV};MVf5us1=+f((vN*EPmo(~UJ$VBhk}eY@}OUw_wKl~c`mDt1zp&t}v)sFT)4TrozZsjxn&B}KQ zQhi3NS4B4OkjMG9i*1vqAUPF%Ol+#{CKL8i9x4b;zEFGhhtJ5THu_V;@}7e9^0&Ac*ZM2+mwpFETuf7tSqrWV zuer4d`K7`sRg5kxNXtzxIyU?E4SCM6&}I4O6=Y~;w}_+lfBZy#$$ISw#U2HDzC!<=UdIyTRzpo{4BxCEnN9ol_uBp& zdGKYwnJrf$rDp=L$G}kwxU><;8BKW!MXkHD3g87*bYH&Pn^dQ|(Tg0KEO>e88WZshYFub)Pa)WbJ*@J15I#zY$K=Fcc#0M(LGxjy(ksD zrKRUP<;b%4?H_YHY_lnj_g&a4%CVOlc^0qJcUTiMnK-pw{NM?TDbST$ngM*QvPXB+hryxx#*h8Nzu zl8U-b9yZu1cxd5#+ZLBx$-9Z`jP*VW{?+R4DASFuWS?b;dYoTnoPSf_M}=+^UCGgr zh7Cg(3GP~JRNSlpS2FVHj?)(R1#jAZ`cO|tS5oVa&eC?S7C28reYdtYA6&?@^=C@# zMhjkR%*r-BF1nDEkmQ~-69o^6wHfhzwF?>Dc~;(}Qo&blsCX!8xC@zm_lEAT)|NQ` zkt4UNWH)djPWk44ELI9`UHjInAVU{&VP}2s)0u+X4hi|?@yMCfN~Is)DXegweLYO) zk-g4j;`qCp_Ky<$Mc~hM2~(WOp3QYveM}I1+qUt2ruaJ(&;3K|slN-}FEOTYxw$iO zII=ml#I_Hc zRK?5&=Xn#l^XiypPQ>Wm?}hWa3Es`5wtXu@C-Nh%@cMwYf>)2~IW+vPBiXU$;PI20 zf?Iq$^QKaqBk^0e=JaAGTb%z_!+mdbgUkMXQex6OV}Rfp2TbF8S9c^4cOy(UZWG)} zJ$-7Wj}D~D`{{i%vIVaZUU_5nlMckU<>u6d?sho8=Z8fj_k}xBP)XUtSoU{oyI`UBPQ3hupuDWYSVnpp?ef;V+NeR+h19eL%@@yO9j zg3o##K6~|DTQcKB?T87bf-kfGdaV8mTcXg}*fX|)Bfjt9jV+&d?qo|gzU}eSdV=6b z>K=;BHnt_lZXSP8yjSpAPv^Q8T(u!_S_y_tp9}u#m!iq!1vbRC>RbC~*7$rj+I8uN z{$b5+NXUXi8Bf{^p5hx|`K!d5{Hm|5UvIA9^;%7O`}>qNap|%<%o7()9;Hx_t43E2FNlxs1>e|7| z1>ZNM$sgX{hnp35nJ9v(hyr z!Z*1x(Muk^ut8VB{VKf<3@`h4E4#j5+2l~T;4eBi+S{*sWs>CbE_2Lo!H2rrR)74| zylmY&7BcCU;CE_|k1JYkPDUAZH>TeNr$fe`NoZkChP647Q@xTKex5r=(pNj=S0Xc- zjO@GCUGQ_(cO4x%RU#LiZkFzAMmHZ?x`D@cy1gWWDO;yAxFK{KS6QDaR|4 zeB1WxhTReTaPy9NZ7NqJE4(+&IR0GlAIncE{Xz^$X6+2Wt%ZW8CEe|^{dxuRXF%Pk zYQF^^TG_;HW|a!$_Q@LhosAXvJ-1r+tK*%y2Bg`ki9;>y1;1$7a{czF`lO_DX^%-& z1b>wnss2)5pL{rCxHzP);1<~pqOPydBaszNPZ$uvFGU{bzL`o(GD%*cIWfq?(+l>o*8s2{qrByoGdf9YD)#5wkj>_>DE%!yN^M0 z)~*$NMM3iYx*dM2EPDQ0vt_H`#jO^Me)0X6>U_0v|B-!y_w#DnCx6FJRr3|I?86fT zuVKWXXhlr*LRu{Q1AP9Ro2S>?&s8k&%9AFZNav0s#D9Sy-7|L{E+3~ z)7{z@sUEyf?KtSB;N5b1XYKs-Mb-T4&6?}d1%L5$i~igVpH;>>qXrDm6kMz7Jl)Ua zlPbLD?Ub9(1izEHrswb5A5?Q5tn$10QgEa5`aAnO7pSIP*?nN>8^O0V-nrFkbiS%# zSl4f>-V3hsTm7fb$vjn;heWUMN5NA|!)IL9eyv)sH&cE3i{JqVT+^O)$WayC9%YsA zP4Hu9tOq#$_gq!2=ABPXiv`d18F?cmB~#UF+5I*{ehKb7y2Y@J!22qbP9-L8zXhLv z!)ws%jGL;8uWI>DE*0D;$SP#Q?~AH;)eoK=_($+=eKPGFS|_MH7rHIa`HN3$Jhe;h z#5HO=;0GWol3;#S>L`@oB7k9s+TVOdbs~j!HxQyS=Qq73w7&z z6)p9C2wu&nR()uAr-kzi z1^+&#MVFf?@6^jn90M2U3$AUkV0+5cLiL;&lkdu0!QXt@Gj#sykLrgNcF*aPE%)ZvOQ*9E_DylPls`|s+)XuS{ZQv{!Mx6m(W#}Bo6+k_~KbAl&W^;G-*C{~BQ zN&G!AQShcG+BKQp@t1mX;J8_%4+`#D*XyIpjuLg}b}{uzb_(u2w%E<(+i&%Tnl^z| zHwwOHQ;600wtv*te@0oqUnzLXt{q1DowTTVgZg{gEfRcCbj{OWzG=~FGiNsUpDnn7 z-*EMut=e?$h~rUo8T8S93L!dsz-lZwWw1&KyddOKf4va)T8cI4tZW_BzX8i=N^w2>r>zD zLGMax2_8CmlV#ig4Ct&u({JB$6@04xu-o|Z|P{O8|S8`6F0J$Kj36MR_fk&o5~SEL`Fth!?NNN}%F z)s|e-Gok~MzkJ%(i%fV)JX=Qth=&pk2R%|=?!5=dkoO=DYp-S*~F%P@d*;I)(n0_IodM&{p zUsomNJDSr4pHJL-XD4{r< z3~tvYsfpm`W4byT`dCo!+rHQ8m)NEQ zI&#s9st*3VQD>Ck(c{8i46JBPee}D$NvkjTPyZtsYdTxgeF=pw=X3?Xb;K40q!Oo3>D+&oIG&$^STS}^KRwcZwa=vgMTeoYe&Hib_YaP`fN)ln77=qD_iKz&$&0W zH~ZMpW8((iUAI^80fFb~=1F#R+fS>~4`T&i`y;2D`7t}X^m@ds!*v7?-0a`_%SSu< zY?t4@`$a-WuWeV&yVZa8)Syo9$3spF-fp3Cr`|Yw8fexvwEs-O%~p9d@{YHs6(1R8 zTx=%z^!3A+rRLhxCe~L*uF(?w-l?@MQr#VBytVV1jAT(CP?}!-QG0*`RiF8EB`Q+z zftw-&OV&Hk9ks_#jBg?MkQ&zwz0)12OJs>2(GfgDyULj>CXSQ_)GKbAEb1GPJNG$V zY~@IM>ttWKHb-#F;P0=TXFJkATk;ENQ^Bp?c1_DX;Ygof41DMPOVnXX4sVG0`reTi z+7N%oB*Bfk1-%~Y;Y2gSFNM#YEcl1#X&ajdIZ@l%1}?+A|Kb&1O_{jNiGIE3eDcaW zQAZlzpnJ1l7oF(UNQ?CNJ%ZD1D^jw)InlVocMQx23;w$1i-)!~oasqjMeYv=!FzQL z>bkF=GmUS!)o{dpQSX}htf1nuWzO{R^#sFyO9el(bHej0YG?Z3K={1OR)UurZ#v%N zqcc7J?5SP$Z&8lmnz4dUXffG0XO;^c zn&kMrrHA01>df~@IN(CNT%O)-L8ho5ZW+}5MA{=4Iy09h)3$ zqHZa8%T+%P zGky5A3rAe(vCf;iwz3gCa!>F2>mIq%C5K{^@z+E>wtAsO!@gQ>)apZY??W>MpV9by zs;QM5J*m?zzRQ1tFB`n&esEPcnv}WZbj=h|Hy+ld?Wdr|Zq&&&XTMr!>Tk&$Kd zT+{aHaYWR^zjc_pX7egH8a;8*BXuRgyQSN<>%Gg3_P+1cY0-F5Z+|i3ZvA74ZZs(N zX#JzNMBP5@eN0LGB{!N;eO#kvwFRG{8rH1EeK%VFR2%KOF{1weCGgnxX)oRAy=~k4 zR=p8@gttAkZ?*gEM)!Hu*&N(T@Pl(4$DP$u(0OLsk9KVr{fK7GE7ZJZrl7Acwcp+U zw*l((14c$CL!A^f*mx{m(Od8?-tiaCR9Dbsr_{XW$3@@7EGgmp=K2b{{=?CylT8Kh z;ZW1Y*$15n^z4g)29lvm%b`w8coPMnA`#c%hV|7Q-k6U{~ zZRUGMK||uhGPAymzTU%Cms><7E2!I-(>2y^|BmyxtQ+MUdR0La#-^EFZYTPK_p6LK z`~8-JPOj3m^Q14L&zSx@=cv^K1ueA8PTIZk3(j-Hx#ibukIUw{{=i{`zvx$*ZRt1I z_Jx8rcqO=i7HD$~iUpUwWX#w`tbn-MvF(Pg@N=xV3we*!4-lzy%9KbP?r^-m6{ zp1Aw}@tyU`xeBfyn8%9EgB;$s6?(X7eE)lqXKJs?Gi|vERVc+y&`zCVOXYaCo7CG$8J=wmDya(UMd$D~S`47JDFJ${Z zayTcnVCMwnaIUap=L+O-4q3s@A;{s}(wv=Jki$90mYs8uFX!i?Syd+jXaW{!?&<=IC3o$c5Xk)&h5z0 z@pFDUJLe+@UGR|U0_310_AwoS9CXKVraO>djl2!l>3x_^N4}oxdOfD=kq7ZQ zz*ANSK)#;W4Psc`0Qn$ZXV}8(49LIoxJ$3X7O>mKb`-2?eW%{mG4 zhMIL1@^L9m$2&k&xfvb*HJU?u7g|uT$+|bt>db zH0xT(H}N`HQ&tBmq-|>%a?H9T<5ouNyaHbz|f?nssL6gEZ^X$fY_qa;ff(T&j~J zm+I=s4K(ZU$o+ZUK9kk$k^5-Y`H_$2eSt#O7eMZ-*+)P=iuWDHu)YKG&Ad-h%=#3_ z=kUIU8S85xpUL|mGg%)5`Fzd33Gy16eHP?5cwgo_>&qZ_((L0PPtff9ApgYsM0%`G zgj`3nuY|mxW*-VU<$bHltZ#*UiDsV*c{Sb_i(`E;VC($iHd!F_C}N?0X`&)a;WYHy*3eS4D26*@s2$eo&)ti@cL& zpBK3|?+XX8zA*9}%|0^neVTn|Xao$h|au zN5~DhPbr%Dl#ollCghS23c2K)LN58NkW0QSPB{?uIbm9rA75r>Dn!ddQ0>Y54k(U*kT& zB<2G|p2&TJDa<#B{4)0$hBKcb@-y6*IGp(sk*90=7?F?QzQ?o7_lP`O(L>_ z)%0N^57qQ-BA0xg$Rjj;p~(N!^pPT$e5c6Aa-V8<=2Jx;%YCh2Cy!k6)guQV{srd4M=tsHkyqnB z|N6}5j~sFV)mSb7`F)-vh-5heW$5hC2OvoYkWXy6;$RQ^cz;aT^Ay?Ib z<*JZF4yy~xVId#Sb6Z_lZVNf&yk4`M7jnpj6|-Cz@|8SC7RqvD$kjY|=F4(t$oKG^ z+8CBoLtc;P+V-$q8}d~=2Up5+aLAYO+*~os%^`=JT`tSnA#cQUd9Ey%ha7TzEm)2Z zIpqFkvfLkX$O$%LIYH!7t`NDDLqra_#i=Z}h#Yc`B`oKN{4~!+7O`9;a>!A(V>wFX zmwE1T8_QiHhn(gLmeWKIxz4F9*NGf*phH;>6glKZ6)ZQ3eJ)F$Gkwc)rpO_e`ibRI zkwcDk2FtM`huo`yYkBS!Ipk!+SWXr>U3Xcl?m$j*&x7*^T9tkwdO|8Ot>zha7ZumV-tv<))EC&N_tU ztdT=5JD%mTkwcEVD$8*rhun7n%Y7qv%J|3nj|Z*SzKT50W8~CorAl%jrSML^mgwWBAFrxDU7{q8o2?R# zYl*)938$hCIX{);X>&*;;F9K(!5yY z%9|aXbe`lX$&vHGf5qi{EmtnhQIcWZU26pwi(L8ap0}$^$|@UgJhXeaX~Msu zmpp83;X@_SI{H5O^;hBda6Q)kr_MbkX=3g@z{NlM2@m_IHQev&Fi z7b}0gq$Iw*Yp2v2Bl6>8A{OZ0J+CAKy6z4($P<2*q)ESr1)Nlp&<*Rd1_le?%Y(<( z{d*i%l1jI%l|ES_KVI5tui3BNO42*#M6IcT!dDa5aq5o7o0VkfpWahqM+zU#-8nuj zELSTDJ!x>zQ7L>o#h<6xFNju>Exm#ocFz_(ug#^NDtzH0Q{p022%BumYY z&q;YExKU_gi+evj$=rq;AKR`K{;@M_?%r_8_axm5XT}zE5I(f;I)6@8c;ZP0FPcJ{ z>k8jm(QM}qOVd2bxl5)Il}-z8R+4+~$XQR4@u}^H9wEXP7u}>_-S`8ZWQ_3{FEbC} zuRAns{hJ$co}_8o;_LdtFi&FkHDPm}i}3qh znfJ5#hT)#1sD_uLW0>H94oz&{r%{>XZa>sT}b>TmpWUOWL z%gd9r9=f6Jqq%~|-%Mj{%rE*Peq9#6I^RTMX-`$37SQ%jhjXu3z^c^UU9+c-9xb&1Q?D{_Hy8K}G}w z`Se{Pc)b&kJQdqKNH?v=^L%pz&#Pj3TNUj=eC;D$k2Dbe(!9X_!?h=Skb_IAhg@4C zc((2N9>e>4kYPrlkLP9!{ysKSe|&%k`5KXaC!?P5yT*61&2#+EgXnDApiW&ZxNprB zBgWZ!kUd?CY`SC#p4YVUk0F2DNm`9RZ6^NZ%XPV7XRMXyP8Qu$U#m7x_;s(i-nVLT z%bipk)S%Lo2ZGxTa*nBS!kxT5GIV8+s=`0)GHCki(;MAMt-~FzUk(#IVwm;jk8|CL z*#_9TjFG@N3j>(CEDia?@$Zrs2(npL}a(&C?soo~N-o;k{OaWrFu? z)OG%paSHNw>9@Q^j|C51AKsy~vx0oKD1P0iitxREu+iV*>!TocIX71xm?-$3t#=IU ztrVnF+{yXV&j~)eWow&*MQ-Hrydiy341{mL->RPd3p3ovj5@h7EHX}#y{nM_fGQ!+Q_ow!zd&+#)TIDZadgW;ET2Fer zk%TSJQUfLlKBncACb13O$n*b9oAo^{_;=&i1|_y`WMQPGTh1TBogS&X`xLs8m;!C* z$E`%pV)`TN$omIeiO1NLdCg-4-?cIS{_~NpB;;=OdUkgNx2s<>rl-=Cn0&}P-`7gy zMBWcrnp;2Hg?v4{)6clK;IW~LYcAO2LfW}~3sG(qynaZJ?ac#Sh|V;X_5LiuRkt%& zj;riKYRu_1{rX?IhHYVIWL2$oji<4glJCR|k zw}X6c3*IvOly-A_CsOpry`J7*Ip2uXbEEBUmpxY@Zp^R+CJOwVNm}>&kD2dCJiAXF zbfSjf$5tgYYtqP(sItS>u4*Ir`v+4(7Z*5?gcV-XUCaEu<(Gq&FWsA(Ff?$x1G#_u zad_`}g2(pmI%`i42lD%7d&_3q1UEXmvqo422QuFx>%{H9a(2f3+CIH`+Me_?neoIV zQ;b(m@zL8q#-6l~+xGiwna>p7H|X+O8)IjCvOH$Ux+rsZF%1bo` zzfk->W7KRrViKrbB|bp#k~UMrMtRwht6Q|}8x0n`_u7r`in45p%cr+0ll{Z; z(i>J6JmYP4kEiXdNqVnd4|Wj21FY}O*!Ine=-AA6KN=|bRBw{1wZn?sUbCg+kPyK` zRy4m}+QEuE*b$wdG+*#Z@oL9S-z~|Z+}F#F#tGg(W#SnyF+UTsHju^^A;CIsde3ZC^e$=a{A1?h0X z-`zyl6TgQaOTXVe`Jpn|pS>*ZmW|-EPo%j#U00br_X`cKQ$z4K>4&x)QdK4f?>0F5 zHWl1$)$XENdFG`0;wcj|ItzZMsQI0WzQw3{n4Y#UcoOFbXc%r zqbcciyNB-nQ-TLu)x0|{*Mw;Kc=~R=BKY!&t0pIVn~;r)m+ND=F zRu(+-L9MaZHW`xU^LyqtaTUCw&(x3_uPTuCMHfa)t0nlGr8>LpYgQnyb?cgSY9#pa z(>-R#O*9})Hbieyw-Wqt?cWm;PU(}Vsu2&M!`XVUwrr zZ=pvTHr_n`_(;Js8d}X<6{SnYnjaXrEL8AeLt|#R-OwSK(MosS*@8QFJN2QvxehVf z)@(uLMS{oVx%R2kLz_&`trT>6rQly)lTy|6_0vzbJ}CIt(7VBkJN{5rcYM*+=(^zT(pEq0a`&5R+v&vn)9wg9pzi+m zO?nim3RT0$&&m+|j`qN+G~G%GcE8BFr2RwiOKz|BT%Yq?wLaa;s^3q+ zCnW?`TN9C~y0`q^@s1^eTjp)|&W%V{656Y+ju6r zq0R+W$ms`)!+!)HoA-UpOYQZlkO}dpOa9^`Qo4;9qJ2!=>85T}&R=}8H<0#8Y*0?X^D# z{}P|4!lI-p>r~n5s)O%q`xOgboUPscT(ehds|Cpe4t*0m`1kSTh^Be!X4hk!Pks^H zDtX?lj3)W&l4IRF^!h0H+kW+09%)pd{&BnJuchw=ujTM@^lZnE>fNDse*NDFt{S*u z#-k0N)T8@V^0@R;@Ex_D9Dio}MIGDo`{j$z1ka7C_-pZ+uj;S2_WA{73Ldwwdi`{Z zZ|ZSl-fW9X7yN3PWyj;ozpE=QoM6`Jrr>)@S6BBiDOT$r&fRe^Rq*;RR-AHN^i#dU z{>7}FYQc~7_ISNUw?sX-(hKU7B>2pLi#wuc|5lr|4qV6a{W!Kr(GDnR&dvz?Y%l4)23Bd{__94MDUH3 z$9lC2(V;c$_qp576FjMD*P*Yg>(WztDT^P33BJfVpj(aCx->YxN__1Q!K>0XGtIZ@ z(Yi6KPr42fe8PH%k*oXY)8~T{Umge)e6wlS&Bv?^Xouo&3aumPKc^NYYeK_=Yx;nXL`Qxzla#px(V}|D} zBI$uKO)0EV-xp&-8+obfPPr}kpHBg6&6}Ijcl-KGO;8KIt6P{)L6IrF z+i}OMMMng`^|92hNxT`2ziF6Ox<&B5kMuu$>tBh^uGK5dVwvDI-!?eBtn7Jt1KL$- zdTFNM3-*6m@L;ex9X=#)jMW&yN$qJpZ=ErxHKHE%)(sN;p2zpe(=02~J{1i6tPBu5 z<(5aY+hZ!zw^?1UT&OEJ$?`V2pHi7FG8*tR#!c|aB&*4MM+>?(Z%e0Q6TxGD%&7Tn zk_BBowX>o^u{e*lef+KQi)$9taM#4UCteADtiy}6MQ)b#bL0)R{vE;F|C!J{W2z;6 zIw!x=_cMYwT)$fD$_-08DaQ6?^lrgh^!YeG$lZ!w^c$E%mkWOUV*dA-X;!q!yXeFX z(*=L#a$KkX4J+zLT~n-w2ySZMw&c2-HI1VErq}Br_y@oF-R@7erbBKvY+39pxJ#(x zlxbJ3>DtYMVh1}5em?ZzgJX_1wAo?b`yu***YHlsn;&9B=hnOut|$=a_C9s9-QH6h zs=Aw5GX1{bkAKYnk!xW~_oY0!5q4JaxsB_*2_0%nUk%I)cG)5LVBg7OR-Uk>m-T3a zF^dE@cwYFYsev7hOuv13V2I$^uKko@z3k}0k;=s--35>Ru&#fzJ$95Vo%pADBfwKJQ)%wc;feZaW4lr(u0-;=MVq&9Fo;rGW*@z??_E2 zZHe9MD8@h7>(IC9S4VoF(aZc!e?9-o?L>r+N3*g#d(fa>yPpYtow3t&USyOLz0%NQ z!~wP7N%k#21eHC1;lr;}34{N7u2-bKm(xBwXIlGOm8(}H#Q2kzLynFa~_E_B1s=~ZW(6MRj*DS;m^xX{w;LoE+X7yP8Tc0zAmS9;R@yhT8L!9x;bJJ)XQ zN(Z_|4e#(-)LR^SS$r8A?OOJnikNQdLxMNDTf2|VO;>uoVQ1%UBL!dI>4?Yb3U2gn zO$+mn?ts+qluR#oF6))(Bqdc!%(uk#6+kN#m@j zPJ(y-eQo&3MP<3^y(x3I8wkG0zmKF4a{|;?-#6Xf3l|FQTGxO-j_Xhri=)P2_JpL4YYZ>zXp!D6O@M*F{cJU&s>2~WH;iZa=t zpuaCiMLioK_?VNs>mEO$pwHAD$K7@iJmhrE1^V|Cv`K~5zx;2Cx~12z6JK=-6?9?J zogGc*2wpY)Q%QV9clux8%7^E@1%Er#!{N85J00+HVx`o-o@2SjvfA|@DtG$f)!K`J zv7+vpP!P6cb02s5XiU(U$$vekVpQeYIr)>_=|_FC(7Eh+9$MuuUs|zZNyit<-07hB z?<;$37yN=N=`m}cJ1vguP|2~2;L8sGn30_9PJNnq%$-y!>d%u0-)MXAu{#}^HDPSO zgM#m!Y}~Ht7k6s4(6{&FAi>|hJ!&x7*n^JUJ^p5af#9!Bj56us;X%W$xCE|0F6!j# zUN3ojr?Cf(YU~sFWT4=i&sO#N*42ZqJGd|_$yo5omw%2?kMy9!)mg=LPKmnx_BS)@ z*P7!&qgQ2&Y7;E@-Vs^uerr5vQoG#b-zI`Tc&2}~Fy4cD&F!+j`zg^!c+hau1%K*6 zS6y&j(0;Js`Eh#N8$a-%6P`J{=NSq9r_B)K^tWa6Ot@iG?XTz7lDS(pnv{A_-AetR zR_ZVM9V0pq`jlJQlV%*rb8J zo`T=0V*ajvdrvxa#L06#evAIj;U^BKR`mCzG0w{uuiGW~W~~!*7mxR(iRpFCW7-Se zuh759ZLTL>Ut!~?!LLNWsc=`x`ChS}^x2t$OzS#=C--^1xW-OT8gydhfAwaGzLksa zL&HtSJ!z#)kJ_bN75%SoX|vzQr+Cuz^z%_EmV*C`9r&hKx+kqR&#z^p0iutV+co{@ z)*Me-5M`2Gb)V?RIYgRhE%@X~y$6j-j{PM1cY~(cE-}zj(x?hL?PoO>{Mpeb=ALFs z8h!Q6>-x{|KQ4g1y4Lh#1@XFXNaSJEVp z-CaipiN566(X&n-Y@wt%EA%GlZ8MLM@bv24-!a3!sr zlz379mk#piapsHGEl|?ywj%;Ms{~grv(XD$uB6p0mp*SkPaDS%*;wt?@%2i&GUBCa zcB&R~&z+}E@7-SZeJ)At=4~W+YJ$rq-*_c`X0~kb@Hp|DhZgI1cIb0VNt*?2?ETnP z^tEdyP1|L6Mo9}V9bA*LOFZWx`dHHJ@ySZ+^=49)>eWTxe4*3CqybmUKL6s%juwZ+ za~?i>Zt%{%rKCekelJ+yE&B3tb!@NdKTy(!O`@4|C_u$XI2juX*oo3$~a`>K4u!n<9hZxB8T6775m+h!*w!( zT_?y#^6Tn6yRMLr=htB%yAF|W}U4{@)&*}nX>x`c`bh5 z-DLM2a=1^O*?o!}?rSx>uaUz#FpaGP1b%l|HzL@&fgIKuJGRaskKyam8MZDV|G?KV zZMKdfkLBy04O{n+d-HX24_haZ!@8Qv)>Y)N4$o!lFmhP8@3D0oc@w_Q?`7*e@+W*> zh+z8yazDP0ykPqX^5=ZtnZfoQHI404$esAU=F0XpA?0) z>{Stlr;jdjHeSyw1ug>}dnS%)Y;fOShdS+^+Pi*-&{ zS?4I9iFHviSr;kqfpyeASw|_ahIQA2Ov&pmscomN!VY058SUDsXKb;^%m9XL_e zfy#fyy0Mt78*SHLPF8*n z>*_wTu2y~t>+sI94p+Vu>-J2tZdX1S>-VlF|7by3rBd$svp*%h6j#g54 zDDQzf<$}~H%D+Zkb3p1EPqD^QHSQ0I#l^>)U6GrZdLvZ>f8ZR=PI9w zy11^?#maZ1jy@xGwDJt7yK_q2t$Y&d^czyAD}NPr{XD7bm1lzva7lCk=4~R z`B~@;n?z?&{yB7s-J(k