From 58117628b8ae8384278c5ef5ef5d32e7c922fb4e Mon Sep 17 00:00:00 2001 From: Andrew Johnson Date: Fri, 26 Jul 2019 09:53:14 -0500 Subject: [PATCH] Fully remove openmc.deplete.cf4 function from API --- openmc/deplete/integrator/cf4.py | 133 ------------------------------- 1 file changed, 133 deletions(-) diff --git a/openmc/deplete/integrator/cf4.py b/openmc/deplete/integrator/cf4.py index a614a42986..7e0b603ae1 100644 --- a/openmc/deplete/integrator/cf4.py +++ b/openmc/deplete/integrator/cf4.py @@ -106,136 +106,3 @@ def _cf4_f4(chain, rates): 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:: - \begin{aligned} - 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 - \end{aligned} - - Parameters - ---------- - operator : openmc.deplete.TransportOperator - The operator object to simulate on. - timesteps : iterable of float - Array of timesteps in units of [s]. Note that values are not cumulative. - power : float or iterable of float, optional - Power of the reactor in [W]. A single value indicates that the power is - constant over all timesteps. An iterable indicates potentially different - power levels for each timestep. For a 2D problem, the power can be given - in [W/cm] as long as the "volume" assigned to a depletion material is - actually an area in [cm^2]. Either `power` or `power_density` must be - specified. - power_density : float or iterable of float, optional - Power density of the reactor in [W/gHM]. It is multiplied by initial - heavy metal inventory to get total power if `power` is not speficied. - print_out : bool, optional - Whether or not to print out time. - """ - if power is None: - if power_density is None: - raise ValueError( - "Neither power nor power density was specified.") - if not isinstance(power_density, Iterable): - power = power_density*operator.heavy_metal - else: - power = [i*operator.heavy_metal for i in power_density] - - if not isinstance(power, Iterable): - power = [power]*len(timesteps) - - # Generate initial conditions - with operator as vec: - # Initialize time and starting index - if operator.prev_res is None: - t = 0.0 - i_res = 0 - else: - t = operator.prev_res[-1].time[-1] - i_res = len(operator.prev_res) - - 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 - # 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) - time_1, x_new = timed_deplete( - chain, x[0], op_results[0].rates, dt, print_out, - matrix_func=_cf4_f1) - x.append(x_new) - op_results.append(operator(x_new, p)) - - # Step 2: deplete with matrix 1/2*A(y1) - time_2, x_new = timed_deplete( - chain, x[0], op_results[1].rates, dt, print_out, - matrix_func=_cf4_f1) - x.append(x_new) - op_results.append(operator(x_new, p)) - - # Step 3: deplete with matrix -1/2*A(y0)+A(y2) - rates = list(zip(op_results[0].rates, op_results[2].rates)) - time_3, x_new = timed_deplete( - chain, x[1], rates, dt, print_out, matrix_func=_cf4_f2) - x.append(x_new) - op_results.append(operator(x_new, p)) - - # Step 4: deplete with two matrix exponentials - rates = list(zip(op_results[0].rates, op_results[1].rates, - op_results[2].rates, op_results[3].rates)) - time_4, x_end = timed_deplete( - chain, x[0], rates, dt, print_out, matrix_func=_cf4_f3) - time_5, x_end = timed_deplete( - chain, x_end, rates, dt, print_out, matrix_func=_cf4_f4) - - # Create results, write to disk - Results.save( - operator, x, op_results, [t, t + dt], p, i_res + i, - time_1 + time_2 + time_3 + time_4 + time_5) - - # Advance time, update vector - t += dt - vec = copy.deepcopy(x_end) - - # Perform one last simulation - x = [copy.deepcopy(vec)] - op_results = [operator(x[0], power[-1])] - - # Create results, write to disk - Results.save(operator, x, op_results, [t, t], p, i_res + len(timesteps))