Datatype DAEMesh for DAEs

pySDC/projects/DAE/misc/HookClass_DAE.py

@@ -3,24 +3,27 @@
 
 class approx_solution_hook(hooks):
     """
-    Hook class to add the approximate solution to the output generated by the sweeper after each time step
+    Hook class to add the approximate solution to the output generated by the sweeper
+    after each time step.
     """
 
     def __init__(self):
-        """
-        Initialization routine for the custom hook
-        """
-        super(approx_solution_hook, self).__init__()
+        """Initialization routine"""
+        super().__init__()
 
     def post_step(self, step, level_number):
-        """
-        Default routine called after each step
-        Args:
-            step: the current step
-            level_number: the current level number
+        r"""
+        Default routine called after each step.
+
+        Parameters
+        ----------
+        step : pySDC.core.Step
+            Current step.
+        level_number : pySDC.core.level
+            Current level number.
         """
 
-        super(approx_solution_hook, self).post_step(step, level_number)
+        super().post_step(step, level_number)
 
         # some abbreviations
         L = step.levels[level_number]
@@ -41,24 +44,27 @@
 
 class error_hook(hooks):
     """
-    Hook class to add the approximate solution to the output generated by the sweeper after each time step
+    Hook class to log the error to the output generated by the sweeper after
+    each time step.
     """
 
     def __init__(self):
-        """
-        Initialization routine for the custom hook
-        """
-        super(error_hook, self).__init__()
+        """Initialization routine"""
+        super().__init__()
 
     def post_step(self, step, level_number):
-        """
-        Default routine called after each step
-        Args:
-            step: the current step
-            level_number: the current level number
+        r"""
+        Default routine called after each step.
+
+        Parameters
+        ----------
+        step : pySDC.core.Step
+            Current step.
+        level_number : pySDC.core.level
+            Current level number.
         """
 
-        super(error_hook, self).post_step(step, level_number)
+        super().post_step(step, level_number)
 
         # some abbreviations
         L = step.levels[level_number]
@@ -70,8 +76,7 @@
         # compute and save errors
         # Note that the component from which the error is measured is specified here
         upde = P.u_exact(step.time + step.dt)
-        err = abs(upde[0] - L.uend[0])
-        # err = abs(upde[4] - L.uend[4])
+        err = abs(upde.diff - L.uend.diff)
 
         self.add_to_stats(
             process=step.status.slot,
@@ -82,3 +87,41 @@
             type='error_post_step',
             value=err,
         )
+
+
+class LogGlobalErrorPostStepAlgebraicVariable(hooks):
+    """
+    Logs the global error in the algebraic variable
+    """
+
+    def post_step(self, step, level_number):
+        r"""
+        Default routine called after each step.
+
+        Parameters
+        ----------
+        step : pySDC.core.Step
+            Current step.
+        level_number : pySDC.core.level
+            Current level number.
+        """
+
+        super().post_step(step, level_number)
+
+        L = step.levels[level_number]
+        P = L.prob
+
+        L.sweep.compute_end_point()
+
+        upde = P.u_exact(step.time + step.dt)
+        e_global_algebraic = abs(upde.alg - L.uend.alg)
+
+        self.add_to_stats(
+            process=step.status.slot,
+            time=L.time + L.dt,
+            level=L.level_index,
+            iter=step.status.iter,
+            sweep=L.status.sweep,
+            type='e_global_algebraic_post_step',
+            value=e_global_algebraic,
+        )

pySDC/projects/DAE/misc/ProblemDAE.py

@@ -2,7 +2,7 @@
 from scipy.optimize import root
 
 from pySDC.core.Problem import ptype, WorkCounter
-from pySDC.implementations.datatype_classes.mesh import mesh
+from pySDC.projects.DAE.misc.dae_mesh import DAEMesh
 
 
 class ptype_dae(ptype):
