Status variables, short switching logics and other improvements, thanks to @pancetta and @brownbaerchen!

Author: lisawim

pySDC/implementations/problem_classes/Battery.py

@@ -8,6 +8,7 @@
 class battery(ptype):
     """
     Example implementing the battery drain model as in the description in the PinTSimE project
+
     Attributes:
         A: system matrix, representing the 2 ODEs
         t_switch: time point of the switch
@@ -17,21 +18,22 @@ class battery(ptype):
     def __init__(self, problem_params, dtype_u=mesh, dtype_f=imex_mesh):
         """
         Initialization routine
+
         Args:
             problem_params (dict): custom parameters for the example
             dtype_u: mesh data type for solution
             dtype_f: mesh data type for RHS
         """
 
-        problem_params['nvars'] = 2
-
         # these parameters will be used later, so assert their existence
-        essential_keys = ['Vs', 'Rs', 'C', 'R', 'L', 'alpha', 'V_ref']
+        essential_keys = ['ncondensators', 'Vs', 'Rs', 'C', 'R', 'L', 'alpha', 'V_ref']
         for key in essential_keys:
             if key not in problem_params:
                 msg = 'need %s to instantiate problem, only got %s' % (key, str(problem_params.keys()))
                 raise ParameterError(msg)
 
+        problem_params['nvars'] = problem_params['ncondensators'] + 1
+
         # invoke super init, passing number of dofs, dtype_u and dtype_f
         super(battery, self).__init__(
             init=(problem_params['nvars'], None, np.dtype('float64')),
@@ -48,56 +50,50 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=imex_mesh):
     def eval_f(self, u, t):
         """
         Routine to evaluate the RHS
+
         Args:
             u (dtype_u): current values
             t (float): current time
+
         Returns:
             dtype_f: the RHS
         """
 
         f = self.dtype_f(self.init, val=0.0)
         f.impl[:] = self.A.dot(u)
 
-        if self.t_switch is not None:
-            if t >= self.t_switch:
-                f.expl[0] = self.params.Vs / self.params.L
+        t_switch = np.inf if self.t_switch is None else self.t_switch
 
-            else:
-                f.expl[0] = 0
+        if u[1] <= self.params.V_ref or t >= t_switch:
+            f.expl[0] = self.params.Vs / self.params.L
 
         else:
-            if u[1] <= self.params.V_ref:
-                f.expl[0] = self.params.Vs / self.params.L
-
-            else:
-                f.expl[0] = 0
+            f.expl[0] = 0
 
         return f
 
     def solve_system(self, rhs, factor, u0, t):
         """
         Simple linear solver for (I-factor*A)u = rhs
+
         Args:
             rhs (dtype_f): right-hand side for the linear system
             factor (float): abbrev. for the local stepsize (or any other factor required)
             u0 (dtype_u): initial guess for the iterative solver
             t (float): current time (e.g. for time-dependent BCs)
+
         Returns:
             dtype_u: solution as mesh
         """
         self.A = np.zeros((2, 2))
 
-        if self.t_switch is not None:
-            if t >= self.t_switch:
-                self.A[0, 0] = -(self.params.Rs + self.params.R) / self.params.L
-            else:
-                self.A[1, 1] = -1 / (self.params.C * self.params.R)
+        t_switch = np.inf if self.t_switch is None else self.t_switch
+
+        if rhs[1] <= self.params.V_ref or t >= t_switch:
+            self.A[0, 0] = -(self.params.Rs + self.params.R) / self.params.L
 
         else:
-            if rhs[1] <= self.params.V_ref:
-                self.A[0, 0] = -(self.params.Rs + self.params.R) / self.params.L
-            else:
-                self.A[1, 1] = -1 / (self.params.C * self.params.R)
+            self.A[1, 1] = -1 / (self.params.C * self.params.R)
 
         me = self.dtype_u(self.init)
         me[:] = np.linalg.solve(np.eye(self.params.nvars) - factor * self.A, rhs)
@@ -106,8 +102,10 @@ def solve_system(self, rhs, factor, u0, t):
     def u_exact(self, t):
         """
         Routine to compute the exact solution at time t
+
         Args:
             t (float): current time
+
         Returns:
             dtype_u: exact solution
         """
@@ -123,12 +121,14 @@ def u_exact(self, t):
     def get_switching_info(self, u, t):
         """
         Provides information about a discrete event for one subinterval.
+
         Args:
             u (dtype_u): current values
             t (float): current time
+
         Returns:
             switch_detected (bool): Indicates if a switch is found or not
