From 9b95a8100094d15bf2c2566192fe29e1c8ff20b8 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Wed, 8 Nov 2023 14:27:15 +0100
Subject: [PATCH 1/5] Inexactness and some other stuff

---
 .../estimate_extrapolation_error.py           |  17 ++-
 .../inexactness.py                            |  49 ++++++++-
 .../problem_classes/AllenCahn_2D_FD.py        |  77 +++++++++++--
 .../implementations/problem_classes/Quench.py |  12 +-
 pySDC/projects/Resilience/sweepers.py         |  44 +++++++-
 .../test_extrapolation_within_Q.py            | 104 ++++++++++++------
 6 files changed, 238 insertions(+), 65 deletions(-)

diff --git a/pySDC/implementations/convergence_controller_classes/estimate_extrapolation_error.py b/pySDC/implementations/convergence_controller_classes/estimate_extrapolation_error.py
index c7bba5fb5e..0052f3b420 100644
--- a/pySDC/implementations/convergence_controller_classes/estimate_extrapolation_error.py
+++ b/pySDC/implementations/convergence_controller_classes/estimate_extrapolation_error.py
@@ -459,6 +459,8 @@ def setup(self, controller, params, description, **kwargs):
         default_params = {
             'Taylor_order': 2 * num_nodes,
             'n': num_nodes,
+            'recompute_coefficients': False,
+            **params,
         }
 
         return {**super().setup(controller, params, description, **kwargs), **default_params}
@@ -485,18 +487,21 @@ def post_iteration_processing(self, controller, S, **kwargs):
         t_eval = S.time + nodes_[-1]
 
         dts = np.append(nodes_[0], nodes_[1:] - nodes_[:-1])
-        self.params.Taylor_order = 2 * len(nodes)
+        self.params.Taylor_order = len(nodes)
         self.params.n = len(nodes)
 
         # compute the extrapolation coefficients
-        # TODO: Maybe this can be reused
-        self.get_extrapolation_coefficients(nodes, dts, t_eval)
+        if None in self.coeff.u or self.params.recompute_coefficients:
+            self.get_extrapolation_coefficients(nodes, dts, t_eval)
 
         # compute the extrapolated solution
+        if lvl.f[0] is None:
+            lvl.f[0] = lvl.prob.eval_f(lvl.u[0], lvl.time)
+
         if type(lvl.f[0]) == imex_mesh:
-            f = [me.impl + me.expl for me in lvl.f]
+            f = [lvl.f[i].impl + lvl.f[i].expl if self.coeff.f[i] and lvl.f[i] else 0.0 for i in range(len(lvl.f) - 1)]
         elif type(lvl.f[0]) == mesh:
-            f = lvl.f
+            f = [lvl.f[i] if self.coeff.f[i] else 0.0 for i in range(len(lvl.f) - 1)]
         else:
             raise DataError(
                 f"Unable to store f from datatype {type(lvl.f[0])}, extrapolation based error estimate only\
@@ -506,7 +511,7 @@ def post_iteration_processing(self, controller, S, **kwargs):
         # compute the error with the weighted sum
         if self.comm:
             idx = (self.comm.rank + 1) % self.comm.size
-            sendbuf = self.coeff.u[idx] * lvl.u[idx] + self.coeff.f[idx] * lvl.f[idx]
+            sendbuf = self.coeff.u[idx] * lvl.u[idx] + self.coeff.f[idx] * f[idx]
             u_ex = lvl.prob.dtype_u(lvl.prob.init, val=0.0) if self.comm.rank == self.comm.size - 1 else None
             self.comm.Reduce(sendbuf, u_ex, op=self.sum, root=self.comm.size - 1)
         else:
diff --git a/pySDC/implementations/convergence_controller_classes/inexactness.py b/pySDC/implementations/convergence_controller_classes/inexactness.py
index 1a6ec9719c..3af8a83a96 100644
--- a/pySDC/implementations/convergence_controller_classes/inexactness.py
+++ b/pySDC/implementations/convergence_controller_classes/inexactness.py
@@ -24,8 +24,39 @@ def setup(self, controller, params, description, **kwargs):
             "ratio": 1e-2,
             "min_tol": 0,
             "max_tol": 1e99,
+            "maxiter": None,
+            "use_e_tol": 'e_tol' in description['level_params'].keys(),
+            "initial_tol": 1e-3,
+            **super().setup(controller, params, description, **kwargs),
         }