@@ -25,8 +25,8 @@ class ptype_dae(ptype):
         in work_counters['rhs']
     """
 
-    dtype_u = mesh
-    dtype_f = mesh
+    dtype_u = DAEMesh
+    dtype_f = DAEMesh
 
     def __init__(self, nvars, newton_tol):
         """Initialization routine"""
@@ -54,14 +54,21 @@ def solve_system(self, impl_sys, u0, t):
         me : dtype_u
             Numerical solution.
         """
-
         me = self.dtype_u(self.init)
+
+        def implSysAsNumpy(unknowns, **kwargs):
+            me.diff[:] = unknowns[: np.size(me.diff)].reshape(me.diff.shape)
+            me.alg[:] = unknowns[np.size(me.diff) :].reshape(me.alg.shape)
+            sys = impl_sys(me, **kwargs)
+            return np.append(sys.diff.flatten(), sys.alg.flatten())  # TODO: more efficient way?
+
         opt = root(
-            impl_sys,
-            u0,
+            implSysAsNumpy,
+            np.append(u0.diff.flatten(), u0.alg.flatten()),
             method='hybr',
             tol=self.newton_tol,
         )
-        me[:] = opt.x
+        me.diff[:] = opt.x[: np.size(me.diff)].reshape(me.diff.shape)
+        me.alg[:] = opt.x[np.size(me.diff) :].reshape(me.alg.shape)
         self.work_counters['newton'].niter += opt.nfev
         return me

pySDC/projects/DAE/misc/dae_mesh.py

@@ -0,0 +1,99 @@
+import numpy as np
+
+from pySDC.implementations.datatype_classes.mesh import mesh
+from pySDC.core.Errors import DataError
+
+
+class DAEMesh(object):
+    r"""
+    Datatype for DAE problems. The solution of the problem can be splitted in the differential part
+    and in an algebraic part.
+
+    This data type can be used for the solution of the problem itself as well as for its derivative.
+
+    Parameters
+    ----------
+    init :
+        Initialization for a mesh. It can either be a tuple (one int per dimension) or a number
+        (if only one dimension is requested) or another imex_mesh object.
+    val : float, optional
+        The value that the mesh is wanted to fill with. Default is ``0.0``.
+
+    Attributes
+    ----------
+    diff : mesh.mesh
+        Differential part.
+    alg : mesh.mesh
+        Algebraic part.
+    """
+
+    def __init__(self, init, val=0.0):
+        """Initialization routine"""
+
+        if isinstance(init, type(self)):
+            self.diff = mesh(init.diff)
+            self.alg = mesh(init.alg)
+        elif (
+            isinstance(init, tuple)
+            and (init[1] is None or str(type(init[1])) == "MPI.Intracomm")
+            and isinstance(init[2], np.dtype)
+        ):
+            self.diff = mesh((init[0][0], init[1], init[2]), val=val)
+            self.alg = mesh((init[0][1], init[1], init[2]), val=val)
+
+        # something is wrong, if none of the ones above hit
+        else:
+            raise DataError('something went wrong during %s initialization' % type(self))
+
+    def __abs__(self):
+        """
+        It chooses the maximum value between the differential part and the algebraic part. If the
+        problem contains no algebraic part, then the maximum values is computed over the differential parts.
+        """
+        return max([abs(self.diff), abs(self.alg)]) if len(self.alg) > 0 else max([abs(self.diff)])
+
+    def __add__(self, other):
+        """
+        Overloading the addition operator for DAE meshes.
+        """
+
+        if isinstance(other, type(self)):
+            me = DAEMesh(self)
+            me.diff[:] = self.diff + other.diff
+            me.alg[:] = self.alg + other.alg
+            return me
+        else:
+            raise DataError("Type error: cannot add %s to %s" % (type(other), type(self)))
+
+    def __sub__(self, other):
+        """
+        Overloading the subtraction operator for DAE meshes.
+        """
+
+        if isinstance(other, type(self)):
+            me = DAEMesh(self)
+            me.diff[:] = self.diff - other.diff
+            me.alg[:] = self.alg - other.alg
+            return me
+        else:
+            raise DataError("Type error: cannot subtract %s to %s" % (type(other), type(self)))
+
+    def __rmul__(self, other):
+        """
+        Overloading the right multiply by factor operator for DAE meshes.
+        """
+
+        if isinstance(other, float):
+            me = DAEMesh(self)
+            me.diff[:] = other * self.diff
+            me.alg[:] = other * self.alg
+            return me
+        else:
+            raise DataError("Type error: cannot multiply %s to %s" % (type(other), type(self)))
+
+    def __str__(self):
+        """
+        Overloading the string operator for DAE Meshes s.t. an output can be separated into the
+        differential part and the algebraic part.
+        """
+        return f"({self.diff}, {self.alg})"