-            m_guess (np.int): Index of where the discrete event would found
+            m_guess (np.int): Index of collocation node inside one subinterval of where the discrete event was found
             vC_switch (list): Contains function values of switching condition (for interpolation)
         """
 
@@ -141,10 +141,7 @@ def get_switching_info(self, u, t):
                 m_guess = m - 1
                 break
 
-        vC_switch = []
-        if switch_detected:
-            for m in range(1, len(u)):
-                vC_switch.append(u[m][1] - self.params.V_ref)
+        vC_switch = [u[m][1] - self.params.V_ref for m in range(1, len(u))] if switch_detected else []
 
         return switch_detected, m_guess, vC_switch
 
@@ -158,21 +155,21 @@ def flip_switches(self):
 
 
 class battery_implicit(battery):
-
     def __init__(self, problem_params, dtype_u=mesh, dtype_f=mesh):
 
-        essential_keys = ['newton_maxiter', 'newton_tol']
+        essential_keys = ['newton_maxiter', 'newton_tol', 'ncondensators', 'Vs', 'Rs', 'C', 'R', 'L', 'alpha', 'V_ref']
         for key in essential_keys:
             if key not in problem_params:
                 msg = 'need %s to instantiate problem, only got %s' % (key, str(problem_params.keys()))
                 raise ParameterError(msg)
 
+        problem_params['nvars'] = problem_params['ncondensators'] + 1
+
         # invoke super init, passing number of dofs, dtype_u and dtype_f
         super(battery_implicit, self).__init__(
-            init=(problem_params['nvars'], None, np.dtype('float64')),
+            problem_params,
             dtype_u=dtype_u,
             dtype_f=dtype_f,
-            params=problem_params,
         )
 
         self.newton_itercount = 0
@@ -183,43 +180,41 @@ def __init__(self, problem_params, dtype_u=mesh, dtype_f=mesh):
     def eval_f(self, u, t):
         """
         Routine to evaluate the RHS
+
         Args:
             u (dtype_u): current values
             t (float): current time
+
         Returns:
             dtype_f: the RHS
         """
 
         f = self.dtype_f(self.init, val=0.0)
         non_f = np.zeros(2)
 
-        if self.t_switch is not None:
-            if t >= self.t_switch:
-                self.A[0, 0] = -(self.params.Rs + self.params.R) / self.params.L
-                non_f[0] = self.params.Vs / self.params.L
-            else:
-                self.A[1, 1] = -1 / (self.params.C * self.params.R)
-                non_f[0] = 0
+        t_switch = np.inf if self.t_switch is None else self.t_switch
+
+        if u[1] <= self.params.V_ref or t >= t_switch:
+            self.A[0, 0] = -(self.params.Rs + self.params.R) / self.params.L
+            non_f[0] = self.params.Vs
 
         else:
-            if u[1] <= self.params.V_ref:
-                self.A[0, 0] = -(self.params.Rs + self.params.R) / self.params.L
-                non_f[0] = self.params.Vs / self.params.L
-            else:
-                self.A[1, 1] = -1 / (self.params.C * self.params.R)
-                non_f[0] = 0
+            self.A[1, 1] = -1 / (self.params.C * self.params.R)
+            non_f[0] = 0
 
         f[:] = self.A.dot(u) + non_f
         return f
 
     def solve_system(self, rhs, factor, u0, t):
         """
         Simple Newton solver
+
         Args:
             rhs (dtype_f): right-hand side for the linear system
             factor (float): abbrev. for the local stepsize (or any other factor required)
             u0 (dtype_u): initial guess for the iterative solver
             t (float): current time (e.g. for time-dependent BCs)
