diff --git a/pySDC/projects/DAE/misc/DAEMesh.py b/pySDC/projects/DAE/misc/DAEMesh.py
new file mode 100644
index 0000000000..cf10e949c3
--- /dev/null
+++ b/pySDC/projects/DAE/misc/DAEMesh.py
@@ -0,0 +1,12 @@
+from pySDC.implementations.datatype_classes.mesh import MultiComponentMesh
+
+
+class DAEMesh(MultiComponentMesh):
+    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.
+    """
+
+    components = ['diff', 'alg']
diff --git a/pySDC/projects/DAE/misc/HookClass_DAE.py b/pySDC/projects/DAE/misc/HookClass_DAE.py
index 4ce9f12e57..df6d08e25f 100644
--- a/pySDC/projects/DAE/misc/HookClass_DAE.py
+++ b/pySDC/projects/DAE/misc/HookClass_DAE.py
@@ -1,77 +1,75 @@
 from pySDC.core.Hooks import hooks
 
 
-class approx_solution_hook(hooks):
+class LogGlobalErrorPostStepDifferentialVariable(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(approx_solution_hook, self).__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]
+        P = L.prob
 
         # TODO: is it really necessary to recompute the end point? Hasn't this been done already?
         L.sweep.compute_end_point()
 
+        # 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)
+        e_global_differential = abs(upde.diff - L.uend.diff)
+
         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='approx_solution',
-            value=L.uend,
+            type='e_global_differential_post_step',
+            value=e_global_differential,
         )
 
 
-class error_hook(hooks):
+class LogGlobalErrorPostStepAlgebraicVariable(hooks):
     """
-    Hook class to add the approximate solution to the output generated by the sweeper after each time step
+    Logs the global error in the algebraic variable
     """
 
-    def __init__(self):
-        """
-        Initialization routine for the custom hook
-        """
-        super(error_hook, self).__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]
         P = L.prob
 
-        # TODO: is it really necessary to recompute the end point? Hasn't this been done already?
         L.sweep.compute_end_point()
 
-        # 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])
+        e_global_algebraic = abs(upde.alg - L.uend.alg)
 
         self.add_to_stats(
             process=step.status.slot,
@@ -79,6 +77,6 @@ def post_step(self, step, level_number):
             level=L.level_index,
             iter=step.status.iter,
             sweep=L.status.sweep,
-            type='error_post_step',
-            value=err,
+            type='e_global_algebraic_post_step',
+            value=e_global_algebraic,
         )
diff --git a/pySDC/projects/DAE/misc/ProblemDAE.py b/pySDC/projects/DAE/misc/ProblemDAE.py
index 08cc16a3a7..931d466b79 100644
--- a/pySDC/projects/DAE/misc/ProblemDAE.py
+++ b/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.DAEMesh 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,18 @@ def solve_system(self, impl_sys, u0, t):
         me : dtype_u
             Numerical solution.
         """
-
         me = self.dtype_u(self.init)
+
+        def implSysFlatten(unknowns, **kwargs):
+            sys = impl_sys(unknowns.reshape(me.shape).view(type(u0)), **kwargs)
+            return sys.flatten()
+
         opt = root(
-            impl_sys,
+            implSysFlatten,
             u0,
             method='hybr',
             tol=self.newton_tol,
         )
-        me[:] = opt.x
+        me[:] = opt.x.reshape(me.shape)
         self.work_counters['newton'].niter += opt.nfev
         return me
diff --git a/pySDC/projects/DAE/problems/DiscontinuousTestDAE.py b/pySDC/projects/DAE/problems/DiscontinuousTestDAE.py
index 8df106a828..3b676d8cec 100644
--- a/pySDC/projects/DAE/problems/DiscontinuousTestDAE.py
+++ b/pySDC/projects/DAE/problems/DiscontinuousTestDAE.py
@@ -1,5 +1,6 @@
 import numpy as np
 
+from pySDC.core.Problem import WorkCounter
 from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae
 
 
@@ -57,14 +58,12 @@ 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)
-        self._makeAttributeAndRegister('newton_tol', localVars=locals())
+        super().__init__(nvars=2, newton_tol=newton_tol)
 
         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,8 +84,8 @@ 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[0], u.alg[0]
+        dy = du.diff[0]
 
         t_switch = np.inf if self.t_switch is None else self.t_switch
 
@@ -94,15 +93,12 @@ def eval_f(self, u, du, t):
         f = self.dtype_f(self.init)
 
         if h >= 0 or t >= t_switch:
-            f[:] = (
-                dy,
-                y**2 - z**2 - 1,
-            )
+            f.diff[0] = dy
+            f.alg[0] = y**2 - z**2 - 1
         else:
-            f[:] = (
-                dy - z,
-                y**2 - z**2 - 1,
-            )
+            f.diff[0] = dy - z
+            f.alg[0] = y**2 - z**2 - 1
+        self.work_counters['rhs']()
         return f
 
     def u_exact(self, t, **kwargs):
@@ -125,9 +121,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[0] = np.cosh(t)
+            me.alg[0] = np.sinh(t)
         else:
-            me[:] = (np.cosh(self.t_switch_exact), np.sinh(self.t_switch_exact))
+            me.diff[0] = np.cosh(self.t_switch_exact)
+            me.alg[0] = np.sinh(self.t_switch_exact)
         return me
 
     def get_switching_info(self, u, t):
@@ -162,14 +160,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[0] - 100
+            h_curr_node = 2 * u[m].diff[0] - 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[0] - 100 for m in range(len(u))]
         return switch_detected, m_guess, state_function
 
     def count_switches(self):
diff --git a/pySDC/projects/DAE/problems/WSCC9BusSystem.py b/pySDC/projects/DAE/problems/WSCC9BusSystem.py
index ad472dc636..3867f0443a 100644
--- a/pySDC/projects/DAE/problems/WSCC9BusSystem.py
+++ b/pySDC/projects/DAE/problems/WSCC9BusSystem.py
@@ -1,6 +1,5 @@
 import numpy as np
 from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae
-from pySDC.implementations.datatype_classes.mesh import mesh
 from pySDC.core.Errors import ParameterError
 
 
@@ -763,16 +762,12 @@ class WSCC9BusSystem(ptype_dae):
        for Power Systems Research and Education. IEEE Transactions on Power Systems. Vol. 26, No. 1, pp. 12–19 (2011).
     """
 
-    dtype_u = mesh
-    dtype_f = mesh
-
-    def __init__(self, nvars=None, newton_tol=1e-10, m=3, n=9):
+    def __init__(self, newton_tol=1e-10):
         """Initialization routine"""
-
+        m, n = 3, 9
         nvars = 11 * m + 2 * m + 2 * n
         # invoke super init, passing number of dofs
-        super().__init__(nvars, newton_tol)
-        self._makeAttributeAndRegister('nvars', 'newton_tol', localVars=locals(), readOnly=True)
+        super().__init__(nvars=nvars, newton_tol=newton_tol)
         self._makeAttributeAndRegister('m', 'n', localVars=locals())
         self.mpc = WSCC9Bus()
 
@@ -993,19 +988,31 @@ def eval_f(self, u, du, t):
             The right-hand side of f (contains two components).
         """
 
-        dEqp, dSi1d, dEdp = du[0 : self.m], du[self.m : 2 * self.m], du[2 * self.m : 3 * self.m]
-        dSi2q, dDelta = du[3 * self.m : 4 * self.m], du[4 * self.m : 5 * self.m]
-        dw, dEfd, dRF = du[5 * self.m : 6 * self.m], du[6 * self.m : 7 * self.m], du[7 * self.m : 8 * self.m]
-        dVR, dTM, dPSV = du[8 * self.m : 9 * self.m], du[9 * self.m : 10 * self.m], du[10 * self.m : 11 * self.m]
+        dEqp, dSi1d, dEdp = du.diff[0 : self.m], du.diff[self.m : 2 * self.m], du.diff[2 * self.m : 3 * self.m]
+        dSi2q, dDelta = du.diff[3 * self.m : 4 * self.m], du.diff[4 * self.m : 5 * self.m]
+        dw, dEfd, dRF = (
+            du.diff[5 * self.m : 6 * self.m],
+            du.diff[6 * self.m : 7 * self.m],
+            du.diff[7 * self.m : 8 * self.m],
+        )
+        dVR, dTM, dPSV = (
+            du.diff[8 * self.m : 9 * self.m],
+            du.diff[9 * self.m : 10 * self.m],
+            du.diff[10 * self.m : 11 * self.m],
+        )
 