pySDC/projects/DAE/problems/DiscontinuousTestDAE.py

@@ -1,6 +1,8 @@
 import numpy as np
 
+from pySDC.core.Problem import WorkCounter
 from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae
+from pySDC.projects.DAE.misc.dae_mesh import DAEMesh
 
 
 class DiscontinuousTestDAE(ptype_dae):
@@ -57,14 +59,13 @@ class DiscontinuousTestDAE(ptype_dae):
 
     def __init__(self, newton_tol=1e-12):
         """Initialization routine"""
-        nvars = 2
-        super().__init__(nvars, newton_tol)
-        self._makeAttributeAndRegister('nvars', localVars=locals(), readOnly=True)
+        super().__init__(nvars=(1, 1), newton_tol=newton_tol)
         self._makeAttributeAndRegister('newton_tol', localVars=locals())
 
         self.t_switch_exact = np.arccosh(50)
         self.t_switch = None
         self.nswitches = 0
+        self.work_counters['rhs'] = WorkCounter()
 
     def eval_f(self, u, du, t):
         r"""
@@ -85,24 +86,21 @@ def eval_f(self, u, du, t):
             The right-hand side of f (contains two components).
         """
 
-        y, z = u[0], u[1]
-        dy = du[0]
+        y, z = u.diff, u.alg
+        dy = du.diff
 
         t_switch = np.inf if self.t_switch is None else self.t_switch
 
         h = 2 * y - 100
         f = self.dtype_f(self.init)
 
         if h >= 0 or t >= t_switch:
-            f[:] = (
-                dy,
-                y**2 - z**2 - 1,
-            )
+            f.diff[:] = dy
+            f.alg[:] = y**2 - z**2 - 1
         else:
-            f[:] = (
-                dy - z,
-                y**2 - z**2 - 1,
-            )
+            f.diff[:] = dy - z
+            f.alg[:] = y**2 - z**2 - 1
+        self.work_counters['rhs']()
         return f
 
     def u_exact(self, t, **kwargs):
@@ -125,9 +123,11 @@ def u_exact(self, t, **kwargs):
 
         me = self.dtype_u(self.init)
         if t <= self.t_switch_exact:
-            me[:] = (np.cosh(t), np.sinh(t))
+            me.diff[:] = np.cosh(t)
+            me.alg[:] = np.sinh(t)
         else:
-            me[:] = (np.cosh(self.t_switch_exact), np.sinh(self.t_switch_exact))
+            me.diff[:] = np.cosh(self.t_switch_exact)
+            me.alg[:] = np.sinh(self.t_switch_exact)
         return me
 
     def get_switching_info(self, u, t):
@@ -162,14 +162,14 @@ def get_switching_info(self, u, t):
         m_guess = -100
 
         for m in range(1, len(u)):
-            h_prev_node = 2 * u[m - 1][0] - 100
-            h_curr_node = 2 * u[m][0] - 100
+            h_prev_node = 2 * u[m - 1].diff - 100
+            h_curr_node = 2 * u[m].diff - 100
             if h_prev_node < 0 and h_curr_node >= 0:
                 switch_detected = True
                 m_guess = m - 1
                 break
 
-        state_function = [2 * u[m][0] - 100 for m in range(len(u))]
+        state_function = [2 * u[m].diff - 100 for m in range(len(u))]
         return switch_detected, m_guess, state_function
 
     def count_switches(self):