+
         Returns:
             dtype_u: solution as mesh
         """
@@ -228,21 +223,15 @@ def solve_system(self, rhs, factor, u0, t):
         non_f = np.zeros(2)
         self.A = np.zeros((2, 2))
 
-        if self.t_switch is not None:
-            if t >= self.t_switch:
-                self.A[0, 0] = -(self.params.Rs + self.params.R) / self.params.L
-                non_f[0] = self.params.Vs / self.params.L
-            else:
-                self.A[1, 1] = -1 / (self.params.C * self.params.R)
-                non_f[0] = 0
+        t_switch = np.inf if self.t_switch is None else self.t_switch
+
+        if rhs[1] <= self.params.V_ref or t >= t_switch:
+            self.A[0, 0] = -(self.params.Rs + self.params.R) / self.params.L
+            non_f[0] = self.params.Vs
 
         else:
-            if rhs[1] <= self.params.V_ref:
-                self.A[0, 0] = -(self.params.Rs + self.params.R) / self.params.L
-                non_f[0] = self.params.Vs / self.params.L
-            else:
-                self.A[1, 1] = -1 / (self.params.C * self.params.R)
-                non_f[0] = 0
+            self.A[1, 1] = -1 / (self.params.C * self.params.R)
+            non_f[0] = 0
 
         # start newton iteration
         n = 0

pySDC/implementations/problem_classes/Battery_2Condensators.py

@@ -164,10 +164,9 @@ def get_switching_info(self, u, t):
                     msg = 'A discrete event is already found! Multiple switching handling in the same interval is not yet implemented!'
                     raise AssertionError(msg)
 
-        vC_switch = []
-        if switch_detected:
-            for m in range(1, len(u)):
-                vC_switch.append(u[m][k_detected] - self.params.V_ref[k_detected - 1])
+        vC_switch = (
+            [u[m][k_detected] - self.params.V_ref[k_detected - 1] for m in range(1, len(u))] if switch_detected else []
+        )
 
         return switch_detected, m_guess, vC_switch
 

pySDC/projects/PinTSimE/battery_2condensators_model.py

@@ -67,6 +67,12 @@ def post_step(self, step, level_number):
 def main(use_switch_estimator=True):
     """
     A simple test program to do SDC/PFASST runs for the battery drain model using 2 condensators
+
+    Args:
+        use_switch_estimator (bool): flag if the switch estimator wants to be used or not
+
+    Returns:
+        description (dict): contains all information for a controller run
     """
 
     # initialize level parameters
@@ -83,6 +89,7 @@ def main(use_switch_estimator=True):
 
     # initialize problem parameters
     problem_params = dict()
+    problem_params['ncondensators'] = 2
     problem_params['Vs'] = 5.0
     problem_params['Rs'] = 0.5
     problem_params['C1'] = 1.0
@@ -98,7 +105,7 @@ def main(use_switch_estimator=True):
 
     # initialize controller parameters
     controller_params = dict()
-    controller_params['logger_level'] = 15
+    controller_params['logger_level'] = 30
     controller_params['hook_class'] = log_data
 
     # convergence controllers
@@ -171,12 +178,18 @@ def main(use_switch_estimator=True):
 
     plot_voltages(description, recomputed, use_switch_estimator)
 
-    return np.mean(niters)
+    return description
 
 
 def plot_voltages(description, recomputed, use_switch_estimator, cwd='./'):
     """
     Routine to plot the numerical solution of the model
+
+    Args:
+        description(dict): contains all information for a controller run
+        recomputed (bool): flag if the values after a restart are used or before
+        use_switch_estimator (bool): flag if the switch estimator wants to be used or not
+        cwd: current working directory
     """
 
     f = open(cwd + 'data/battery_2condensators.dat', 'rb')
@@ -200,7 +213,7 @@ def plot_voltages(description, recomputed, use_switch_estimator, cwd='./'):
         switches = get_recomputed(stats, type='switch', sortby='time')
 
         if recomputed is not None:
-            assert len(switches) >= 1 and len(switches) >= 2, "No switches found"
+            assert len(switches) >= 2, f"Expected at least 2 switches, got {len(switches)}!"
         t_switches = [v[1] for v in switches]
 
         for i in range(len(t_switches)):
@@ -218,6 +231,10 @@ def plot_voltages(description, recomputed, use_switch_estimator, cwd='./'):
 def proof_assertions_description(description, use_switch_estimator):
     """
     Function to proof the assertions (function to get cleaner code)
+
+    Args:
+        description(dict): contains all information for a controller run
+        use_switch_estimator (bool): flag if the switch estimator wants to be used or not
     """
 
     assert (
@@ -227,7 +244,13 @@ def proof_assertions_description(description, use_switch_estimator):
         description['problem_params']['alpha'] > description['problem_params']['V_ref'][1]
     ), 'Please set "alpha" greater than "V_ref2"'
 
-    assert type(description['problem_params']['V_ref']) == np.ndarray, '"V_ref" needs to be an array (of type float)'
+    if description['problem_params']['ncondensators'] > 1:
+        assert (
+            type(description['problem_params']['V_ref']) == np.ndarray
+        ), '"V_ref" needs to be an array (of type float)'
+        assert (
+            description['problem_params']['ncondensators'] == np.shape(description['problem_params']['V_ref'])[0]
+        ), 'Number of reference values needs to be equal to number of condensators'
 
     assert description['problem_params']['V_ref'][0] > 0, 'Please set "V_ref1" greater than 0'
     assert description['problem_params']['V_ref'][1] > 0, 'Please set "V_ref2" greater than 0'