Brand-new switch estimator + new discontinuous ODE model! (#338)

* Added more reasonable variable names and changes in descriptions

* Discontinuous test problem + variable and description changes for SE

* Update for switch estimator

* Start with playing around with DisTestDAE

* Big update for the switch estimator!

* Discontinuous ODE model + tests

* read only Parameters changed

* Black

* Typo

* Events on boundary or with high resolution should also count

* Fix for battery problem with residual tolerance

* Changes for PR

* Lint

* Added warning for computing local/global error

* Warning only for discontinuous test ODE problem

* Requested changes from @brownbaerchen

* Removed warning for piline

* Change to LogGlobalErrorPostStep hook

---------

Co-authored-by: Thomas Baumann <t.baumann@fz-juelich.de>

Author: lisawim

pySDC/implementations/problem_classes/Auzinger_implicit.py

@@ -6,7 +6,7 @@
 
 # noinspection PyUnusedLocal
 class auzinger(ptype):
-    """
+    r"""
     This class implements the Auzinger equation as initial value problem. It can be found in doi.org/10.2140/camcos.2015.10.1.
     The system of two ordinary differential equations (ODEs) is given by
 

pySDC/implementations/problem_classes/Battery.py

@@ -77,8 +77,9 @@ def __init__(self, ncapacitors=2, Vs=5.0, Rs=0.5, C=None, R=1.0, L=1.0, alpha=1.
             'nvars', 'ncapacitors', 'Vs', 'Rs', 'C', 'R', 'L', 'alpha', 'V_ref', localVars=locals(), readOnly=True
         )
 
-        self.A = np.zeros((n + 1, n + 1))
         self.switch_A, self.switch_f = self.get_problem_dict()
+        self.A = self.switch_A[0]
+
         self.t_switch = None
         self.nswitches = 0
 
@@ -222,17 +223,18 @@ def get_switching_info(self, u, t):
             Indicates if a switch is found or not.
         m_guess : 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).
+        state_function : list
+            Contains values of the state function (for interpolation).
         """
 
         switch_detected = False
         m_guess = -100
         break_flag = False
-
         for m in range(1, len(u)):
             for k in range(1, self.nvars):
-                if u[m][k] - self.V_ref[k - 1] <= 0:
+                h_prev_node = u[m - 1][k] - self.V_ref[k - 1]
+                h_curr_node = u[m][k] - self.V_ref[k - 1]
+                if h_prev_node > 0 and h_curr_node <= 0:
                     switch_detected = True
                     m_guess = m - 1
                     k_detected = k
@@ -242,9 +244,11 @@ def get_switching_info(self, u, t):
             if break_flag:
                 break
 
-        vC_switch = [u[m][k_detected] - self.V_ref[k_detected - 1] for m in range(1, len(u))] if switch_detected else []
+        state_function = (
+            [u[m][k_detected] - self.V_ref[k_detected - 1] for m in range(len(u))] if switch_detected else []
+        )
 
-        return switch_detected, m_guess, vC_switch
+        return switch_detected, m_guess, state_function
 
     def count_switches(self):
         """
@@ -262,7 +266,6 @@ def get_problem_dict(self):
         n = self.ncapacitors
         v = np.zeros(n + 1)
         v[0] = 1
-
         A, f = dict(), dict()
         A = {k: np.diag(-1 / (self.C[k] * self.R) * np.roll(v, k + 1)) for k in range(n)}
         A.update({n: np.diag(-(self.Rs + self.R) / self.L * v)})
