Clean-up for PinTSimE + Update for SwitchEstimator (Preparation for PR for paper stuff)

pySDC/implementations/problem_classes/Battery.py

@@ -60,16 +60,24 @@ class battery_n_capacitors(ptype):
     dtype_u = mesh
     dtype_f = imex_mesh
 
-    def __init__(self, ncapacitors=2, Vs=5.0, Rs=0.5, C=None, R=1.0, L=1.0, alpha=1.2, V_ref=None):
+    def __init__(self, ncapacitors=1, Vs=5.0, Rs=0.5, C=None, R=1.0, L=1.0, alpha=1.2, V_ref=None):
         """Initialization routine"""
         n = ncapacitors
         nvars = n + 1
 
-        if C is None:
-            C = np.array([1.0, 1.0])
+        if C is None and ncapacitors == 1:
+            C = np.array([1.0])
+        else:
+            assert type(C) == np.ndarray, '"C" needs to be an np.ndarray'
+            assert np.shape(C)[0] == n, 'Number of capacitance values needs to be equal to number of condensators'
 
-        if V_ref is None:
-            V_ref = np.array([1.0, 1.0])
+        if V_ref is None and ncapacitors == 1:
+            V_ref = np.array([1.0])
+        else:
+            assert (alpha > V_ref[k] for k in range(n)), 'Please set "alpha" greater than values of "V_ref"'
+            assert type(V_ref) == np.ndarray, '"V_ref" needs to be an np.ndarray'
+            assert np.shape(V_ref)[0] == n, 'Number of reference values needs to be equal to number of condensators'
+            assert (V_ref[k] > 0 for k in range(n)), 'Please set values of "V_ref" greater than 0'
 
         # invoke super init, passing number of dofs, dtype_u and dtype_f
         super().__init__(init=(nvars, None, np.dtype('float64')))
@@ -230,6 +238,7 @@ def get_switching_info(self, u, t):
         switch_detected = False
         m_guess = -100
         break_flag = False
+        k_detected = 1
         for m in range(1, len(u)):
             for k in range(1, self.nvars):
                 h_prev_node = u[m - 1][k] - self.V_ref[k - 1]
@@ -244,9 +253,7 @@ def get_switching_info(self, u, t):
             if break_flag:
                 break
 
-        state_function = (
-            [u[m][k_detected] - self.V_ref[k_detected - 1] for m in range(len(u))] if switch_detected else []
-        )
+        state_function = [u[m][k_detected] - self.V_ref[k_detected - 1] for m in range(len(u))]
 
         return switch_detected, m_guess, state_function
 
@@ -443,11 +450,20 @@ def __init__(
         newton_maxiter=100,
         newton_tol=1e-11,
     ):
+        n = ncapacitors
         if C is None:
             C = np.array([1.0])
+        else:
+            assert type(C) == np.ndarray, '"C" needs to be an np.ndarray'
+            assert np.shape(C)[0] == n, 'Number of capacitance values needs to be equal to number of condensators'
 
         if V_ref is None:
             V_ref = np.array([1.0])
+        else:
+            assert alpha > V_ref[0], 'Please set "alpha" greater than values of "V_ref"'
+            assert type(V_ref) == np.ndarray, '"V_ref" needs to be an np.ndarray'
+            assert np.shape(V_ref)[0] == n, 'Number of reference values needs to be equal to number of condensators'
+            assert V_ref[0] > 0, 'Please "V_ref" greater than 0'
 
         super().__init__(ncapacitors, Vs, Rs, C, R, L, alpha, V_ref)
         self._makeAttributeAndRegister('newton_maxiter', 'newton_tol', localVars=locals(), readOnly=True)

pySDC/playgrounds/EnergyGrids/battery_adaptivity.py

@@ -2,71 +2,85 @@
 import dill
 
 from pySDC.helpers.stats_helper import get_sorted
-from pySDC.implementations.collocations import Collocation
-from pySDC.projects.PinTSimE.switch_controller_nonMPI import switch_controller_nonMPI
+from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
 from pySDC.implementations.problem_classes.Battery import battery
 from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
 
 from pySDC.playgrounds.EnergyGrids.log_data_battery import log_data_battery
-from pySDC.projects.PinTSimE.piline_model import setup_mpl
 import pySDC.helpers.plot_helper as plt_helper
 
+from pySDC.implementations.hooks.log_solution import LogSolution
+
+from pySDC.projects.PinTSimE.switch_estimator import SwitchEstimator
+from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
+from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI
+
 
 def main():
     """
     A simple test program to do SDC/PFASST runs for the battery drain model
     """
 
     # initialize level parameters
-    level_params = dict()
-    level_params['restol'] = -1e-10
-    level_params['e_tol'] = 1e-5
-    level_params['dt'] = 2e-2
+    level_params = {
+        'restol': -1,
+        'dt': 2e-2,
+    }
 
     # initialize sweeper parameters
-    sweeper_params = dict()
-
-    sweeper_params['quad_type'] = 'LOBATTO'
-    sweeper_params['num_nodes'] = 5
-    sweeper_params['QI'] = 'LU'  # For the IMEX sweeper, the LU-trick can be activated for the implicit part
-    sweeper_params['initial_guess'] = 'spread'
-
-    # initialize problem parameters
-    problem_params = dict()
-    problem_params['Vs'] = 5.0
-    problem_params['Rs'] = 0.5
-    problem_params['C'] = 1
-    problem_params['R'] = 1
-    problem_params['L'] = 1
-    problem_params['alpha'] = 3  # 10
-    problem_params['V_ref'] = 1
+    sweeper_params = {
+        'quad_type': 'LOBATTO',
+        'num_nodes': 5,
+        'QI': 'LU',
+        'initial_guess': 'spread',
+    }
 
     # initialize step parameters