-        Eqp, Si1d, Edp = u[0 : self.m], u[self.m : 2 * self.m], u[2 * self.m : 3 * self.m]
-        Si2q, Delta = u[3 * self.m : 4 * self.m], u[4 * self.m : 5 * self.m]
-        w, Efd, RF = u[5 * self.m : 6 * self.m], u[6 * self.m : 7 * self.m], u[7 * self.m : 8 * self.m]
-        VR, TM, PSV = u[8 * self.m : 9 * self.m], u[9 * self.m : 10 * self.m], u[10 * self.m : 11 * self.m]
+        Eqp, Si1d, Edp = u.diff[0 : self.m], u.diff[self.m : 2 * self.m], u.diff[2 * self.m : 3 * self.m]
+        Si2q, Delta = u.diff[3 * self.m : 4 * self.m], u.diff[4 * self.m : 5 * self.m]
+        w, Efd, RF = u.diff[5 * self.m : 6 * self.m], u.diff[6 * self.m : 7 * self.m], u.diff[7 * self.m : 8 * self.m]
+        VR, TM, PSV = (
+            u.diff[8 * self.m : 9 * self.m],
+            u.diff[9 * self.m : 10 * self.m],
+            u.diff[10 * self.m : 11 * self.m],
+        )
 
-        Id, Iq = u[11 * self.m : 11 * self.m + self.m], u[11 * self.m + self.m : 11 * self.m + 2 * self.m]
-        V = u[11 * self.m + 2 * self.m : 11 * self.m + 2 * self.m + self.n]
-        TH = u[11 * self.m + 2 * self.m + self.n : 11 * self.m + 2 * self.m + 2 * self.n]
+        Id, Iq = u.alg[0 : self.m], u.alg[self.m : 2 * self.m]
+        V = u.alg[2 * self.m : 2 * self.m + self.n]
+        TH = u.alg[2 * self.m + self.n : 2 * self.m + 2 * self.n]
 
         # line outage disturbance:
         if t >= 0.05:
@@ -1137,7 +1144,8 @@ def eval_f(self, u, du, t):
         eqs.append(-QL2[self.m : self.n] - sum4)  # (17)
         eqs_flatten = [item for sublist in eqs for item in sublist]
 
-        f[:] = eqs_flatten
+        f.diff[: 11 * self.m] = eqs_flatten[0 : 11 * self.m]
+        f.alg[: 2 * self.n + 2 * self.m] = eqs_flatten[11 * self.m :]
         return f
 
     def u_exact(self, t):
@@ -1157,24 +1165,24 @@ def u_exact(self, t):
         assert t == 0, 'ERROR: u_exact only valid for t=0'
 
         me = self.dtype_u(self.init)
-        me[0 : self.m] = self.Eqp0
-        me[self.m : 2 * self.m] = self.Si1d0
-        me[2 * self.m : 3 * self.m] = self.Edp0
-        me[3 * self.m : 4 * self.m] = self.Si2q0
-        me[4 * self.m : 5 * self.m] = self.D0
-        me[5 * self.m : 6 * self.m] = self.ws_vector
-        me[6 * self.m : 7 * self.m] = self.Efd0
-        me[7 * self.m : 8 * self.m] = self.RF0
-        me[8 * self.m : 9 * self.m] = self.VR0
-        me[9 * self.m : 10 * self.m] = self.TM0
-        me[10 * self.m : 11 * self.m] = self.PSV0
-        me[11 * self.m : 11 * self.m + self.m] = self.Id0
-        me[11 * self.m + self.m : 11 * self.m + 2 * self.m] = self.Iq0
-        me[11 * self.m + 2 * self.m : 11 * self.m + 2 * self.m + self.n] = self.V0
-        me[11 * self.m + 2 * self.m + self.n : 11 * self.m + 2 * self.m + 2 * self.n] = self.TH0
+        me.diff[0 : self.m] = self.Eqp0
+        me.diff[self.m : 2 * self.m] = self.Si1d0
+        me.diff[2 * self.m : 3 * self.m] = self.Edp0
+        me.diff[3 * self.m : 4 * self.m] = self.Si2q0
+        me.diff[4 * self.m : 5 * self.m] = self.D0
+        me.diff[5 * self.m : 6 * self.m] = self.ws_vector
+        me.diff[6 * self.m : 7 * self.m] = self.Efd0
+        me.diff[7 * self.m : 8 * self.m] = self.RF0
+        me.diff[8 * self.m : 9 * self.m] = self.VR0
+        me.diff[9 * self.m : 10 * self.m] = self.TM0
+        me.diff[10 * self.m : 11 * self.m] = self.PSV0
+        me.alg[0 : self.m] = self.Id0
+        me.alg[self.m : 2 * self.m] = self.Iq0
+        me.alg[2 * self.m : 2 * self.m + self.n] = self.V0
+        me.alg[2 * self.m + self.n : 2 * self.m + 2 * self.n] = self.TH0
         return me
 
-    def get_switching_info(self, u, t, du=None):
+    def get_switching_info(self, u, t):
         r"""
         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
@@ -1208,14 +1216,14 @@ def get_switching_info(self, u, t, du=None):
         switch_detected = False
         m_guess = -100
         for m in range(1, len(u)):
-            h_prev_node = u[m - 1][10 * self.m] - self.psv_max
-            h_curr_node = u[m][10 * self.m] - self.psv_max
+            h_prev_node = u[m - 1].diff[10 * self.m] - self.psv_max
+            h_curr_node = u[m].diff[10 * self.m] - self.psv_max
             if h_prev_node < 0 and h_curr_node >= 0:
                 switch_detected = True
                 m_guess = m - 1
                 break
 
-        state_function = [u[m][10 * self.m] - self.psv_max for m in range(len(u))]
+        state_function = [u[m].diff[10 * self.m] - self.psv_max for m in range(len(u))]
         return switch_detected, m_guess, state_function
 
     def count_switches(self):
diff --git a/pySDC/projects/DAE/problems/simple_DAE.py b/pySDC/projects/DAE/problems/simple_DAE.py
index 343128624a..31be74e809 100644
--- a/pySDC/projects/DAE/problems/simple_DAE.py
+++ b/pySDC/projects/DAE/problems/simple_DAE.py
@@ -3,6 +3,7 @@
 from scipy.interpolate import interp1d
 
 from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae
+from pySDC.implementations.datatype_classes.mesh import mesh
 
 
 class pendulum_2d(ptype_dae):
@@ -11,15 +12,22 @@ class pendulum_2d(ptype_dae):
     The DAE system is given by the equations
 
     .. math::
-        x' = u,
+        \frac{dp}{dt} = u,
 
     .. math::
-        \frac{d}{dt} \frac{\partial}{\partial u} L = \frac{\partial L}{\partial x} + f + G^{T} \lambda,
+        \frac{dq}{dt} = v,
 
     .. math::
-        0 = \phi.
+        m\frac{du}{dt} = -p \lambda,
 
-    The pendulum is used in most introductory literature on DAEs, for example on page 8 of [1]_.
+    .. math::
+        m\frac{dv}{dt} = -q \lambda - g,
+
+    .. math::
+        0 = p^2 + q^2 - l^2
+
+    for :math:`l=1` and :math:`m=1`. The pendulum is used in most introductory literature on DAEs, for example on page 8
+    of [1]_.
 
     Parameters
     ----------
@@ -39,9 +47,9 @@ class pendulum_2d(ptype_dae):
         Lect. Notes Math. (1989).
     """
 
-    def __init__(self, nvars, newton_tol):
+    def __init__(self, newton_tol):
         """Initialization routine"""
-        super().__init__(nvars, newton_tol)
+        super().__init__(nvars=5, newton_tol=newton_tol)
         # load reference solution
         # data file must be generated and stored under misc/data and self.t_end = t[-1]
         # data = np.load(r'pySDC/projects/DAE/misc/data/pendulum.npy')
@@ -72,7 +80,13 @@ def eval_f(self, u, du, t):
         # The last element of u is a Lagrange multiplier. Not sure if this needs to be time dependent, but must model the
         # weight somehow
         f = self.dtype_f(self.init)
-        f[:] = (du[0] - u[2], du[1] - u[3], du[2] + u[4] * u[0], du[3] + u[4] * u[1] + g, u[0] ** 2 + u[1] ** 2 - 1)
+        f.diff[:4] = (
+            du.diff[0] - u.diff[2],
+            du.diff[1] - u.diff[3],
+            du.diff[2] + u.alg[0] * u.diff[0],
+            du.diff[3] + u.alg[0] * u.diff[1] + g,
+        )
+        f.alg[0] = u.diff[0] ** 2 + u.diff[1] ** 2 - 1
         self.work_counters['rhs']()
         return f
 