@@ -273,23 +276,28 @@ def get_problem_dict(self):
 
 class battery(battery_n_capacitors):
     r"""
-    Example implementing the battery drain model with :math:`N=1` capacitor, inherits from battery_n_capacitors. The ODE system
-    of this model is given by the following equations. If :math:`v_1 > V_{ref, 0}:`
+    Example implementing the battery drain model with :math:`N=1` capacitor, inherits from battery_n_capacitors. This model is an example
+    of a discontinuous problem. The state function :math:`decides` which differential equation is solved. When the state function has
+    a sign change the dynamics of the solution changes by changing the differential equation. The ODE system of this model is given by
+    the following equations:
+
+    If :math:`h(v_1) := v_1 - V_{ref, 0} > 0:`
 
     .. math::
         \frac{d i_L (t)}{dt} = 0,
 
     .. math::
         \frac{d v_1 (t)}{dt} = -\frac{1}{CR}v_1 (t),
 
-    where :math:`i_L` denotes the function of the current over time :math:`t`.
-    If :math:`v_1 \leq V_{ref, 0}:`
+    else:
 
     .. math::
         \frac{d i_L(t)}{dt} = -\frac{R_s + R}{L}i_L (t) + \frac{1}{L} V_s,
 
     .. math::
-        \frac{d v_1(t)}{dt} = 0.
+        \frac{d v_1(t)}{dt} = 0,
+
+    where :math:`i_L` denotes the function of the current over time :math:`t`.
 
     Note
     ----
@@ -320,7 +328,7 @@ def eval_f(self, u, t):
 
         t_switch = np.inf if self.t_switch is None else self.t_switch
 
-        if u[1] <= self.V_ref[0] or t >= t_switch:
+        if u[1] - self.V_ref[0] <= 0 or t >= t_switch:
             f.expl[0] = self.Vs / self.L
 
         else:
@@ -352,7 +360,7 @@ def solve_system(self, rhs, factor, u0, t):
 
         t_switch = np.inf if self.t_switch is None else self.t_switch
 
-        if rhs[1] <= self.V_ref[0] or t >= t_switch:
+        if rhs[1] - self.V_ref[0] <= 0 or t >= t_switch:
             self.A[0, 0] = -(self.Rs + self.R) / self.L
 
         else:
@@ -432,8 +440,8 @@ def __init__(
         L=1.0,
         alpha=1.2,
         V_ref=None,
-        newton_maxiter=200,
-        newton_tol=1e-8,
+        newton_maxiter=100,
+        newton_tol=1e-11,
     ):
         if C is None:
             C = np.array([1.0])
@@ -469,7 +477,7 @@ def eval_f(self, u, t):
 
         t_switch = np.inf if self.t_switch is None else self.t_switch
 
-        if u[1] <= self.V_ref[0] or t >= t_switch:
+        if u[1] - self.V_ref[0] <= 0 or t >= t_switch:
             self.A[0, 0] = -(self.Rs + self.R) / self.L
             non_f[0] = self.Vs
 
@@ -507,7 +515,7 @@ def solve_system(self, rhs, factor, u0, t):
 
         t_switch = np.inf if self.t_switch is None else self.t_switch
 
-        if rhs[1] <= self.V_ref[0] or t >= t_switch:
+        if rhs[1] - self.V_ref[0] <= 0 or t >= t_switch:
             self.A[0, 0] = -(self.Rs + self.R) / self.L
             non_f[0] = self.Vs
 

pySDC/implementations/problem_classes/DiscontinuousTestODE.py

@@ -0,0 +1,226 @@
+import numpy as np
+
+from pySDC.core.Errors import ParameterError, ProblemError
+from pySDC.core.Problem import ptype
+from pySDC.implementations.datatype_classes.mesh import mesh
+
+
+class DiscontinuousTestODE(ptype):
+    r"""
+    This class implements a very simple test case of a ordinary differential equation consisting of one discrete event. The dynamics of
+    the solution changes when the state function :math:`h(u) := u - 5` changes the sign. The problem is defined by:
+
+    if :math:`u - 5 < 0:`
+
+        .. math::
+            \fra{d u}{dt} = u
+
+    else:
+
+        .. math::
+            \frac{d u}{dt} = \frac{4}{t^*},
+
+    where :math:`t^* = \log(5) \approx 1.6094379`. For :math:`h(u) < 0`, i.e., :math:`t \leq t^*` the exact solution is
+    :math:`u(t) = exp(t)`; for :math:`h(u) \geq 0`, i.e., :math:`t \geq t^*` the exact solution is :math:`u(t) = \frac{4 t}{t^*} + 1`.
+
+    Attributes
+    ----------
+    t_switch_exact : float
+        Exact event time with :math:`t^* = \log(5)`.
+    t_switch: float
+        Time point of the discrete event found by switch estimation.
+    nswitches: int
+        Number of switches found by switch estimation.
+    newton_itercount: int
+        Counts the number of Newton iterations.
+    newton_ncalls: int
+        Counts the number of how often Newton is called in the simulation of the problem.
+    """
+
+    dtype_u = mesh
+    dtype_f = mesh
+
+    def __init__(self, newton_maxiter=100, newton_tol=1e-8):
+        """Initialization routine"""
+        nvars = 1
+        super().__init__(init=(nvars, None, np.dtype('float64')))
+        self._makeAttributeAndRegister('nvars', localVars=locals(), readOnly=True)
+        self._makeAttributeAndRegister('newton_maxiter', 'newton_tol', localVars=locals())
+
+        if self.nvars != 1:
+            raise ParameterError('nvars has to be equal to 1!')
+
+        self.t_switch_exact = np.log(5)
+        self.t_switch = None
+        self.nswitches = 0
+        self.newton_itercount = 0
+        self.newton_ncalls = 0
+
+    def eval_f(self, u, t):
+        """
+        Routine to evaluate the right-hand side of the problem.
+
+        Parameters
+        ----------
+        u : dtype_u
+            Current values of the numerical solution.
+        t : float
+            Current time of the numerical solution is computed.
+
+        Returns
+        -------
+        f : dtype_f
+            The right-hand side of the problem.
+        """
+
+        t_switch = np.inf if self.t_switch is None else self.t_switch
+
+        f = self.dtype_f(self.init, val=0.0)
+        h = u[0] - 5
+        if h >= 0 or t >= t_switch:
+            f[:] = 4 / self.t_switch_exact
+        else:
+            f[:] = u
+        return f
+
+    def solve_system(self, rhs, dt, u0, t):
+        r"""
+        Simple Newton solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`.
+
+        Parameters
+        ----------
+        rhs : dtype_f
+            Right-hand side for the linear system.
+        dt : 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
+        -------
+        me : dtype_u
+            The solution as mesh.
+        """
+
+        t_switch = np.inf if self.t_switch is None else self.t_switch
+
+        h = rhs[0] - 5
+        u = self.dtype_u(u0)
+
+        n = 0
+        res = 99
+        while n < self.newton_maxiter:
+            # form function g with g(u) = 0
+            if h >= 0 or t >= t_switch:
+                g = u - dt * (4 / self.t_switch_exact) - rhs
+            else:
+                g = u - dt * u - rhs
+
+            # if g is close to 0, then we are done
+            res = np.linalg.norm(g, np.inf)
+
+            if res < self.newton_tol:
+                break
+
+            if h >= 0 or t >= t_switch:
+                dg = 1
+            else:
+                dg = 1 - dt
+
+            # newton update
+            u -= 1.0 / dg * g
+
+            n += 1
+
+        if np.isnan(res) and self.stop_at_nan:
+            raise ProblemError('Newton got nan after %i iterations, aborting...' % n)
+        elif np.isnan(res):
+            self.logger.warning('Newton got nan after %i iterations...' % n)
+
+        if n == self.newton_maxiter:
+            self.logger.warning('Newton did not converge after %i iterations, error is %s' % (n, res))
+
+        self.newton_ncalls += 1
+        self.newton_itercount += n
+
+        me = self.dtype_u(self.init)
+        me[:] = u[:]
+
+        return me
+
+    def u_exact(self, t, u_init=None, t_init=None):
+        """
+        Routine to compute the exact solution at time t.
+
+        Parameters
+        ----------
+        t : float
+            Time of the exact solution.
+        u_init : pySDC.problem.DiscontinuousTestODE.dtype_u
+            Initial conditions for getting the exact solution.
+        t_init : float
+            The starting time.
+
+        Returns
+        -------
+        me : dtype_u
+            The exact solution.
+        """
+
+        if t_init is not None and u_init is not None:
+            self.logger.warning(
+                f'{type(self).__name__} uses an analytic exact solution from t=0. If you try to compute the local error, you will get the global error instead!'
+            )
+
+        me = self.dtype_u(self.init)
+        if t <= self.t_switch_exact:
+            me[:] = np.exp(t)
+        else:
+            me[:] = (4 * t) / self.t_switch_exact + 1
+        return me
+
+    def get_switching_info(self, u, t):
+        """
+        Provides information about the state function of the problem. When the state function changes its sign,
+        typically an event occurs. So the check for an event should be done in the way that the state function
+        is checked for a sign change. If this is the case, the intermediate value theorem states a root in this
+        step.
+
+        Parameters
+        ----------
+        u : dtype_u
+            Current values of the numerical solution at time t.
+        t : float
+            Current time of the numerical solution.
+
+        Returns
+        -------
+        switch_detected : bool
+            Indicates whether a discrete event is found or not.
+        m_guess : int
+            The index before the sign changes.
+        state_function : list
+            Defines the values of the state function at collocation nodes where it changes the sign.
+        """
+
+        switch_detected = False
+        m_guess = -100
+
+        for m in range(1, len(u)):
+            h_prev_node = u[m - 1][0] - 5
+            h_curr_node = u[m][0] - 5
+            if h_prev_node < 0 and h_curr_node >= 0:
+                switch_detected = True
+                m_guess = m - 1
+                break
+
+        state_function = [u[m][0] - 5 for m in range(len(u))] if switch_detected else []
+        return switch_detected, m_guess, state_function
+
+    def count_switches(self):
+        """
+        Setter to update the number of switches if one is found.
+        """
+        self.nswitches += 1

pySDC/implementations/problem_classes/Lorenz.py

@@ -192,6 +192,7 @@ def u_exact(self, t, u_init=None, t_init=None):
         me : dtype_u
             The approximated exact solution.
         """
+
         me = self.dtype_u(self.init)
 
         if t > 0:

pySDC/implementations/problem_classes/TestEquation_0D.py

@@ -125,6 +125,7 @@ def u_exact(self, t, u_init=None, t_init=None):
         me : dtype_u
             The exact solution.
         """
+
         u_init = (self.u0 if u_init is None else u_init) * 1.0
         t_init = 0.0 if t_init is None else t_init * 1.0