-    step_params = dict()
-    step_params['maxiter'] = 5
+    step_params = {
+        'maxiter': 10,
+    }
 
     # initialize controller parameters
-    controller_params = dict()
-    controller_params['use_adaptivity'] = True
-    controller_params['use_switch_estimator'] = True
-    controller_params['logger_level'] = 20
-    controller_params['hook_class'] = log_data_battery
+    controller_params = {
+        'logger_level': 30,
+        'hook_class': [LogSolution],
+        'mssdc_jac': False,
+    }
+
+    # convergence controllers
+    convergence_controllers = {}
+    switch_estimator_params = {
+        'tol': 1e-10,
+        'alpha': 1.0,
+    }
+    convergence_controllers.update({SwitchEstimator: switch_estimator_params})
+    adaptivity_params = {
+        'e_tol': 1e-7,
+    }
+    convergence_controllers.update({Adaptivity: adaptivity_params})
+    restarting_params = {
+        'max_restarts': 50,
+        'crash_after_max_restarts': False,
+    }
+    convergence_controllers.update({BasicRestartingNonMPI: restarting_params})
 
     # fill description dictionary for easy step instantiation
-    description = dict()
-    description['problem_class'] = battery
-    description['problem_params'] = problem_params
-    description['sweeper_class'] = imex_1st_order
-    description['sweeper_params'] = sweeper_params
-    description['level_params'] = level_params
-    description['step_params'] = step_params
+    description = {
+        'problem_class': battery,
+        'problem_params': {},  # use default problem params
+        'sweeper_class': imex_1st_order,
+        'sweeper_params': sweeper_params,
+        'level_params': level_params,
+        'step_params': step_params,
+        'convergence_controllers': convergence_controllers,
+    }
 
     # set time parameters
     t0 = 0.0
-    Tend = 4
+    Tend = 0.5
 
     # instantiate controller
-    controller = switch_controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
+    controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
 
     # get initial values on finest level
     P = controller.MS[0].levels[0].prob
@@ -91,19 +105,17 @@ def plot_voltages(cwd='./'):
     stats = dill.load(f)
     f.close()
 
-    # convert filtered statistics to list of iterations count, sorted by process
-    cL = get_sorted(stats, type='current L', sortby='time', recomputed=False)
-    vC = get_sorted(stats, type='voltage C', sortby='time', recomputed=False)
-
-    times = [v[0] for v in cL]
+    u = get_sorted(stats, type='u', sortby='time', recomputed=False)
+    t = np.array([me[0] for me in u])
+    iL = np.array([me[1][0] for me in u])
+    vC = np.array([me[1][1] for me in u])
 
     dt = np.array(get_sorted(stats, type='dt', recomputed=False))
     list_gs = get_sorted(stats, type='restart')
 
-    setup_mpl()
     fig, ax = plt_helper.plt.subplots(1, 1)
-    ax.plot(times, [v[1] for v in cL], label='$i_L$')
-    ax.plot(times, [v[1] for v in vC], label='$v_C$')
+    ax.plot(t, iL, label='$i_L$')
+    ax.plot(t, vC, label='$v_C$')
     ax.set_xlabel('Time', fontsize=20)
     ax.set_ylabel('Energy', fontsize=20)
     for element in list_gs:

pySDC/playgrounds/EnergyGrids/playground.py

@@ -2,12 +2,11 @@
 import dill
 
 from pySDC.helpers.stats_helper import get_sorted
-from pySDC.implementations.collocations import Collocation
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
 from pySDC.implementations.problem_classes.Piline import piline
 from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
 
-from pySDC.playgrounds.EnergyGrids.log_data import log_data
+from pySDC.implementations.hooks.log_solution import LogSolution
 import pySDC.helpers.plot_helper as plt_helper
 
 
@@ -45,7 +44,7 @@ def main():
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 20
-    controller_params['hook_class'] = log_data
+    controller_params['hook_class'] = [LogSolution]
 
     # fill description dictionary for easy step instantiation
     description = dict()
@@ -82,16 +81,16 @@ def plot_voltages(cwd='./'):
     f.close()
 
     # convert filtered statistics to list of iterations count, sorted by process
-    v1 = get_sorted(stats, type='v1', sortby='time')
-    v2 = get_sorted(stats, type='v2', sortby='time')
-    p3 = get_sorted(stats, type='p3', sortby='time')
-
-    times = [v[0] for v in v1]
+    u_val = get_sorted(stats, type='u', sortby='time')
+    t = np.array([me[0] for me in u_val])
+    v1 = np.array([me[1][0] for me in u_val])
+    v2 = np.array([me[1][1] for me in u_val])
+    p3 = np.array([me[1][2] for me in u_val])
 
     # plt_helper.setup_mpl()
-    plt_helper.plt.plot(times, [v[1] for v in v1], label='v1')
-    plt_helper.plt.plot(times, [v[1] for v in v2], label='v2')
-    plt_helper.plt.plot(times, [v[1] for v in p3], label='p3')
+    plt_helper.plt.plot(t, v1, label='v1')
+    plt_helper.plt.plot(t, v2, label='v2')
+    plt_helper.plt.plot(t, p3, label='p3')
     plt_helper.plt.legend()
 
     plt_helper.plt.show()