@@ -92,12 +106,16 @@ def u_exact(self, t):
         """
         me = self.dtype_u(self.init)
         if t == 0:
-            me[:] = (-1, 0, 0, 0, 0)
+            me.diff[:4] = (-1, 0, 0, 0)
+            me.alg[0] = 0
         elif t < self.t_end:
-            me[:] = self.u_ref(t)
+            u_ref = self.u_ref(t)
+            me.diff[:4] = u_ref[:4]
+            me.alg[0] = u_ref[5]
         else:
             self.logger.warning("Requested time exceeds domain of the reference solution. Returning zero.")
-            me[:] = (0, 0, 0, 0, 0)
+            me.diff[:4] = (0, 0, 0, 0)
+            me.alg[0] = 0
         return me
 
 
@@ -139,6 +157,10 @@ class simple_dae_1(ptype_dae):
         equations. Society for Industrial and Applied Mathematics (1998).
     """
 
+    def __init__(self, newton_tol=1e-10):
+        """Initialization routine"""
+        super().__init__(nvars=3, newton_tol=newton_tol)
+
     def eval_f(self, u, du, t):
         r"""
         Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`.
@@ -160,11 +182,12 @@ def eval_f(self, u, du, t):
         # Smooth index-2 DAE pg. 267 Ascher and Petzold (also the first example in KDC Minion paper)
         a = 10.0
         f = self.dtype_f(self.init)
-        f[:] = (
-            -du[0] + (a - 1 / (2 - t)) * u[0] + (2 - t) * a * u[2] + np.exp(t) * (3 - t) / (2 - t),
-            -du[1] + (1 - a) / (t - 2) * u[0] - u[1] + (a - 1) * u[2] + 2 * np.exp(t),
-            (t + 2) * u[0] + (t**2 - 4) * u[1] - (t**2 + t - 2) * np.exp(t),
+
+        f.diff[:2] = (
+            -du.diff[0] + (a - 1 / (2 - t)) * u.diff[0] + (2 - t) * a * u.alg[0] + (3 - t) / (2 - t) * np.exp(t),
+            -du.diff[1] + (1 - a) / (t - 2) * u.diff[0] - u.diff[1] + (a - 1) * u.alg[0] + 2 * np.exp(t),
         )
+        f.alg[0] = (t + 2) * u.diff[0] + (t**2 - 4) * u.diff[1] - (t**2 + t - 2) * np.exp(t)
         self.work_counters['rhs']()
         return f
 
@@ -183,7 +206,8 @@ def u_exact(self, t):
             The reference solution as mesh object containing three components.
         """
         me = self.dtype_u(self.init)
-        me[:] = (np.exp(t), np.exp(t), -np.exp(t) / (2 - t))
+        me.diff[:2] = (np.exp(t), np.exp(t))
+        me.alg[0] = -np.exp(t) / (2 - t)
         return me
 
 
@@ -193,10 +217,10 @@ class problematic_f(ptype_dae):
     numerically solvable for certain choices of the parameter :math:`\eta`. The DAE system is given by
 
     .. math::
-        y (t) + \eta t z (t) = f(t),
+        \frac{d y(t)}{dt} + \eta t \frac{d z(t)}{dt} + (1 + \eta) z (t) = g (t).
 
     .. math::
-        \frac{d y(t)}{dt} + \eta t \frac{d z(t)}{dt} + (1 + \eta) z (t) = g (t).
+        y (t) + \eta t z (t) = f(t),
 
     See, for example, page 264 of [1]_.
 
@@ -218,9 +242,12 @@ class problematic_f(ptype_dae):
         equations. Society for Industrial and Applied Mathematics (1998).
     """
 
-    def __init__(self, nvars, newton_tol, eta=1):
+    dtype_u = mesh
+    dtype_f = mesh
+
+    def __init__(self, newton_tol, eta=1):
         """Initialization routine"""
-        super().__init__(nvars, newton_tol)
+        super().__init__(nvars=2, newton_tol=newton_tol)
         self._makeAttributeAndRegister('eta', localVars=locals())
 
     def eval_f(self, u, du, t):
diff --git a/pySDC/projects/DAE/problems/synchronous_machine.py b/pySDC/projects/DAE/problems/synchronous_machine.py
index 5510f0dd46..9042325303 100644
--- a/pySDC/projects/DAE/problems/synchronous_machine.py
+++ b/pySDC/projects/DAE/problems/synchronous_machine.py
@@ -149,8 +149,8 @@ class synchronous_machine_infinite_bus(ptype_dae):
     .. [1] P. Kundur, N. J. Balu, M. G. Lauby. Power system stability and control. The EPRI power system series (1994).
     """
 
-    def __init__(self, nvars, newton_tol):
-        super(synchronous_machine_infinite_bus, self).__init__(nvars, newton_tol)
+    def __init__(self, newton_tol):
+        super().__init__(nvars=14, newton_tol=newton_tol)
         # load reference solution
         # data file must be generated and stored under misc/data and self.t_end = t[-1]
         # data = np.load(r'pySDC/projects/DAE/misc/data/synch_gen.npy')
@@ -217,16 +217,22 @@ def eval_f(self, u, du, t):
 
         # extract variables for readability
         # algebraic components
-        psi_d, psi_q, psi_F, psi_D, psi_Q1, psi_Q2 = u[0], u[1], u[2], u[3], u[4], u[5]
-        i_d, i_q, i_F, i_D, i_Q1, i_Q2 = u[6], u[7], u[8], u[9], u[10], u[11]
-        delta_r = u[12]
-        omega_m = u[13]
+        psi_d, psi_q, psi_F, psi_D, psi_Q1, psi_Q2 = u.diff[0], u.diff[1], u.diff[2], u.diff[3], u.diff[4], u.diff[5]
+        delta_r, omega_m = u.diff[6], u.diff[7]
+        i_d, i_q, i_F, i_D, i_Q1, i_Q2 = u.alg[0], u.alg[1], u.alg[2], u.alg[3], u.alg[4], u.alg[5]
 
         # differential components
         # these result directly from the voltage equations, introduced e.g. pg. 145 Krause
-        dpsi_d, dpsi_q, dpsi_F, dpsi_D, dpsi_Q1, dpsi_Q2 = du[0], du[1], du[2], du[3], du[4], du[5]
-        ddelta_r = du[12]
-        domega_m = du[13]
+        dpsi_d, dpsi_q, dpsi_F, dpsi_D, dpsi_Q1, dpsi_Q2 = (
+            du.diff[0],
+            du.diff[1],
+            du.diff[2],
+            du.diff[3],
+            du.diff[4],
+            du.diff[5],
+        )
+        ddelta_r, domega_m = du.diff[6], du.diff[7]
+
         # Network current
         I_Re = i_d * np.sin(delta_r) + i_q * np.cos(delta_r)
         I_Im = -i_d * np.cos(delta_r) + i_q * np.sin(delta_r)
@@ -237,9 +243,7 @@ def eval_f(self, u, du, t):
         v_d = np.real(V_comp) * np.sin(delta_r) - np.imag(V_comp) * np.cos(delta_r)
         v_q = np.real(V_comp) * np.cos(delta_r) + np.imag(V_comp) * np.sin(delta_r)
 
