moved more algorithms from Colin: epc-rk4, si-celi, si-leqi

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

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)

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)