-        return {**defaults, **super().setup(controller, params, description, **kwargs)}
+        if defaults['maxiter']:
+            self.set_maxiter(description, defaults['maxiter'])
+        return defaults
+
+    def dependencies(self, controller, description, **kwargs):
+        """
+        Load the embedded error estimator if needed.
+
+        Args:
+            controller (pySDC.Controller): The controller
+            description (dict): The description object used to instantiate the controller
+
+        Returns:
+            None
+        """
+        super().dependencies(controller, description)
+
+        if self.params.use_e_tol:
+            from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import (
+                EstimateEmbeddedError,
+            )
+
+            controller.add_convergence_controller(
+                EstimateEmbeddedError,
+                description=description,
+            )
+
+        return None
 
     def post_iteration_processing(self, controller, step, **kwargs):
         """
@@ -39,8 +70,18 @@ def post_iteration_processing(self, controller, step, **kwargs):
             None
         """
         for lvl in step.levels:
-            lvl.prob.newton_tol = max(
-                [min([lvl.status.residual * self.params.ratio, self.params.max_tol]), self.params.min_tol]
+            SDC_accuracy = (
+                lvl.status.get('error_embedded_estimate', lvl.status.residual)
+                if self.params.use_e_tol
+                else lvl.status.residual
             )
+            SDC_accuracy = self.params.initial_tol if SDC_accuracy is None else SDC_accuracy
+            tol = max([min([SDC_accuracy * self.params.ratio, self.params.max_tol]), self.params.min_tol])
+            self.set_tolerance(lvl, tol)
+            self.log(f'Changed tolerance to {tol:.2e}', step)
+
+    def set_tolerance(self, lvl, tol):
+        lvl.prob.newton_tol = tol
 
-            self.log(f'Changed Newton tolerance to {lvl.prob.newton_tol:.2e}', step)
+    def set_maxiter(self, description, maxiter):
+        description['problem_params']['newton_maxiter'] = maxiter
diff --git a/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py b/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py
index a391dd5bec..04260f30ec 100644
--- a/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py
+++ b/pySDC/implementations/problem_classes/AllenCahn_2D_FD.py
@@ -3,7 +3,7 @@
 from scipy.sparse.linalg import cg
 
 from pySDC.core.Errors import ParameterError, ProblemError
-from pySDC.core.Problem import ptype
+from pySDC.core.Problem import ptype, WorkCounter
 from pySDC.helpers import problem_helper
 from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh, comp2_mesh
 
@@ -78,6 +78,7 @@ def __init__(
         newton_tol=1e-12,
         lin_tol=1e-8,
         lin_maxiter=100,
+        inexact_linear_ratio=None,
         radius=0.25,
         order=2,
     ):
@@ -96,14 +97,19 @@ def __init__(
             'nvars',
             'nu',
             'eps',
+            'radius',
+            'order',
+            localVars=locals(),
+            readOnly=True,
+        )
+        self._makeAttributeAndRegister(
             'newton_maxiter',
             'newton_tol',
             'lin_tol',
             'lin_maxiter',
-            'radius',
-            'order',
+            'inexact_linear_ratio',
             localVars=locals(),
-            readOnly=True,
+            readOnly=False,
         )
 
         # compute dx and get discretization matrix A
@@ -124,6 +130,10 @@ def __init__(
         self.newton_ncalls = 0
         self.lin_ncalls = 0
 
+        self.work_counters['newton'] = WorkCounter()
+        self.work_counters['rhs'] = WorkCounter()
+        self.work_counters['linear'] = WorkCounter()
+
     @staticmethod
     def __get_A(N, dx):
         """
@@ -198,6 +208,10 @@ def solve_system(self, rhs, factor, u0, t):
             # if g is close to 0, then we are done
             res = np.linalg.norm(g, np.inf)
 
+            # do inexactness in the linear solver
+            if self.inexact_linear_ratio:
+                self.lin_tol = res * self.inexact_linear_ratio
+
             if res < self.newton_tol:
                 break
 
@@ -206,11 +220,15 @@ def solve_system(self, rhs, factor, u0, t):
 
             # newton update: u1 = u0 - g/dg
             # u -= spsolve(dg, g)
-            u -= cg(dg, g, x0=z, tol=self.lin_tol, atol=0)[0]
+            u -= cg(
+                dg, g, x0=z, tol=self.lin_tol, maxiter=self.lin_maxiter, atol=0, callback=self.work_counters['linear']
+            )[0]
             # increase iteration count
             n += 1
             # print(n, res)
 
+            self.work_counters['newton']()
+
         # if n == self.newton_maxiter:
         #     raise ProblemError('Newton did not converge after %i iterations, error is %s' % (n, res))
 
@@ -242,9 +260,10 @@ def eval_f(self, u, t):
         v = u.flatten()
         f[:] = (self.A.dot(v) + 1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
 
+        self.work_counters['rhs']()
         return f
 
-    def u_exact(self, t):
+    def u_exact(self, t, u_init=None, t_init=None):
         r"""
         Routine to compute the exact solution at time :math:`t`.
 
@@ -258,13 +277,19 @@ def u_exact(self, t):
         me : dtype_u
             The exact solution.
         """
-
-        assert t == 0, 'ERROR: u_exact only valid for t=0'
         me = self.dtype_u(self.init, val=0.0)
-        for i in range(self.nvars[0]):
-            for j in range(self.nvars[1]):
-                r2 = self.xvalues[i] ** 2 + self.xvalues[j] ** 2
-                me[i, j] = np.tanh((self.radius - np.sqrt(r2)) / (np.sqrt(2) * self.eps))
+        if t > 0:
+
+            def eval_rhs(t, u):
+                return self.eval_f(u.reshape(self.init[0]), t).flatten()
+
+            me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
+
+        else:
+            for i in range(self.nvars[0]):
+                for j in range(self.nvars[1]):
+                    r2 = self.xvalues[i] ** 2 + self.xvalues[j] ** 2
+                    me[i, j] = np.tanh((self.radius - np.sqrt(r2)) / (np.sqrt(2) * self.eps))
 
         return me
 
@@ -310,6 +335,7 @@ def eval_f(self, u, t):
         f.impl[:] = self.A.dot(v).reshape(self.nvars)
         f.expl[:] = (1.0 / self.eps**2 * v * (1.0 - v**self.nu)).reshape(self.nvars)
 
+        self.work_counters['rhs']()
         return f
 
     def solve_system(self, rhs, factor, u0, t):
@@ -338,6 +364,7 @@ class context:
 
         def callback(xk):
             context.num_iter += 1
+            self.work_counters['linear']()
             return context.num_iter
 
         me = self.dtype_u(self.init)
@@ -359,6 +386,32 @@ def callback(xk):
 
         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.
+
+        Returns
+        -------
+        me : dtype_u
+            The exact solution.
+        """
+        me = self.dtype_u(self.init, val=0.0)
+        if t > 0:
+
+            def eval_rhs(t, u):
+                f = self.eval_f(u.reshape(self.init[0]), t)
+                return (f.impl + f.expl).flatten()
+
+            me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init)
+        else:
+            me[:] = super().u_exact(t, u_init, t_init)
+        return me
+
 
 # noinspection PyUnusedLocal
 class allencahn_semiimplicit_v2(allencahn_fullyimplicit):
diff --git a/pySDC/implementations/problem_classes/Quench.py b/pySDC/implementations/problem_classes/Quench.py
index 8527d29d8a..2175d197b6 100644
--- a/pySDC/implementations/problem_classes/Quench.py
+++ b/pySDC/implementations/problem_classes/Quench.py
@@ -60,6 +60,10 @@ class Quench(ptype):
         Maximum number of linear iterations inside the Newton solver.
     direct_solver : bool, optional
         Indicates if a direct solver should be used.
+    inexact_linear_ratio : float, optional
+        Ratio of tolerance of linear solver to the Newton residual, overrides `lintol`
+    min_lintol : float, optional
+        Minimal tolerance for the linear solver
     reference_sol_type : str, optional
         Indicates which method should be used to compute a reference solution.
         Choose between ``'scipy'``, ``'SDC'``, or ``'DIRK'``.
@@ -106,6 +110,7 @@ def __init__(
         liniter=99,
         direct_solver=True,
         inexact_linear_ratio=None,
+        min_lintol=1e-12,
         reference_sol_type='scipy',
     ):
         """
@@ -142,6 +147,7 @@ def __init__(
             'lintol',
             'liniter',
             'inexact_linear_ratio',
+            'min_lintol',
             localVars=locals(),
             readOnly=False,
         )
@@ -307,11 +313,11 @@ def solve_system(self, rhs, factor, u0, t):
         u = self.dtype_u(u0)
         res = np.inf
         delta = self.dtype_u(self.init, val=0.0)
-        z = self.dtype_u(self.init, val=0.0)
 
         # construct a preconditioner for the space solver
         if not self.direct_solver:
             M = inv((self.Id - factor * self.A).toarray())
+            zero = self.dtype_u(self.init, val=0.0)
 
         for n in range(0, self.newton_maxiter):
             # assemble G such that G(u) = 0 at the solution of the step
@@ -325,7 +331,7 @@ def solve_system(self, rhs, factor, u0, t):
                 break
 
             if self.inexact_linear_ratio:
-                self.lintol = max([res * self.inexact_linear_ratio, 1e-12])
+                self.lintol = max([res * self.inexact_linear_ratio, self.min_lintol])
 
             # assemble Jacobian J of G
             J = self.Id - factor * (self.A + self.get_non_linear_Jacobian(u))
@@ -337,7 +343,7 @@ def solve_system(self, rhs, factor, u0, t):
                 delta, info = gmres(
                     J,
                     G,
-                    x0=z,
+                    x0=zero,
                     M=M,
                     tol=self.lintol,
                     maxiter=self.liniter,
diff --git a/pySDC/projects/Resilience/sweepers.py b/pySDC/projects/Resilience/sweepers.py
index e144264ba3..f14ba0120d 100644
--- a/pySDC/projects/Resilience/sweepers.py
+++ b/pySDC/projects/Resilience/sweepers.py
@@ -127,10 +127,25 @@ def update_nodes(self):
         return None
 
     def compute_residual(self, stage=''):
-        if stage not in self.params.skip_residual_computation:
-            lvl = self.level
-            lvl.f[-1] = lvl.prob.eval_f(lvl.u[-1], lvl.time + lvl.dt * self.coll.nodes[-1])
-        super().compute_residual(stage)
+        lvl = self.level
+
+        if lvl.params.residual_type[:4] == 'full' or stage in self.params.skip_residual_computation:
+            super().compute_residual(stage=stage)
+            return None
+
+        # compute residual of last rank without involvement of all of them
+        lvl.f[-1] = lvl.prob.eval_f(lvl.u[-1], lvl.time + lvl.dt * self.coll.nodes[-1])
+
+        res = lvl.u[0] - lvl.u[-1]
+        for m in range(1, self.coll.num_nodes + 1):
+            res += lvl.dt * self.coll.Qmat[-1, m] * lvl.f[m]
+
+        if lvl.params.residual_type[-3:] == 'abs':
+            lvl.status.residual = abs(res)
+        else:
+            lvl.status.residual = abs(res) / abs(lvl.u[0])
+
+        lvl.status.updated = False
 
 
 class imex_1st_order_efficient(efficient_sweeper, imex_1st_order):
@@ -218,3 +233,24 @@ def update_nodes(self):
         L.status.updated = True
 
         return None
+
+    def compute_residual(self, stage=''):
+        lvl = self.level
+
+        if lvl.params.residual_type[:4] == 'full' or stage in self.params.skip_residual_computation:
+            super().compute_residual(stage=stage)
+            return None
+
+        # compute residual of last rank without involvement of all of them
+        lvl.f[-1] = lvl.prob.eval_f(lvl.u[-1], lvl.time + lvl.dt * self.coll.nodes[-1])
+
+        res = lvl.u[0] - lvl.u[-1]
+        for m in range(1, self.coll.num_nodes + 1):
+            res += lvl.dt * self.coll.Qmat[-1, m] * (lvl.f[m].impl + lvl.f[m].expl)
+
+        if lvl.params.residual_type[-3:] == 'abs':
+            lvl.status.residual = abs(res)
+        else:
+            lvl.status.residual = abs(res) / abs(lvl.u[0])
+
+        lvl.status.updated = False
diff --git a/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py b/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py
index 73a3297d3d..eae18798b5 100644
--- a/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py
+++ b/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py
@@ -1,25 +1,22 @@
 import pytest
 
 
-def single_run(dt, Tend, num_nodes, quad_type, QI, useMPI):
+def get_controller(dt, num_nodes, quad_type, useMPI, **kwargs):
     """
-    Runs a single advection problem with certain parameters
+    Gets a controller setup for the polynomial test problem.
 
     Args:
         dt (float): Step size
-        Tend (float): Final time
         num_nodes (int): Number of nodes
         quad_type (str): Type of quadrature
-        QI (str): Preconditioner
         useMPI (bool): Whether or not to use MPI
 
     Returns:
        (dict): Stats object generated during the run
        (pySDC.Controller.controller): Controller used in the run
     """
-    from pySDC.implementations.problem_classes.AdvectionEquation_ND_FD import advectionNd
+    from pySDC.implementations.problem_classes.polynomial_test_problem import polynomial_testequation
     from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
-    from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostStep
     from pySDC.implementations.convergence_controller_classes.estimate_extrapolation_error import (
         EstimateExtrapolationErrorWithinQ,
     )
@@ -37,30 +34,28 @@ def single_run(dt, Tend, num_nodes, quad_type, QI, useMPI):
     # initialize level parameters
     level_params = {}
     level_params['dt'] = dt
-    level_params['restol'] = 1e-10
+    level_params['restol'] = 1.0
 
     # initialize sweeper parameters
     sweeper_params = {}
     sweeper_params['quad_type'] = quad_type
     sweeper_params['num_nodes'] = num_nodes
-    sweeper_params['QI'] = QI
     sweeper_params['comm'] = comm
 
-    problem_params = {'freq': 2, 'nvars': 2**9, 'c': 1.0, 'stencil_type': 'center', 'order': 6, 'bc': 'periodic'}
+    problem_params = {'degree': 20}
 
     # initialize step parameters
     step_params = {}
-    step_params['maxiter'] = 99
+    step_params['maxiter'] = 0
 
     # initialize controller parameters
     controller_params = {}
     controller_params['logger_level'] = 30
-    controller_params['hook_class'] = LogGlobalErrorPostStep
     controller_params['mssdc_jac'] = False
 
     # fill description dictionary for easy step instantiation
     description = {}
-    description['problem_class'] = advectionNd
+    description['problem_class'] = polynomial_testequation
     description['problem_params'] = problem_params
     description['sweeper_class'] = sweeper_class
     description['sweeper_params'] = sweeper_params
@@ -68,19 +63,57 @@ def single_run(dt, Tend, num_nodes, quad_type, QI, useMPI):
     description['step_params'] = step_params
     description['convergence_controllers'] = {EstimateExtrapolationErrorWithinQ: {}}
 
-    # set time parameters
-    t0 = 0.0
-
-    # instantiate controller
     controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
+    return controller
 
-    # get initial values on finest level
-    P = controller.MS[0].levels[0].prob
-    uinit = P.u_exact(t0)
 
-    # call main function to get things done...
-    uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
-    return stats, controller
+def single_test(**kwargs):
+    """
+    Run a single test where the solution is replaced by a polynomial and the nodes are changed.
+    Because we know the polynomial going in, we can check if the interpolation based change was
+    exact. If the solution is not a polynomial or a polynomial of higher degree then the number
+    of nodes, the change in nodes does add some error, of course, but here it is on the order of
+    machine precision.
+    """
+    import numpy as np
+
+    args = {
+        'num_nodes': 3,
+        'quad_type': 'RADAU-RIGHT',
+        'useMPI': False,
+        'dt': 1.0,
+        **kwargs,
+    }
+
+    # prepare variables
+    controller = get_controller(**args)
+    step = controller.MS[0]
+    level = step.levels[0]
+    prob = level.prob
+    cont = controller.convergence_controllers[
+        np.arange(len(controller.convergence_controllers))[
+            [type(me).__name__ == 'EstimateExtrapolationErrorWithinQ' for me in controller.convergence_controllers]
+        ][0]
+    ]
+    nodes = np.append([0], level.sweep.coll.nodes)
+
+    # initialize variables
+    step.status.slot = 0
+    step.status.iter = 1
+    level.status.time = 0.0
+    level.status.residual = 0.0
+    level.u[0] = prob.u_exact(t=0)
+    level.sweep.predict()
+
+    for i in range(len(level.u)):
+        if level.u[i] is not None:
+            level.u[i][:] = prob.u_exact(nodes[i] * level.dt)
+
+    # perform the interpolation
+    cont.post_iteration_processing(controller, step)
+    error = level.status.error_extrapolation_estimate
+
+    return error
 
 
 def multiple_runs(dts, **kwargs):
@@ -101,14 +134,10 @@ def multiple_runs(dts, **kwargs):
     res = {}
 
     for dt in dts:
-        stats, controller = single_run(Tend=5.0 * dt, dt=dt, **kwargs)
-
         res[dt] = {}
-        res[dt]['e_loc'] = max([me[1] for me in get_sorted(stats, type='e_global_post_step')])
-        res[dt]['e_ex'] = max([me[1] for me in get_sorted(stats, type='error_extrapolation_estimate')])
+        res[dt]['e'] = single_test(dt=dt, **kwargs)
 
-    coll_order = controller.MS[0].levels[0].sweep.coll.order
-    return res, coll_order
+    return res
 
 
 def check_order(dts, **kwargs):
@@ -122,19 +151,18 @@ def check_order(dts, **kwargs):
     """
     import numpy as np
 
-    res, coll_order = multiple_runs(dts, **kwargs)
+    res = multiple_runs(dts, **kwargs)
     dts = np.array(list(res.keys()))
     keys = list(res[dts[0]].keys())
 
     expected_order = {
-        'e_loc': coll_order + 1,
-        'e_ex': kwargs['num_nodes'] + 1,
+        'e': kwargs['num_nodes'],
     }
 
     for key in keys:
         errors = np.array([res[dt][key] for dt in dts])
 
-        mask = np.logical_and(errors < 1e-0, errors > 1e-10)
+        mask = np.logical_and(errors < 1e-0, errors > 1e-14)
         order = np.log(errors[mask][1:] / errors[mask][:-1]) / np.log(dts[mask][1:] / dts[mask][:-1])
 
         assert np.isclose(
@@ -143,7 +171,7 @@ def check_order(dts, **kwargs):
 
 
 @pytest.mark.base
-@pytest.mark.parametrize('num_nodes', [2, 3])
+@pytest.mark.parametrize('num_nodes', [2, 3, 4])
 @pytest.mark.parametrize('quad_type', ['RADAU-RIGHT', 'GAUSS'])
 def test_extrapolation_within_Q(num_nodes, quad_type):
     kwargs = {
@@ -152,12 +180,16 @@ def test_extrapolation_within_Q(num_nodes, quad_type):
         'useMPI': False,
         'QI': 'MIN',
     }
-    check_order([5e-1, 1e-1, 8e-2, 5e-2], **kwargs)
+
+    import numpy as np
+
+    steps = np.logspace(1, -4, 20)
+    check_order(steps, **kwargs)
 
 
 @pytest.mark.mpi4py
-@pytest.mark.parametrize('num_nodes', [2, 3])
-@pytest.mark.parametrize('quad_type', ['RADAU-RIGHT'])
+@pytest.mark.parametrize('num_nodes', [2, 4])
+@pytest.mark.parametrize('quad_type', ['RADAU-RIGHT', 'GAUSS'])
 def test_extrapolation_within_Q_MPI(num_nodes, quad_type):
     import subprocess
     import os

From dba1d8f890f19c85ed51f807103dd44e1d5bd682 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Wed, 8 Nov 2023 14:58:14 +0100
Subject: [PATCH 2/5] Fixes

---
 pySDC/projects/Resilience/strategies.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py
index 3e63878bee..654d1492d5 100644
--- a/pySDC/projects/Resilience/strategies.py
+++ b/pySDC/projects/Resilience/strategies.py
@@ -1457,8 +1457,8 @@ def get_reference_value(self, problem, key, op, num_procs=1):
         """
         if problem.__name__ == "run_vdp":
             if key == 'work_newton' and op == sum:
-                return 2677
+                return 3443
             elif key == 'e_global_post_run' and op == max:
-                return 4.375184403937471e-06
+                return 4.929282266752377e-06
 
         raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!')

From 76c3ba381364d7a0fe4d1ca5aac8bf6d507d7542 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Wed, 8 Nov 2023 16:33:50 +0100
Subject: [PATCH 3/5] Fixes

---
 pySDC/projects/Resilience/sweepers.py         |  6 -----
 .../test_efficient_sweepers.py                |  5 ++--
 .../test_extrapolated_error_within_Q.py       | 23 -------------------
 3 files changed, 3 insertions(+), 31 deletions(-)
 delete mode 100644 pySDC/tests/test_projects/test_resilience/test_extrapolated_error_within_Q.py

diff --git a/pySDC/projects/Resilience/sweepers.py b/pySDC/projects/Resilience/sweepers.py
index f14ba0120d..f75b988bb0 100644
--- a/pySDC/projects/Resilience/sweepers.py
+++ b/pySDC/projects/Resilience/sweepers.py
@@ -133,9 +133,6 @@ def compute_residual(self, stage=''):
             super().compute_residual(stage=stage)
             return None
 
-        # compute residual of last rank without involvement of all of them
-        lvl.f[-1] = lvl.prob.eval_f(lvl.u[-1], lvl.time + lvl.dt * self.coll.nodes[-1])
-
         res = lvl.u[0] - lvl.u[-1]
         for m in range(1, self.coll.num_nodes + 1):
             res += lvl.dt * self.coll.Qmat[-1, m] * lvl.f[m]
@@ -241,9 +238,6 @@ def compute_residual(self, stage=''):
             super().compute_residual(stage=stage)
             return None
 
-        # compute residual of last rank without involvement of all of them
-        lvl.f[-1] = lvl.prob.eval_f(lvl.u[-1], lvl.time + lvl.dt * self.coll.nodes[-1])
-
         res = lvl.u[0] - lvl.u[-1]
         for m in range(1, self.coll.num_nodes + 1):
             res += lvl.dt * self.coll.Qmat[-1, m] * (lvl.f[m].impl + lvl.f[m].expl)
diff --git a/pySDC/tests/test_projects/test_resilience/test_efficient_sweepers.py b/pySDC/tests/test_projects/test_resilience/test_efficient_sweepers.py
index c03ac8cca3..3dab46f6d1 100644
--- a/pySDC/tests/test_projects/test_resilience/test_efficient_sweepers.py
+++ b/pySDC/tests/test_projects/test_resilience/test_efficient_sweepers.py
@@ -12,6 +12,7 @@ def run_Lorenz(efficient, skip_residual_computation, num_procs=1):
     # initialize level parameters
     level_params = {}
     level_params['dt'] = 1e-1
+    level_params['residual_type'] = 'last_rel'
 
     # initialize sweeper parameters
     sweeper_params = {}
@@ -77,6 +78,7 @@ def run_Schroedinger(efficient=False, num_procs=1, skip_residual_computation=Fal
     level_params['restol'] = 1e-8
     level_params['dt'] = 2e-01
     level_params['nsweeps'] = 1
+    level_params['residual_type'] = 'last_rel'
 
     # initialize sweeper parameters
     sweeper_params = {}
@@ -133,7 +135,7 @@ def run_Schroedinger(efficient=False, num_procs=1, skip_residual_computation=Fal
 
 
 @pytest.mark.base
-def test_generic_implicit_efficient(skip_residual_computation=True):
+def test_generic_implicit_efficient(skip_residual_computation=False):
     stats_normal = run_Lorenz(efficient=False, skip_residual_computation=skip_residual_computation)
     stats_efficient = run_Lorenz(efficient=True, skip_residual_computation=skip_residual_computation)
     assert_sameness(stats_normal, stats_efficient, 'generic_implicit')
@@ -145,7 +147,6 @@ def test_residual_skipping():
     stats_normal = run_Lorenz(efficient=True, skip_residual_computation=False)
     stats_efficient = run_Lorenz(efficient=True, skip_residual_computation=True)
     assert_sameness(stats_normal, stats_efficient, 'generic_implicit', check_residual=False)
-    assert_benefit(stats_normal, stats_efficient)
 
 
 @pytest.mark.mpi4py
diff --git a/pySDC/tests/test_projects/test_resilience/test_extrapolated_error_within_Q.py b/pySDC/tests/test_projects/test_resilience/test_extrapolated_error_within_Q.py
deleted file mode 100644
index f090379c60..0000000000
--- a/pySDC/tests/test_projects/test_resilience/test_extrapolated_error_within_Q.py
+++ /dev/null
@@ -1,23 +0,0 @@
-import pytest
-
-
-@pytest.mark.base
-@pytest.mark.parametrize("prob_name", ['advection', 'piline'])
-@pytest.mark.parametrize('num_nodes', [2, 3])
-@pytest.mark.parametrize('quad_type', ['RADAU-RIGHT', 'GAUSS'])
-def test_order_extrapolation_estimate_within_Q(prob_name, num_nodes, quad_type):
-    from pySDC.projects.Resilience.extrapolation_within_Q import check_order
-
-    if prob_name == 'advection':
-        from pySDC.projects.Resilience.advection import run_advection
-
-        prob = run_advection
-    elif prob_name == 'piline':
-        from pySDC.projects.Resilience.piline import run_piline
-
-        prob = run_piline
-
-    else:
-        raise NotImplementedError(f'Problem \"{prob_name}\" not implemented in this test!')
-
-    check_order(None, prob=prob, dts=[5e-1, 1e-1, 5e-2, 1e-2], num_nodes=num_nodes, quad_type=quad_type)

From 9976f09f9825fa1496fd8f39b15d6997795b8f4c Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Thu, 9 Nov 2023 07:22:49 +0100
Subject: [PATCH 4/5] Hopefully fixed a flaky test

---
 .../test_extrapolation_within_Q.py                            | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py b/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py
index eae18798b5..b2e16a6532 100644
--- a/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py
+++ b/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py
@@ -162,7 +162,7 @@ def check_order(dts, **kwargs):
     for key in keys:
         errors = np.array([res[dt][key] for dt in dts])
 
-        mask = np.logical_and(errors < 1e-0, errors > 1e-14)
+        mask = np.logical_and(errors < 1e-1, errors > 1e-12)
         order = np.log(errors[mask][1:] / errors[mask][:-1]) / np.log(dts[mask][1:] / dts[mask][:-1])
 
         assert np.isclose(
@@ -183,7 +183,7 @@ def test_extrapolation_within_Q(num_nodes, quad_type):
 
     import numpy as np
 
-    steps = np.logspace(1, -4, 20)
+    steps = np.logspace(-1, -3, 10)
     check_order(steps, **kwargs)
 
 

From 375bc5e01fa3a0163584c495c243badefbe835ce Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Thu, 9 Nov 2023 09:29:05 +0100
Subject: [PATCH 5/5] IMEX version of the polynomial test problem

---
 .../polynomial_test_problem.py                | 34 ++++++++++++++++++-
 .../test_extrapolation_within_Q.py            | 17 ++++++++--
 2 files changed, 47 insertions(+), 4 deletions(-)

diff --git a/pySDC/implementations/problem_classes/polynomial_test_problem.py b/pySDC/implementations/problem_classes/polynomial_test_problem.py
index 3878227b24..4d813630f4 100644
--- a/pySDC/implementations/problem_classes/polynomial_test_problem.py
+++ b/pySDC/implementations/problem_classes/polynomial_test_problem.py
@@ -1,7 +1,7 @@
 import numpy as np
 
 from pySDC.core.Problem import ptype
-from pySDC.implementations.datatype_classes.mesh import mesh
+from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh
 
 
 class polynomial_testequation(ptype):
@@ -88,3 +88,35 @@ def u_exact(self, t, **kwargs):
         me = self.dtype_u(self.init)
         me[:] = self.poly(t)
         return me
+
+
+class polynomial_testequation_IMEX(polynomial_testequation):
+    """
+    IMEX version of the polynomial test problem that assigns half the derivative to the implicit part and the other half to the explicit part.
+    Keep in mind that you still cannot Really perform any solves.
+    """
+
+    dtype_f = imex_mesh
+
+    def eval_f(self, u, t):
+        """
+        Derivative of the polynomial.
+
+        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.
+        """
+
+        f = self.dtype_f(self.init)
+        derivative = self.poly.deriv(m=1)(t)
+        f.impl[:] = derivative / 2
+        f.expl[:] = derivative / 2
+        return f
diff --git a/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py b/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py
index b2e16a6532..964a6e84b3 100644
--- a/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py
+++ b/pySDC/tests/test_convergence_controllers/test_extrapolation_within_Q.py
@@ -1,7 +1,7 @@
 import pytest
 
 
-def get_controller(dt, num_nodes, quad_type, useMPI, **kwargs):
+def get_controller(dt, num_nodes, quad_type, useMPI, imex, **kwargs):
     """
     Gets a controller setup for the polynomial test problem.
 
@@ -10,17 +10,26 @@ def get_controller(dt, num_nodes, quad_type, useMPI, **kwargs):
         num_nodes (int): Number of nodes
         quad_type (str): Type of quadrature
         useMPI (bool): Whether or not to use MPI
+        imex (bool): Use IMEX version of the test problem
 
     Returns:
        (dict): Stats object generated during the run
        (pySDC.Controller.controller): Controller used in the run
     """
-    from pySDC.implementations.problem_classes.polynomial_test_problem import polynomial_testequation
     from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
     from pySDC.implementations.convergence_controller_classes.estimate_extrapolation_error import (
         EstimateExtrapolationErrorWithinQ,
     )
 
+    if imex:
+        from pySDC.implementations.problem_classes.polynomial_test_problem import (
+            polynomial_testequation_IMEX as problem_class,
+        )
+    else:
+        from pySDC.implementations.problem_classes.polynomial_test_problem import (
+            polynomial_testequation as problem_class,
+        )
+
     if useMPI:
         from pySDC.implementations.sweeper_classes.generic_implicit_MPI import generic_implicit_MPI as sweeper_class
         from mpi4py import MPI
@@ -55,7 +64,7 @@ def get_controller(dt, num_nodes, quad_type, useMPI, **kwargs):
 
     # fill description dictionary for easy step instantiation
     description = {}
-    description['problem_class'] = polynomial_testequation
+    description['problem_class'] = problem_class
     description['problem_params'] = problem_params
     description['sweeper_class'] = sweeper_class
     description['sweeper_params'] = sweeper_params
@@ -179,6 +188,7 @@ def test_extrapolation_within_Q(num_nodes, quad_type):
         'quad_type': quad_type,
         'useMPI': False,
         'QI': 'MIN',
+        'imex': False,
     }
 
     import numpy as np
@@ -219,5 +229,6 @@ def test_extrapolation_within_Q_MPI(num_nodes, quad_type):
             'quad_type': sys.argv[2],
             'useMPI': True,
             'QI': 'MIN',
+            'imex': True,
         }
         check_order([5e-1, 1e-1, 8e-2, 5e-2], **kwargs)