-        # algebraic variables are i_d, i_q, i_F, i_D, i_Q1, i_Q2, il_d, il_q
-        f[:] = (
-            # differential generator
+        f.diff[:8] = (
             -dpsi_d + self.omega_b * (v_d - self.R_s * i_d + omega_m * psi_q),
             -dpsi_q + self.omega_b * (v_q - self.R_s * i_q - omega_m * psi_d),
             -dpsi_F + self.omega_b * (self.v_F - self.R_F * i_F),
@@ -249,7 +253,8 @@ def eval_f(self, u, du, t):
             -ddelta_r + self.omega_b * (omega_m - 1),
             -domega_m
             + 1 / (2 * self.H_) * (self.T_m - (psi_q * i_d - psi_d * i_q) - self.K_D * self.omega_b * (omega_m - 1)),
-            # algebraic generator
+        )
+        f.alg[:6] = (
             -psi_d + self.L_d * i_d + self.L_md * i_F + self.L_md * i_D,
             -psi_q + self.L_q * i_q + self.L_mq * i_Q1 + self.L_mq * i_Q2,
             -psi_F + self.L_md * i_d + self.L_F * i_F + self.L_md * i_D,
@@ -296,12 +301,16 @@ def u_exact(self, t):
             omega_b = 2 * np.pi * 60
             omega_m = omega_0 / omega_b  # = omega_r since pf = 2 i.e. two pole machine
 
-            me[:] = (psi_d, psi_q, psi_F, psi_D, psi_Q1, psi_Q2, i_d, i_q, i_F, i_D, i_Q1, i_Q2, delta_r, omega_m)
+            me.diff[:8] = (psi_d, psi_q, psi_F, psi_D, psi_Q1, psi_Q2, delta_r, omega_m)
+            me.alg[:6] = (i_d, i_q, i_F, i_D, i_Q1, i_Q2)
         elif t < self.t_end:
-            me[:] = self.u_ref(t)
+            u_ref = self.u_ref(t)
+            me.diff[:8] = u_ref[:8]
+            me.alg[:6] = u_ref[8:]
         else:
             self.logger.warning("Requested time exceeds domain of the reference solution. Returning zero.")
-            me[:] = (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
+            me.diff[:8] = (0, 0, 0, 0, 0, 0, 0, 0)
+            me.alg[:6] = (0, 0, 0, 0, 0, 0)
 
         return me
 
diff --git a/pySDC/projects/DAE/problems/transistor_amplifier.py b/pySDC/projects/DAE/problems/transistor_amplifier.py
index ae2adeb9cc..bc0bb9cf88 100644
--- a/pySDC/projects/DAE/problems/transistor_amplifier.py
+++ b/pySDC/projects/DAE/problems/transistor_amplifier.py
@@ -3,6 +3,7 @@
 from scipy.interpolate import interp1d
 
 from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae
+from pySDC.implementations.datatype_classes.mesh import mesh
 
 
 # Helper function
@@ -12,7 +13,7 @@ def _transistor(u_in):
 
 class one_transistor_amplifier(ptype_dae):
     r"""
-    The one transistor amplifier example from pg. 404 in [1]_. The problem is an index-1 differential-algebraic equation
+    The one transistor amplifier example from pg. 377 in [1]_. The problem is an index-1 differential-algebraic equation
     (DAE) having the equations
 
     .. math::
@@ -61,8 +62,11 @@ class one_transistor_amplifier(ptype_dae):
         Springer (2009).
     """
 
-    def __init__(self, nvars, newton_tol):
-        super().__init__(nvars, newton_tol)
+    dtype_u = mesh
+    dtype_f = mesh
+
+    def __init__(self, newton_tol):
+        super().__init__(nvars=5, newton_tol=newton_tol)
         # load reference solution
         # data file must be generated and stored under misc/data and self.t_end = t[-1]
         # data = np.load(r'pySDC/projects/DAE/misc/data/one_trans_amp.npy')
@@ -195,8 +199,11 @@ class two_transistor_amplifier(ptype_dae):
         Lect. Notes Math. (1989).
     """
 
-    def __init__(self, nvars, newton_tol):
-        super().__init__(nvars, newton_tol)
+    dtype_u = mesh
+    dtype_f = mesh
+
+    def __init__(self, newton_tol):
+        super().__init__(nvars=8, newton_tol=newton_tol)
         # load reference solution
         # data file must be generated and stored under misc/data and self.t_end = t[-1]
         # data = np.load(r'pySDC/projects/DAE/misc/data/two_trans_amp.npy')
diff --git a/pySDC/projects/DAE/run/fully_implicit_dae_playground.py b/pySDC/projects/DAE/run/fully_implicit_dae_playground.py
index e4f0918855..72c50311b5 100644
--- a/pySDC/projects/DAE/run/fully_implicit_dae_playground.py
+++ b/pySDC/projects/DAE/run/fully_implicit_dae_playground.py
@@ -5,9 +5,9 @@
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
 from pySDC.projects.DAE.problems.simple_DAE import problematic_f
 from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE
-from pySDC.projects.DAE.misc.HookClass_DAE import approx_solution_hook
-from pySDC.projects.DAE.misc.HookClass_DAE import error_hook
+from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostStep
 from pySDC.helpers.stats_helper import get_sorted
+from pySDC.implementations.hooks.log_solution import LogSolution
 
 
 def main():
@@ -26,8 +26,7 @@ def main():
 
     # initialize problem parameters
     problem_params = dict()
-    problem_params['newton_tol'] = 1e-12  # tollerance for implicit solver
-    problem_params['nvars'] = 2
+    problem_params['newton_tol'] = 1e-12
 
     # initialize step parameters
     step_params = dict()
@@ -36,7 +35,7 @@ def main():
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 30
-    controller_params['hook_class'] = [approx_solution_hook, error_hook]
+    controller_params['hook_class'] = [LogSolution, LogGlobalErrorPostStep]
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -64,15 +63,15 @@ def main():
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
     # check error
-    err = get_sorted(stats, type='error_post_step', sortby='time')
+    err = get_sorted(stats, type='e_global_post_step', sortby='time')
     err = np.linalg.norm([err[i][1] for i in range(len(err))], np.inf)
     print(f"Error is {err}")
     assert np.isclose(err, 0.0, atol=1e-4), "Error too large."
 
     # store results
-    sol = get_sorted(stats, type='approx_solution_hook', sortby='time')
+    sol = get_sorted(stats, type='u', sortby='time')
     sol_dt = np.array([sol[i][0] for i in range(len(sol))])
-    sol_data = np.array([[sol[j][1][i] for j in range(len(sol))] for i in range(problem_params['nvars'])])
+    sol_data = np.array([[sol[j][1][i] for j in range(len(sol))] for i in range(P.nvars)])
 
     data = dict()
     data['dt'] = sol_dt
diff --git a/pySDC/projects/DAE/run/run_convergence_test.py b/pySDC/projects/DAE/run/run_convergence_test.py
index 3229d4f7f7..b8a23671fc 100644
--- a/pySDC/projects/DAE/run/run_convergence_test.py
+++ b/pySDC/projects/DAE/run/run_convergence_test.py
@@ -5,7 +5,7 @@
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
 from pySDC.projects.DAE.problems.simple_DAE import simple_dae_1
 from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE
-from pySDC.projects.DAE.misc.HookClass_DAE import error_hook
+from pySDC.projects.DAE.misc.HookClass_DAE import LogGlobalErrorPostStepDifferentialVariable
 from pySDC.helpers.stats_helper import get_sorted
 from pySDC.helpers.stats_helper import filter_stats
 
@@ -25,7 +25,6 @@ def setup():
     # This comes as read-in for the problem class
     problem_params = dict()
     problem_params['newton_tol'] = 1e-3  # tollerance for implicit solver
-    problem_params['nvars'] = 3
 
     # This comes as read-in for the step class
     step_params = dict()
@@ -34,7 +33,7 @@ def setup():
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 30
-    controller_params['hook_class'] = error_hook
+    controller_params['hook_class'] = LogGlobalErrorPostStepDifferentialVariable
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -90,7 +89,7 @@ def run(description, controller_params, run_params):
                 uend, stats = controller.run(u0=uinit, t0=run_params['t0'], Tend=run_params['tend'])
 
                 # compute exact solution and compare
-                err = get_sorted(stats, type='error_post_step', sortby='time')
+                err = get_sorted(stats, type='e_global_differential_post_step', sortby='time')
                 niter = filter_stats(stats, type='niter')
 
                 conv_data[qd_type][num_nodes]['error'][j] = np.linalg.norm([err[j][1] for j in range(len(err))], np.inf)
diff --git a/pySDC/projects/DAE/run/run_iteration_test.py b/pySDC/projects/DAE/run/run_iteration_test.py
index 986895b044..bf6c0d1fdb 100644
--- a/pySDC/projects/DAE/run/run_iteration_test.py
+++ b/pySDC/projects/DAE/run/run_iteration_test.py
@@ -6,7 +6,7 @@
 from pySDC.projects.DAE.problems.simple_DAE import simple_dae_1
 from pySDC.projects.DAE.problems.transistor_amplifier import one_transistor_amplifier
 from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE
-from pySDC.projects.DAE.misc.HookClass_DAE import error_hook
+from pySDC.projects.DAE.misc.HookClass_DAE import LogGlobalErrorPostStepDifferentialVariable
 from pySDC.helpers.stats_helper import get_sorted
 from pySDC.helpers.stats_helper import filter_stats
 
@@ -32,7 +32,6 @@ def setup():
     # Absolute termination tollerance for implicit solver
     # Exactly how this is used can be adjusted in update_nodes() in the fully implicit sweeper
     problem_params['newton_tol'] = 1e-7
-    problem_params['nvars'] = 3
 
     # This comes as read-in for the step class
     step_params = dict()
@@ -40,7 +39,7 @@ def setup():
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 30
-    controller_params['hook_class'] = error_hook
+    controller_params['hook_class'] = LogGlobalErrorPostStepDifferentialVariable
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -98,7 +97,7 @@ def run(description, controller_params, run_params):
                 uend, stats = controller.run(u0=uinit, t0=run_params['t0'], Tend=run_params['tend'])
 
                 # compute exact solution and compare
-                err = get_sorted(stats, type='error_post_step', sortby='time')
+                err = get_sorted(stats, type='e_global_differential_post_step', sortby='time')
                 residual = get_sorted(stats, type='residual_post_step', sortby='time')
                 niter = filter_stats(stats, type='niter')
 
diff --git a/pySDC/projects/DAE/run/synchronous_machine_playground.py b/pySDC/projects/DAE/run/synchronous_machine_playground.py
index 32e8b5a4fc..f658966224 100644
--- a/pySDC/projects/DAE/run/synchronous_machine_playground.py
+++ b/pySDC/projects/DAE/run/synchronous_machine_playground.py
@@ -6,10 +6,10 @@
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
 from pySDC.projects.DAE.problems.synchronous_machine import synchronous_machine_infinite_bus
 from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE
-from pySDC.projects.DAE.misc.HookClass_DAE import approx_solution_hook
-from pySDC.projects.DAE.misc.HookClass_DAE import error_hook
+from pySDC.projects.DAE.misc.HookClass_DAE import LogGlobalErrorPostStepDifferentialVariable
 from pySDC.helpers.stats_helper import get_sorted
 from pySDC.helpers.stats_helper import filter_stats
+from pySDC.implementations.hooks.log_solution import LogSolution
 
 
 def main():
@@ -29,8 +29,7 @@ def main():
 
     # initialize problem parameters
     problem_params = dict()
-    problem_params['newton_tol'] = 1e-3  # tollerance for implicit solver
-    problem_params['nvars'] = 14
+    problem_params['newton_tol'] = 1e-3
 
     # initialize step parameters
     step_params = dict()
@@ -39,7 +38,7 @@ def main():
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 30
-    controller_params['hook_class'] = [error_hook, approx_solution_hook]
+    controller_params['hook_class'] = [LogGlobalErrorPostStepDifferentialVariable, LogSolution]
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -67,33 +66,41 @@ def main():
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
     # check error (only available if reference solution was provided)
-    # err = get_sorted(stats, type='error_post_step', sortby='time')
+    # err = get_sorted(stats, type='e_global_differential_post_step', sortby='time')
     # err = np.linalg.norm([err[i][1] for i in range(len(err))], np.inf)
     # print(f"Error is {err}")
 
-    uend_ref = [
+    uend_ref = P.dtype_u(P.init)
+    uend_ref.diff[:8] = (
         8.30823565e-01,
         -4.02584174e-01,
         1.16966755e00,
         9.47592808e-01,
         -3.68076863e-01,
         -3.87492326e-01,
+        3.10281509e-01,
+        9.94039645e-01,
+    )
+    uend_ref.alg[:6] = (
         -7.77837831e-01,
         -1.67347611e-01,
         1.34810867e00,
         5.46223705e-04,
         1.29690691e-02,
         -8.00823474e-02,
-        3.10281509e-01,
-        9.94039645e-01,
-    ]
-    err = np.linalg.norm(uend - uend_ref, np.inf)
+    )
+    err = abs(uend.diff - uend_ref.diff)
     assert np.isclose(err, 0, atol=1e-4), "Error too large."
 
     # store results
-    sol = get_sorted(stats, type='approx_solution', sortby='time')
+    sol = get_sorted(stats, type='u', sortby='time')
     sol_dt = np.array([sol[i][0] for i in range(len(sol))])
-    sol_data = np.array([[sol[j][1][i] for j in range(len(sol))] for i in range(problem_params['nvars'])])
+    sol_data = np.array(
+        [
+            [(sol[j][1].diff[id], sol[j][1].alg[ia]) for j in range(len(sol))]
+            for id, ia in zip(range(len(uend.diff)), range(len(uend.alg)))
+        ]
+    )
     niter = filter_stats(stats, type='niter')
     niter = np.fromiter(niter.values(), int)
 
diff --git a/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py b/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py
index aab811ffb9..8efeff6263 100644
--- a/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py
+++ b/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py
@@ -3,6 +3,7 @@
 
 from pySDC.core.Errors import ParameterError
 from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+from pySDC.projects.DAE.misc.DAEMesh import DAEMesh
 
 
 class fully_implicit_DAE(generic_implicit):
@@ -46,18 +47,13 @@ class fully_implicit_DAE(generic_implicit):
     """
 
     def __init__(self, params):
-        """
-        Initialization routine for the custom sweeper
-
-        Args:
-            params: parameters for the sweeper
-        """
+        """Initialization routine"""
 
         if 'QI' not in params:
             params['QI'] = 'IE'
 
         # call parent's initialization routine
-        super(fully_implicit_DAE, self).__init__(params)
+        super().__init__(params)
 
         msg = f"Quadrature type {self.params.quad_type} is not implemented yet. Use 'RADAU-RIGHT' instead!"
         if self.coll.left_is_node:
@@ -71,34 +67,27 @@ def update_nodes(self):
         preconditioned Richardson iteration in **"ordinary"** SDC.
         """
 
-        # get current level and problem description
         L = self.level
-        # in the fully implicit case L.prob.eval_f() evaluates the function F(u, u', t)
         P = L.prob
 
         # only if the level has been touched before
         assert L.status.unlocked
 
-        # get number of collocation nodes for easier access
         M = self.coll.num_nodes
         u_0 = L.u[0]
 
         # get QU^k where U = u'
-        # note that for multidimensional functions the required Kronecker product is achieved since
-        # e.g. L.f[j] is a mesh object and multiplication with a number distributes over the mesh
         integral = self.integrate()
         # build the rest of the known solution u_0 + del_t(Q - Q_del)U_k
         for m in range(1, M + 1):
             for j in range(1, M + 1):
                 integral[m - 1] -= L.dt * self.QI[m, j] * L.f[j]
-            # add initial value
             integral[m - 1] += u_0
 
         # do the sweep
         for m in range(1, M + 1):
-            # build implicit function, consisting of the known values from above and new values from previous nodes (at k+1)
-            u_approx = P.dtype_u(integral[m - 1])
             # add the known components from current sweep del_t*Q_del*U_k+1
+            u_approx = P.dtype_u(integral[m - 1])
             for j in range(1, m):
                 u_approx += L.dt * self.QI[m, j] * L.f[j]
 
@@ -118,11 +107,10 @@ def implSystem(params):
                     System to be solved as implicit function.
                 """
 
-                params_mesh = P.dtype_f(P.init)
-                params_mesh[:] = params
+                params_mesh = P.dtype_f(params)
 
                 # build parameters to pass to implicit function
-                local_u_approx = u_approx
+                local_u_approx = P.dtype_f(u_approx)
 
                 # note that derivatives of algebraic variables are taken into account here too
                 # these do not directly affect the output of eval_f but rather indirectly via QI
@@ -131,22 +119,14 @@ def implSystem(params):
                 sys = P.eval_f(local_u_approx, params_mesh, L.time + L.dt * self.coll.nodes[m - 1])
                 return sys
 
-            # get U_k+1
-            # note: not using solve_system here because this solve step is the same for any problem
-            # See link for how different methods use the default tol parameter
-            # https://github.com/scipy/scipy/blob/8a6f1a0621542f059a532953661cd43b8167fce0/scipy/optimize/_root.py#L220
-            # options['xtol'] = P.params.newton_tol
-            # options['eps'] = 1e-16
-
-            u_new = P.solve_system(implSystem, L.f[m], L.time + L.dt * self.coll.nodes[m - 1])
-
             # update gradient (recall L.f is being used to store the gradient)
-            L.f[m][:] = u_new
+            L.f[m] = P.solve_system(implSystem, L.f[m], L.time + L.dt * self.coll.nodes[m - 1])
 
         # Update solution approximation
         integral = self.integrate()
         for m in range(M):
             L.u[m + 1] = u_0 + integral[m]
+
         # indicate presence of new values at this level
         L.status.updated = True
 
@@ -213,11 +193,7 @@ def compute_residual(self, stage=None):
             L.status.residual = 0.0 if L.status.residual is None else L.status.residual
             return None
 
-        # check if there are new values (e.g. from a sweep)
-        # assert L.status.updated
-
         # compute the residual for each node
-
         res_norm = []
         for m in range(self.coll.num_nodes):
             # use abs function from data type here
@@ -242,3 +218,18 @@ def compute_residual(self, stage=None):
         L.status.updated = False
 
         return None
+
+    def compute_end_point(self):
+        """
+        Compute u at the right point of the interval
+
+        The value uend computed here is a full evaluation of the Picard formulation unless do_full_update==False
+
+        Returns:
+            None
+        """
+
+        if not self.coll.right_is_node or self.params.do_coll_update:
+            raise NotImplementedError()
+
+        super().compute_end_point()
diff --git a/pySDC/tests/test_projects/test_DAE/test_DAEMesh.py b/pySDC/tests/test_projects/test_DAE/test_DAEMesh.py
new file mode 100644
index 0000000000..297d35c800
--- /dev/null
+++ b/pySDC/tests/test_projects/test_DAE/test_DAEMesh.py
@@ -0,0 +1,76 @@
+import pytest
+
+
+@pytest.mark.base
+@pytest.mark.parametrize('shape', [(6,), (4, 6), (4, 6, 8)])
+def testInitialization(shape):
+    """
+    Tests for a random init if initialization results in desired shape of mesh.
+    """
+
+    import numpy as np
+    from pySDC.projects.DAE.misc.DAEMesh import DAEMesh
+
+    init = (shape, None, np.dtype('float64'))
+    mesh = DAEMesh(init)
+
+    assert np.shape(mesh.diff) == shape, f'ERROR: Component diff does not have the desired length!'
+    assert np.shape(mesh.alg) == shape, f'ERROR: Component alg does not have the desired length!'
+
+    assert len(mesh.components) == len(mesh), 'ERROR: Mesh does not contain two component arrays!'
+
+
+@pytest.mark.base
+def testInitializationGivenMesh():
+    """
+    Tests if for a given mesh the initialization results in the same mesh.
+    """
+
+    import numpy as np
+    from pySDC.projects.DAE.misc.DAEMesh import DAEMesh
+
+    nvars_1d = 6
+    init = (nvars_1d, None, np.dtype('float64'))
+    mesh1 = DAEMesh(init)
+    mesh1.diff[:] = np.arange(6)
+    mesh1.alg[:] = np.arange(6, 12)
+
+    mesh2 = DAEMesh(mesh1)
+
+    assert np.allclose(mesh1.diff, mesh2.diff) and np.allclose(
+        mesh1.alg, mesh2.alg
+    ), 'ERROR: Components in initialized meshes do not match!'
+
+
+@pytest.mark.base
+@pytest.mark.parametrize('shape', [(6,), (4, 6), (4, 6, 8)])
+def testArrayUFuncOperator(shape):
+    """
+    Test if overloaded __array_ufunc__ operator of datatype does what it is supposed to do.
+    """
+
+    import numpy as np
+    from pySDC.projects.DAE.misc.DAEMesh import DAEMesh
+
+    init = (shape, None, np.dtype('float64'))
+    mesh = DAEMesh(init)
+    mesh2 = DAEMesh(mesh)
+
+    randomArr = np.random.random(shape)
+    mesh.diff[:] = randomArr
+    mesh2.diff[:] = 2 * randomArr
+
+    subMesh = mesh - mesh2
+    assert type(subMesh) == DAEMesh
+    assert np.allclose(subMesh.diff, randomArr - 2 * randomArr)
+    assert np.allclose(subMesh.alg, 0)
+
+    addMesh = mesh + mesh2
+    assert type(addMesh) == DAEMesh
+    assert np.allclose(addMesh.diff, randomArr + 2 * randomArr)
+    assert np.allclose(addMesh.alg, 0)
+
+    sinMesh = np.sin(mesh)
+    assert type(sinMesh) == DAEMesh
+    assert np.allclose(sinMesh.diff, np.sin(randomArr))
+    assert np.allclose(sinMesh.alg, 0)
diff --git a/pySDC/tests/test_projects/test_DAE/test_HookClass_DAE.py b/pySDC/tests/test_projects/test_DAE/test_HookClass_DAE.py
new file mode 100644
index 0000000000..65f583b3ff
--- /dev/null
+++ b/pySDC/tests/test_projects/test_DAE/test_HookClass_DAE.py
@@ -0,0 +1,73 @@
+import numpy as np
+import pytest
+
+
+@pytest.mark.base
+@pytest.mark.parametrize('M', [2, 3, 4])
+def testHookClassDiffAlgComps(M):
+    """
+    Test if the hook class returns the correct errors.
+    """
+
+    from pySDC.helpers.stats_helper import get_sorted
+    from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE
+    from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE
+    from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
+    from pySDC.projects.DAE.misc.HookClass_DAE import (
+        LogGlobalErrorPostStepDifferentialVariable,
+        LogGlobalErrorPostStepAlgebraicVariable,
+    )
+
+    dt = 1e-2
+    level_params = {
+        'restol': 1e-13,
+        'dt': dt,
+    }
+
+    problem_params = {
+        'newton_tol': 1e-6,
+    }
+
+    sweeper_params = {
+        'quad_type': 'RADAU-RIGHT',
+        'num_nodes': M,
+        'QI': 'IE',
+    }
+
+    step_params = {
+        'maxiter': 45,
+    }
+
+    controller_params = {
+        'logger_level': 30,
+        'hook_class': [LogGlobalErrorPostStepDifferentialVariable, LogGlobalErrorPostStepAlgebraicVariable],
+    }
+
+    description = {
+        'problem_class': DiscontinuousTestDAE,
+        'problem_params': problem_params,
+        'sweeper_class': fully_implicit_DAE,
+        'sweeper_params': sweeper_params,
+        'level_params': level_params,
+        'step_params': step_params,
+    }
+
+    controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
+
+    t0 = 1.0
+    Tend = t0 + dt
+
+    P = controller.MS[0].levels[0].prob
+    uinit = P.u_exact(t0)
+    uex = P.u_exact(Tend)
+
+    uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
+
+    errHookDiff = np.array(get_sorted(stats, type='e_global_differential_post_step', sortby='time'))[:, 1]
+    errHookAlg = np.array(get_sorted(stats, type='e_global_algebraic_post_step', sortby='time'))[:, 1]
+
+    errRunDiff = abs(uex.diff[0] - uend.diff[0])
+    errRunAlg = abs(uex.alg[0] - uend.alg[0])
+
+    assert np.isclose(errHookDiff, errRunDiff), 'ERROR: Error in differential component does not match!'
+    assert np.isclose(errHookAlg, errRunAlg), 'ERROR: Error in algebraic component does not match!'
diff --git a/pySDC/tests/test_projects/test_DAE/test_misc.py b/pySDC/tests/test_projects/test_DAE/test_misc.py
deleted file mode 100644
index 0ab713951e..0000000000
--- a/pySDC/tests/test_projects/test_DAE/test_misc.py
+++ /dev/null
@@ -1,20 +0,0 @@
-import pytest
-
-
-#
-# Tests that problem class enforces parameter requirements
-@pytest.mark.base
-def test_problem_class_main():
-    from pySDC.projects.DAE.problems.simple_DAE import simple_dae_1
-
-    # initialize problem parameters
-    problem_params = dict()
-
-    # instantiate problem
-    try:
-        simple_dae_1(**problem_params)
-    # ensure error thrown is correct
-    except Exception as error:
-        assert type(error) == TypeError, "Parameter error was not thrown correctly"
-    else:
-        raise Exception("Parameter error was not thrown correctly")
diff --git a/pySDC/tests/test_projects/test_DAE/test_problems.py b/pySDC/tests/test_projects/test_DAE/test_problems.py
index f291a70fdc..ce03041562 100644
--- a/pySDC/tests/test_projects/test_DAE/test_problems.py
+++ b/pySDC/tests/test_projects/test_DAE/test_problems.py
@@ -10,15 +10,12 @@ def test_pendulum_u_exact_main():
     # initialize problem parameters
     problem_params = dict()
     problem_params['newton_tol'] = 1e-3  # tollerance for implicit solver
-    problem_params['nvars'] = 5
 
     # instantiate problem
     prob = pendulum_2d(**problem_params)
 
     u_test = prob.u_exact(5.0)
-    assert np.array_equal(u_test, np.zeros(5))
-
-    u_test = prob.u_exact(5.0)
+    assert np.isclose(abs(u_test), 0.0)
 
 
 @pytest.mark.base
@@ -28,15 +25,12 @@ def test_one_transistor_amplifier_u_exact_main():
     # initialize problem parameters
     problem_params = dict()
     problem_params['newton_tol'] = 1e-12  # tollerance for implicit solver
-    problem_params['nvars'] = 5
 
     # instantiate problem
     prob = one_transistor_amplifier(**problem_params)
 
     u_test = prob.u_exact(5.0)
-    assert np.array_equal(u_test, np.zeros(5))
-
-    u_test = prob.u_exact(5.0)
+    assert np.array_equal(abs(u_test), 0.0)
 
 
 @pytest.mark.base
@@ -46,15 +40,12 @@ def test_two_transistor_amplifier_u_exact_main():
     # initialize problem parameters
     problem_params = dict()
     problem_params['newton_tol'] = 1e-3  # tollerance for implicit solver
-    problem_params['nvars'] = 8
 
     # instantiate problem
     prob = two_transistor_amplifier(**problem_params)
 
     u_test = prob.u_exact(5.0)
-    assert np.array_equal(u_test, np.zeros(8))
-
-    u_test = prob.u_exact(5.0)
+    assert np.isclose(abs(u_test), 0.0)
 
 
 #
@@ -65,7 +56,6 @@ def test_pendulum_main():
     from pySDC.projects.DAE.problems.simple_DAE import pendulum_2d
     from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
     from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE
-    from pySDC.projects.DAE.misc.HookClass_DAE import error_hook
 
     # initialize level parameters
     level_params = dict()
@@ -80,7 +70,6 @@ def test_pendulum_main():
     # initialize problem parameters
     problem_params = dict()
     problem_params['newton_tol'] = 1e-3  # tollerance for implicit solver
-    problem_params['nvars'] = 5
 
     # initialize step parameters
     step_params = dict()
@@ -89,7 +78,6 @@ def test_pendulum_main():
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 30
-    # controller_params['hook_class'] = error_hook
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -113,10 +101,11 @@ def test_pendulum_main():
 
     # call main function to get things done...
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
-    uend_ref = [0.98613917, -0.16592027, 0.29956023, 1.77825875, 4.82500525]
-
+    uend_ref = P.dtype_u(P.init)
+    uend_ref.diff[:4] = (0.98613917, -0.16592027, 0.29956023, 1.77825875)
+    uend_ref.alg[0] = 4.82500525
     # check error
-    err = np.linalg.norm(uend - uend_ref, np.inf)
+    err = abs(uend.diff - uend_ref.diff)
     assert np.isclose(err, 0.0, atol=1e-4), "Error too large."
 
 
@@ -125,7 +114,6 @@ def test_one_transistor_amplifier_main():
     from pySDC.projects.DAE.problems.transistor_amplifier import one_transistor_amplifier
     from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
     from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE
-    from pySDC.projects.DAE.misc.HookClass_DAE import error_hook
 
     # initialize level parameters
     level_params = dict()
@@ -140,7 +128,6 @@ def test_one_transistor_amplifier_main():
     # initialize problem parameters
     problem_params = dict()
     problem_params['newton_tol'] = 1e-3  # tollerance for implicit solver
-    problem_params['nvars'] = 5
 
     # initialize step parameters
     step_params = dict()
@@ -149,7 +136,6 @@ def test_one_transistor_amplifier_main():
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 30
-    # controller_params['hook_class'] = error_hook
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -174,10 +160,11 @@ def test_one_transistor_amplifier_main():
     # call main function to get things done...
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
-    uend_ref = [-0.02182035, 3.06674603, 2.89634691, 2.45212382, -2.69727238]
+    uend_ref = P.dtype_u(P.init)
+    uend_ref[:] = (-0.02182035, 3.06674603, 2.89634691, 2.45212382, -2.69727238)
 
     # check error
-    err = np.linalg.norm(uend - uend_ref, np.inf)
+    err = abs(uend - uend_ref)
     assert np.isclose(err, 0.0, atol=1e-4), "Error too large."
 
 
@@ -186,7 +173,6 @@ def test_two_transistor_amplifier_main():
     from pySDC.projects.DAE.problems.transistor_amplifier import two_transistor_amplifier
     from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
     from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE
-    from pySDC.projects.DAE.misc.HookClass_DAE import error_hook
 
     # initialize level parameters
     level_params = dict()
@@ -201,7 +187,6 @@ def test_two_transistor_amplifier_main():
     # initialize problem parameters
     problem_params = dict()
     problem_params['newton_tol'] = 1e-3  # tollerance for implicit solver
-    problem_params['nvars'] = 8
 
     # initialize step parameters
     step_params = dict()
@@ -210,7 +195,6 @@ def test_two_transistor_amplifier_main():
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 30
-    # controller_params['hook_class'] = error_hook
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -235,7 +219,8 @@ def test_two_transistor_amplifier_main():
     # call main function to get things done...
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
-    uend_ref = [
+    uend_ref = P.dtype_u(P.init)
+    uend_ref[:] = (
         -5.52721527e-03,
         3.00630407e00,
         2.84974338e00,
@@ -244,10 +229,10 @@ def test_two_transistor_amplifier_main():
         2.19430889e00,
         5.89240699e00,
         9.99531182e-02,
-    ]
+    )
 
     # check error
-    err = np.linalg.norm(uend - uend_ref, np.inf)
+    err = abs(uend - uend_ref)
     assert np.isclose(err, 0.0, atol=1e-4), "Error too large."
 
 
@@ -256,7 +241,6 @@ def test_synchgen_infinite_bus_main():
     from pySDC.projects.DAE.problems.synchronous_machine import synchronous_machine_infinite_bus
     from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
     from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE
-    from pySDC.projects.DAE.misc.HookClass_DAE import error_hook
 
     # initialize level parameters
     level_params = dict()
@@ -271,7 +255,6 @@ def test_synchgen_infinite_bus_main():
     # initialize problem parameters
     problem_params = dict()
     problem_params['newton_tol'] = 1e-3  # tollerance for implicit solver
-    problem_params['nvars'] = 14
 
     # initialize step parameters
     step_params = dict()
@@ -280,7 +263,6 @@ def test_synchgen_infinite_bus_main():
     # initialize controller parameters
     controller_params = dict()
     controller_params['logger_level'] = 30
-    # controller_params['hook_class'] = error_hook
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -305,25 +287,29 @@ def test_synchgen_infinite_bus_main():
     # call main function to get things done...
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
-    uend_ref = [
+    uend_ref = P.dtype_u(P.init)
+    uend_ref.diff[:8] = (
         8.30823565e-01,
         -4.02584174e-01,
         1.16966755e00,
         9.47592808e-01,
         -3.68076863e-01,
         -3.87492326e-01,
+        3.10281509e-01,
+        9.94039645e-01,
+    )
+
+    uend_ref.alg[:6] = (
         -7.77837831e-01,
         -1.67347611e-01,
         1.34810867e00,
         5.46223705e-04,
         1.29690691e-02,
         -8.00823474e-02,
-        3.10281509e-01,
-        9.94039645e-01,
-    ]
+    )
 
     # check error
-    err = np.linalg.norm(uend - uend_ref, np.inf)
+    err = abs(uend.diff - uend_ref.diff)
     assert np.isclose(err, 0.0, atol=1e-4), "Error too large."
 
 
@@ -342,31 +328,37 @@ def test_DiscontinuousTestDAE_singularity():
     eps = 1e-3
     t_before_event = t_event - eps
     u_before_event = disc_test_DAE.u_exact(t_before_event)
-    du_before_event = (np.sinh(t_before_event), np.cosh(t_before_event))
+    du_before_event = disc_test_DAE.dtype_f(disc_test_DAE.init)
+    du_before_event.diff[0] = np.sinh(t_before_event)
+    du_before_event.alg[0] = np.cosh(t_before_event)
     f_before_event = disc_test_DAE.eval_f(u_before_event, du_before_event, t_before_event)
 
-    assert np.isclose(f_before_event[0], 0.0) and np.isclose(
-        f_before_event[1], 0.0
+    assert np.isclose(
+        abs(f_before_event), 0.0
     ), f"ERROR: Right-hand side after event does not match! Expected {(0.0, 0.0)}, got {f_before_event}"
 
     # test for t <= t^*
     u_event = disc_test_DAE.u_exact(t_event)
-    du_event = (np.sinh(t_event), np.cosh(t_event))
+    du_event = disc_test_DAE.dtype_f(disc_test_DAE.init)
+    du_event.diff[0] = np.sinh(t_event)
+    du_event.alg[0] = np.cosh(t_event)
     f_event = disc_test_DAE.eval_f(u_event, du_event, t_event)
 
-    assert np.isclose(f_event[0], 7 * np.sqrt(51.0)) and np.isclose(
-        f_event[1], 0.0
-    ), f"ERROR: Right-hand side at event does not match! Expected {(7 * np.sqrt(51), 0.0)}, got {f_event}"
+    assert np.isclose(f_event.diff[0], 7 * np.sqrt(51.0)) and np.isclose(
+        f_event.alg[0], 0.0
+    ), f"ERROR: Right-hand side at event does not match! Expected {(7 * np.sqrt(51), 0.0)}, got {(f_event.diff[0], f_event.alg[0])}"
 
     # test for t > t^* by setting t^* = t^* + eps
     t_after_event = t_event + eps
     u_after_event = disc_test_DAE.u_exact(t_after_event)
-    du_after_event = (np.sinh(t_event), np.cosh(t_event))
+    du_after_event = disc_test_DAE.dtype_f(disc_test_DAE.init)
+    du_after_event.diff[0] = np.sinh(t_event)
+    du_after_event.alg[0] = np.cosh(t_event)
     f_after_event = disc_test_DAE.eval_f(u_after_event, du_after_event, t_after_event)
 
-    assert np.isclose(f_after_event[0], 7 * np.sqrt(51.0)) and np.isclose(
-        f_after_event[1], 0.0
-    ), f"ERROR: Right-hand side after event does not match! Expected {(7 * np.sqrt(51), 0.0)}, got {f_after_event}"
+    assert np.isclose(f_after_event.diff[0], 7 * np.sqrt(51.0)) and np.isclose(
+        f_after_event.alg[0], 0.0
+    ), f"ERROR: Right-hand side after event does not match! Expected {(7 * np.sqrt(51), 0.0)}, got {(f_after_event.diff[0], f_after_event.alg[0])}"
 
 
 @pytest.mark.base
@@ -431,7 +423,7 @@ def test_DiscontinuousTestDAE_SDC(M):
 
     uend, _ = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
-    err = abs(uex[0] - uend[0])
+    err = abs(uex.diff[0] - uend.diff[0])
     assert err < err_tol[M], f"ERROR: Error is too large! Expected {err_tol[M]}, got {err}"
 
 
@@ -458,8 +450,8 @@ def test_DiscontinuousTestDAE_SDC_detection(M):
     }
 
     event_err_tol = {
-        2: 3.6968e-5,
-        3: 1.3496e-8,
+        2: 0.0011,
+        3: 0.01,
         4: 0.02,
         5: 0.0101,
     }
@@ -523,7 +515,7 @@ def test_DiscontinuousTestDAE_SDC_detection(M):
 
     uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
-    err = abs(uex[0] - uend[0])
+    err = abs(uex.diff[0] - uend.diff[0])
     assert err < err_tol[M], f"ERROR for M={M}: Error is too large! Expected {err_tol[M]}, got {err}"
 
     switches = get_sorted(stats, type='switch', sortby='time', recomputed=False)
@@ -557,21 +549,11 @@ def test_WSCC9_evaluation():
     # test if right-hand side of does have the correct length
     t0 = 0.0
     u0 = WSCC9.u_exact(t0)
-    du0 = np.zeros(len(u0))
+    du0 = WSCC9.dtype_f(WSCC9.init, val=0.0)
 
     f = WSCC9.eval_f(u0, du0, t0)
 
-    assert len(f) == nvars, 'Shape of f does not match with shape it is supposed to be!'
-
-    # test if ParameterError is raised if m != 3 or n != 9 is set
-    problem_params.update(
-        {
-            'm': 4,
-            'n': 8,
-        }
-    )
-    with pytest.raises(ParameterError):
-        WSCC9_test = WSCC9BusSystem(**problem_params)
+    assert len(f.diff) == nvars and len(f.alg) == nvars, 'Shape of f does not match with shape it is supposed to be!'
 
 
 @pytest.mark.base
@@ -630,7 +612,7 @@ def test_WSCC9_update_YBus():
 
     assert np.allclose(YBus_initial, YBus_initial_ref), 'YBus does not match with the YBus at initialization!'
 
-    uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
+    _, _ = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
     YBus_line_outage = P.YBus
     YBus_line6_8_outage = get_event_Ybus()
@@ -664,7 +646,7 @@ def test_WSCC9_SDC_detection():
 
     sweeper_params = {
         'quad_type': 'RADAU-RIGHT',
-        'num_nodes': 2,
+        'num_nodes': 3,
         'QI': 'LU',
     }
 
@@ -709,12 +691,12 @@ def test_WSCC9_SDC_detection():
     P = controller.MS[0].levels[0].prob
     uinit = P.u_exact(t0)
 
-    uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
+    _, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
 
     switches = get_sorted(stats, type='switch', sortby='time', recomputed=False)
     assert len(switches) >= 1, 'ERROR: No events found!'
     t_switch = [me[1] for me in switches][0]
-    assert np.isclose(t_switch, 0.6103290792685618, atol=1e-3), 'Found event does not match a threshold!'
+    assert np.isclose(t_switch, 0.5335289411812812, atol=1e-3), 'Found event does not match a threshold!'
 
 
 # @pytest.mark.base
diff --git a/pySDC/tests/test_projects/test_DAE/test_sweeper.py b/pySDC/tests/test_projects/test_DAE/test_sweeper.py
index e2dbd43b21..ff1c0bfb27 100644
--- a/pySDC/tests/test_projects/test_DAE/test_sweeper.py
+++ b/pySDC/tests/test_projects/test_DAE/test_sweeper.py
@@ -21,7 +21,6 @@ def test_predict_main():
     # initialize problem parameters
     problem_params = dict()
     problem_params['newton_tol'] = 1e-12  # tollerance for implicit solver
-    problem_params['nvars'] = 3
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -41,10 +40,10 @@ def test_predict_main():
     # call prediction function to initialise nodes
     L.sweep.predict()
     # check correct initialisation
-    assert np.array_equal(L.f[0], np.zeros(3))
+    assert np.allclose(abs(L.f[0]), 0.0)
     for i in range(sweeper_params['num_nodes']):
-        assert np.array_equal(L.u[i + 1], np.zeros(3))
-        assert np.array_equal(L.f[i + 1], np.zeros(3))
+        assert np.allclose(abs(L.u[i + 1]), 0.0)
+        assert np.allclose(abs(L.f[i + 1]), 0.0)
 
     # rerun check for random initialisation
     # expecting that random initialisation does not initialise to zero
@@ -58,10 +57,10 @@ def test_predict_main():
     # compute initial value (using the exact function here)
     L.u[0] = P.u_exact(L.time)
     L.sweep.predict()
-    assert np.array_equal(L.f[0], np.zeros(3))
+    assert abs(L.f[0]) == 0.0
     for i in range(sweeper_params['num_nodes']):
-        assert np.not_equal(L.u[i + 1], np.zeros(3)).any()
-        assert np.not_equal(L.f[i + 1], np.zeros(3)).any()
+        assert abs(L.u[i + 1]) > 0.0
+        assert abs(L.f[i + 1]) > 0.0
 
 
 @pytest.mark.base
@@ -83,7 +82,6 @@ def test_residual_main():
     # initialize problem parameters
     problem_params = dict()
     problem_params['newton_tol'] = 1e-12  # tollerance for implicit solver
-    problem_params['nvars'] = 3
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -101,8 +99,10 @@ def test_residual_main():
     # set reference values
     u = P.dtype_u(P.init)
     du = P.dtype_u(P.init)
-    u[:] = (5, 5, 5)
-    du[:] = (0, 0, 0)
+    u.diff[:2] = (5, 5)
+    u.alg[0] = 5
+    du.diff[:2] = (0, 0)
+    du.alg[0] = 0
     # set initial time in the status of the level
     L.status.time = 0.0
     L.u[0] = u
@@ -168,7 +168,6 @@ def test_compute_end_point_main():
     # initialize problem parameters
     problem_params = dict()
     problem_params['newton_tol'] = 1e-12  # tollerance for implicit solver
-    problem_params['nvars'] = 3
 
     # Fill description dictionary for easy hierarchy creation
     description = dict()
@@ -192,4 +191,4 @@ def test_compute_end_point_main():
     L.sweep.compute_end_point()
 
     for m in range(1, L.sweep.coll.num_nodes):
-        assert np.array_equal(L.u[m], L.uend), "ERROR: end point not computed correctly"
+        assert np.allclose(abs(L.u[m] - L.uend), 0.0), "ERROR: end point not computed correctly"