From 4d10733d8ff0d7e097df7d4988dbca49af09762a Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Fri, 11 Oct 2024 14:40:37 +0200
Subject: [PATCH 01/14] Added basic infrastructure for Rayleigh-Benard in
 resilience project

---
 .../step_size_limiter.py                      |   2 +-
 pySDC/projects/Resilience/RBC.py              | 165 ++++++++++++++++++
 2 files changed, 166 insertions(+), 1 deletion(-)
 create mode 100644 pySDC/projects/Resilience/RBC.py

diff --git a/pySDC/implementations/convergence_controller_classes/step_size_limiter.py b/pySDC/implementations/convergence_controller_classes/step_size_limiter.py
index cb090eebd3..7551ad2cfa 100644
--- a/pySDC/implementations/convergence_controller_classes/step_size_limiter.py
+++ b/pySDC/implementations/convergence_controller_classes/step_size_limiter.py
@@ -146,7 +146,7 @@ def get_new_step_size(self, controller, S, **kwargs):
                         S,
                     )
                     L.status.dt_new = dt_new
-                elif abs(L.status.dt_new / L.params.dt - 1) < self.params.dt_rel_min_slope:
+                elif abs(L.status.dt_new / L.params.dt - 1) < self.params.dt_rel_min_slope and not S.status.restart:
                     L.status.dt_new = L.params.dt
                     self.log(
                         f"Step size did not change sufficiently to warrant step size change, keeping {L.status.dt_new:.2e}",
diff --git a/pySDC/projects/Resilience/RBC.py b/pySDC/projects/Resilience/RBC.py
new file mode 100644
index 0000000000..280283f557
--- /dev/null
+++ b/pySDC/projects/Resilience/RBC.py
@@ -0,0 +1,165 @@
+# script to run a Rayleigh-Benard Convection problem
+from pySDC.implementations.problem_classes.generic_spectral import compute_residual_DAE
+from pySDC.implementations.problem_classes.RayleighBenard import RayleighBenard
+from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
+from pySDC.projects.Resilience.hook import hook_collection, LogData
+from pySDC.projects.Resilience.strategies import merge_descriptions
+from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
+from pySDC.core.convergence_controller import ConvergenceController
+
+from pySDC.core.errors import ConvergenceError
+
+
+class ReachTendExactly(ConvergenceController):
+
+    def setup(self, controller, params, description, **kwargs):
+        defaults = {
+            "control_order": +50,
+            "Tend": None,
+        }
+        return {**defaults, **super().setup(controller, params, description, **kwargs)}
+
+    def get_new_step_size(self, controller, step, **kwargs):
+        L = step.levels[0]
+        L.status.dt_new = min([L.params.dt, self.params.Tend - L.time - L.dt])
+
+
+def run_RBC(
+    custom_description=None,
+    num_procs=1,
+    Tend=14.0,
+    hook_class=LogData,
+    fault_stuff=None,
+    custom_controller_params=None,
+    u0=None,
+    t0=10,
+    use_MPI=False,
+    **kwargs,
+):
+    """
+    Args:
+        custom_description (dict): Overwrite presets
+        num_procs (int): Number of steps for MSSDC
+        Tend (float): Time to integrate to
+        hook_class (pySDC.Hook): A hook to store data
+        fault_stuff (dict): A dictionary with information on how to add faults
+        custom_controller_params (dict): Overwrite presets
+        u0 (dtype_u): Initial value
+        t0 (float): Starting time
+        use_MPI (bool): Whether or not to use MPI
+
+    Returns:
+        dict: The stats object
+        controller: The controller
+        bool: If the code crashed
+    """
+    level_params = {}
+    level_params['dt'] = 5e-4
+    level_params['restol'] = -1
+
+    sweeper_params = {}
+    sweeper_params['quad_type'] = 'RADAU-RIGHT'
+    sweeper_params['num_nodes'] = 3
+    sweeper_params['QI'] = 'LU'
+    sweeper_params['QE'] = 'PIC'
+
+    problem_params = {}
+
+    step_params = {}
+    step_params['maxiter'] = 5
+
+    convergence_controllers = {}
+    convergence_controllers[ReachTendExactly] = {'Tend': Tend}
+
+    controller_params = {}
+    controller_params['logger_level'] = 30
+    controller_params['hook_class'] = hook_collection + (hook_class if type(hook_class) == list else [hook_class])
+    controller_params['mssdc_jac'] = False
+
+    if custom_controller_params is not None:
+        controller_params = {**controller_params, **custom_controller_params}
+
+    imex_1st_order.compute_residual = compute_residual_DAE
+
+    RayleighBenard._u_exact = RayleighBenard.u_exact
+    RayleighBenard.u_exact = u_exact
+
+    description = {}
+    description['problem_class'] = RayleighBenard
+    description['problem_params'] = problem_params
+    description['sweeper_class'] = imex_1st_order
+    description['sweeper_params'] = sweeper_params
+    description['level_params'] = level_params
+    description['step_params'] = step_params
+
+    if custom_description is not None:
+        description = merge_descriptions(description, custom_description)
+
+    controller_args = {
+        'controller_params': controller_params,
+        'description': description,
+    }
+    if use_MPI:
+        from mpi4py import MPI
+        from pySDC.implementations.controller_classes.controller_MPI import controller_MPI
+
+        comm = kwargs.get('comm', MPI.COMM_WORLD)
+        controller = controller_MPI(**controller_args, comm=comm)
+        P = controller.S.levels[0].prob
+    else:
+        controller = controller_nonMPI(**controller_args, num_procs=num_procs)
+        P = controller.MS[0].levels[0].prob
+
+    t0 = 0.0 if t0 is None else t0
+    uinit = P.u_exact(t0) if u0 is None else u0
+
+    if fault_stuff is not None:
+        from pySDC.projects.Resilience.fault_injection import prepare_controller_for_faults
+
+        prepare_controller_for_faults(controller, fault_stuff)
+
+    crash = False
+    try:
+        uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
+    except ConvergenceError as e:
+        print(f'Warning: Premature termination!: {e}')
+        stats = controller.return_stats()
+        crash = True
+    return stats, controller, crash
+
+
+def u_exact(self, t, u_init=None, t_init=None, recompute=False):
+    import pickle
+    import os
+
+    path = f'data/stats/RBC-u_init-{t}.pickle'
+    if os.path.exists(path) and not recompute and t_init is None:
+        with open(path, 'rb') as file:
+            data = pickle.load(file)
+    else:
+        from pySDC.helpers.stats_helper import get_sorted
+        from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
+
+        convergence_controllers = {
+            Adaptivity: {'e_tol': 1e-8, 'dt_rel_min_slope': 0.25},
+        }
+        desc = {'convergence_controllers': convergence_controllers}
+
+        u0 = RayleighBenard()._u_exact(0) if u_init is None else u_init
+        t0 = 0 if t_init is None else t_init
+        stats, _, _ = run_RBC(Tend=t, u0=u0, t0=t0, custom_description=desc)
+
+        u = get_sorted(stats, type='u', recomputed=False)[-1]
+        data = u[1]
+
+        if t0 == 0:
+            with open(path, 'wb') as file:
+                pickle.dump(data, file)
+
+    return data
+
+
+if __name__ == '__main__':
+    from pySDC.implementations.hooks.log_errors import LogLocalErrorPostStep
+
+    stats, _, _ = run_RBC(hook_class=LogLocalErrorPostStep)

From 31eb63daee68f1fab9a97f9c05bf39d6e4b78c80 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Sat, 12 Oct 2024 14:13:44 +0200
Subject: [PATCH 02/14] Verified order and added configurations for
 dt-adaptivity and base strategy

---
 .../problem_classes/generic_spectral.py       |  22 +++
 pySDC/projects/Resilience/RBC.py              | 152 ++++++++++++++----
 pySDC/projects/Resilience/strategies.py       |  18 ++-
 3 files changed, 160 insertions(+), 32 deletions(-)

diff --git a/pySDC/implementations/problem_classes/generic_spectral.py b/pySDC/implementations/problem_classes/generic_spectral.py
index 04f88823fc..556ff46760 100644
--- a/pySDC/implementations/problem_classes/generic_spectral.py
+++ b/pySDC/implementations/problem_classes/generic_spectral.py
@@ -387,3 +387,25 @@ def compute_residual_DAE_MPI(self, stage=None):
     L.status.updated = False
 
     return None
+
+
+def get_extrapolated_error_DAE(self, S, **kwargs):
+    """
+    The extrapolation estimate combines values of u and f from multiple steps to extrapolate and compare to the
+    solution obtained by the time marching scheme. This function can be used in `EstimateExtrapolationError`.
+
+    Args:
+        S (pySDC.Step): The current step
+
+    Returns:
+        None
+    """
+    u_ex = self.get_extrapolated_solution(S)
+    diff_mask = S.levels[0].prob.diff_mask
+    if u_ex is not None:
+        S.levels[0].status.error_extrapolation_estimate = (
+            abs((u_ex - S.levels[0].u[-1])[diff_mask]) * self.coeff.prefactor
+        )
+        # print([abs(me) for me in (u_ex - S.levels[0].u[-1]) * self.coeff.prefactor])
+    else:
+        S.levels[0].status.error_extrapolation_estimate = None
diff --git a/pySDC/projects/Resilience/RBC.py b/pySDC/projects/Resilience/RBC.py
index 280283f557..a72bd8e5c3 100644
--- a/pySDC/projects/Resilience/RBC.py
+++ b/pySDC/projects/Resilience/RBC.py
@@ -1,14 +1,59 @@
 # script to run a Rayleigh-Benard Convection problem
-from pySDC.implementations.problem_classes.generic_spectral import compute_residual_DAE
+from pySDC.implementations.problem_classes.generic_spectral import compute_residual_DAE, get_extrapolated_error_DAE
 from pySDC.implementations.problem_classes.RayleighBenard import RayleighBenard
 from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
 from pySDC.projects.Resilience.hook import hook_collection, LogData
 from pySDC.projects.Resilience.strategies import merge_descriptions
 from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
 from pySDC.core.convergence_controller import ConvergenceController
+from pySDC.implementations.convergence_controller_classes.estimate_extrapolation_error import (
+    EstimateExtrapolationErrorNonMPI,
+)
 
 from pySDC.core.errors import ConvergenceError
 
+import numpy as np
+
+
+def u_exact(self, t, u_init=None, t_init=None, recompute=False):
+    import pickle
+    import os
+
+    path = f'data/stats/RBC-u_init-{t:.8f}.pickle'
+    if os.path.exists(path) and not recompute and t_init is None:
+        with open(path, 'rb') as file:
+            data = pickle.load(file)
+    else:
+        from pySDC.helpers.stats_helper import get_sorted
+        from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
+
+        convergence_controllers = {
+            Adaptivity: {'e_tol': 1e-8, 'dt_rel_min_slope': 0.25},
+        }
+        desc = {'convergence_controllers': convergence_controllers}
+
+        u0 = self._u_exact(0) if u_init is None else u_init
+        t0 = 0 if t_init is None else t_init
+
+        if t == t0:
+            return u0
+        else:
+            stats, _, _ = run_RBC(Tend=t, u0=u0, t0=t0, custom_description=desc)
+
+        u = get_sorted(stats, type='u', recomputed=False)[-1]
+        data = u[1]
+
+        if t0 == 0:
+            with open(path, 'wb') as file:
+                pickle.dump(data, file)
+
+    return data
+
+
+RayleighBenard._u_exact = RayleighBenard.u_exact
+RayleighBenard.u_exact = u_exact
+EstimateExtrapolationErrorNonMPI.get_extrapolated_error = get_extrapolated_error_DAE
+
 
 class ReachTendExactly(ConvergenceController):
 
@@ -21,7 +66,9 @@ def setup(self, controller, params, description, **kwargs):
 
     def get_new_step_size(self, controller, step, **kwargs):
         L = step.levels[0]
-        L.status.dt_new = min([L.params.dt, self.params.Tend - L.time - L.dt])
+        dt = L.status.dt_new if L.status.dt_new else L.params.dt
+        if self.params.Tend - L.time - L.dt < dt:
+            L.status.dt_new = min([dt, self.params.Tend - L.time - L.dt])
 
 
 def run_RBC(
@@ -32,7 +79,7 @@ def run_RBC(
     fault_stuff=None,
     custom_controller_params=None,
     u0=None,
-    t0=10,
+    t0=11,
     use_MPI=False,
     **kwargs,
 ):
@@ -54,7 +101,7 @@ def run_RBC(
         bool: If the code crashed
     """
     level_params = {}
-    level_params['dt'] = 5e-4
+    level_params['dt'] = 1e-3
     level_params['restol'] = -1
 
     sweeper_params = {}
@@ -63,7 +110,9 @@ def run_RBC(
     sweeper_params['QI'] = 'LU'
     sweeper_params['QE'] = 'PIC'
 
-    problem_params = {}
+    from mpi4py import MPI
+
+    problem_params = {'comm': MPI.COMM_SELF}
 
     step_params = {}
     step_params['maxiter'] = 5
@@ -81,9 +130,6 @@ def run_RBC(
 
     imex_1st_order.compute_residual = compute_residual_DAE
 
-    RayleighBenard._u_exact = RayleighBenard.u_exact
-    RayleighBenard.u_exact = u_exact
-
     description = {}
     description['problem_class'] = RayleighBenard
     description['problem_params'] = problem_params
@@ -91,6 +137,7 @@ def run_RBC(
     description['sweeper_params'] = sweeper_params
     description['level_params'] = level_params
     description['step_params'] = step_params
+    description['convergence_controllers'] = convergence_controllers
 
     if custom_description is not None:
         description = merge_descriptions(description, custom_description)
@@ -128,38 +175,81 @@ def run_RBC(
     return stats, controller, crash
 
 
-def u_exact(self, t, u_init=None, t_init=None, recompute=False):
+def generate_data_for_fault_stats():
+    prob = RayleighBenard()
+    for t in [11.0, 14.0]:
+        prob.u_exact(t)
+
+
+def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recompute=False):
+    from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun
+    from pySDC.implementations.hooks.log_work import LogSDCIterations, LogWork
+    from pySDC.helpers.stats_helper import get_sorted
     import pickle
     import os
 
-    path = f'data/stats/RBC-u_init-{t}.pickle'
-    if os.path.exists(path) and not recompute and t_init is None:
-        with open(path, 'rb') as file:
-            data = pickle.load(file)
-    else:
-        from pySDC.helpers.stats_helper import get_sorted
-        from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
+    sweeper_params = {'num_nodes': num_nodes}
+    level_params = {'e_tol': e_tol, 'restol': restol}
+    step_params = {'maxiter': 99}
 
-        convergence_controllers = {
-            Adaptivity: {'e_tol': 1e-8, 'dt_rel_min_slope': 0.25},
-        }
-        desc = {'convergence_controllers': convergence_controllers}
+    errors = []
+    dts = []
+    for i in range(steps):
+        dts += [dt / 2**i]
 
-        u0 = RayleighBenard()._u_exact(0) if u_init is None else u_init
-        t0 = 0 if t_init is None else t_init
-        stats, _, _ = run_RBC(Tend=t, u0=u0, t0=t0, custom_description=desc)
+        path = f'data/stats/RBC-u_order-{t:.8f}-{dts[-1]:.8e}-{num_nodes}-{e_tol:.2e}-{restol:.2e}.pickle'
 
-        u = get_sorted(stats, type='u', recomputed=False)[-1]
-        data = u[1]
+        if os.path.exists(path) and not recompute:
+            with open(path, 'rb') as file:
+                stats = pickle.load(file)
+        else:
+
+            level_params['dt'] = dts[-1]
+
+            desc = {'sweeper_params': sweeper_params, 'level_params': level_params, 'step_params': step_params}
+
+            stats, _, _ = run_RBC(
+                Tend=t + dt,
+                t0=t,
+                custom_description=desc,
+                hook_class=[LogGlobalErrorPostRun, LogSDCIterations, LogWork],
+            )
 
-        if t0 == 0:
             with open(path, 'wb') as file:
-                pickle.dump(data, file)
+                pickle.dump(stats, file)
 
-    return data
+        e = get_sorted(stats, type='e_global_post_run')
+        # k = get_sorted(stats, type='k')
 
+        errors += [e[-1][1]]
 
-if __name__ == '__main__':
-    from pySDC.implementations.hooks.log_errors import LogLocalErrorPostStep
+    errors = np.array(errors)
+    dts = np.array(dts)
+    mask = np.isfinite(errors)
+    max_error = np.nanmax(errors)
+
+    errors = errors[mask]
+    dts = dts[mask]
+    ax.loglog(dts, errors, label=f'{num_nodes} nodes')
+    ax.loglog(
+        dts, [max_error * (me / dts[0]) ** (2 * num_nodes - 1) for me in dts], ls='--', label=f'order {2*num_nodes-1}'
+    )
+    ax.set_xlabel(r'$\Delta t$')
+    ax.set_ylabel(r'global error')
+    ax.legend(frameon=False)
+
+
+def test_order(t=14, dt=1e-1, steps=6):
 
-    stats, _, _ = run_RBC(hook_class=LogLocalErrorPostStep)
+    import matplotlib.pyplot as plt
+
+    fig, ax = plt.subplots(1, 1)
+    for num_nodes in [1, 2, 3, 4]:
+        plot_order(t=t, dt=dt, steps=steps, num_nodes=num_nodes, ax=ax)
+    plt.show()
+
+
+if __name__ == '__main__':
+    generate_data_for_fault_stats()
+    test_order()
+    # stats, _, _ = run_RBC()
diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py
index 2d4edaa689..f7c2e67fe5 100644
--- a/pySDC/projects/Resilience/strategies.py
+++ b/pySDC/projects/Resilience/strategies.py
@@ -130,6 +130,8 @@ def get_fault_args(self, problem, num_procs):
             args['time'] = 0.3
         elif problem.__name__ == "run_AC":
             args['time'] = 1e-2
+        elif problem.__name__ == "run_RBC":
+            args['time'] = 11.5
 
         return args
 
@@ -150,7 +152,7 @@ def get_random_params(self, problem, num_procs):
         rnd_params['iteration'] = base_params['step_params']['maxiter']
         rnd_params['rank'] = num_procs
 
-        if problem.__name__ in ['run_Schroedinger', 'run_quench', 'run_AC']:
+        if problem.__name__ in ['run_Schroedinger', 'run_quench', 'run_AC', 'run_RBC']:
             rnd_params['min_node'] = 1
 
         if problem.__name__ == "run_quench":
@@ -206,6 +208,8 @@ def get_Tend(cls, problem, num_procs=1):
             return 500.0
         elif problem.__name__ == "run_AC":
             return 0.025
+        elif problem.__name__ == "run_RBC":
+            return 14
         else:
             raise NotImplementedError('I don\'t have a final time for your problem!')
 
@@ -269,6 +273,9 @@ def get_base_parameters(self, problem, num_procs=1):
                 'nu': 2,
             }
             custom_description['level_params'] = {'restol': -1, 'dt': 0.1 * eps**2}
+        elif problem.__name__ == 'run_RBC':
+            custom_description['level_params']['dt'] = 5e-3
+            custom_description['step_params'] = {'maxiter': 5}
 
         custom_description['convergence_controllers'] = {
             # StepSizeLimiter: {'dt_min': self.get_Tend(problem=problem, num_procs=num_procs) / self.max_steps}
@@ -287,6 +294,7 @@ def get_base_parameters(self, problem, num_procs=1):
             'run_Schroedinger': 150,
             'run_quench': 150,
             'run_AC': 150,
+            'run_RBC': 500,
         }
 
         custom_description['convergence_controllers'][StopAtMaxRuntime] = {
@@ -521,6 +529,7 @@ def get_custom_description(self, problem, num_procs):
 
         dt_max = np.inf
         dt_slope_max = np.inf
+        dt_slope_min = 0
 
         if problem.__name__ == "run_piline":
             e_tol = 1e-7
@@ -545,6 +554,9 @@ def get_custom_description(self, problem, num_procs):
         elif problem.__name__ == "run_AC":
             e_tol = 1e-7
             # dt_max = 0.1 * base_params['problem_params']['eps'] ** 2
+        elif problem.__name__ == 'run_RBC':
+            e_tol = 1e-3
+            dt_slope_min = 0.25
 
         else:
             raise NotImplementedError(
@@ -555,6 +567,7 @@ def get_custom_description(self, problem, num_procs):
         custom_description['convergence_controllers'][Adaptivity] = {
             'e_tol': e_tol,
             'dt_slope_max': dt_slope_max,
+            'dt_rel_min_slope': dt_slope_min,
         }
         custom_description['convergence_controllers'][StepSizeLimiter] = {
             'dt_max': dt_max,
@@ -928,6 +941,9 @@ def get_custom_description(self, problem, num_procs):
         elif problem.__name__ == 'run_AC':
             HotRod_tol = 9.564437e-06
             maxiter = 6
+        elif problem.__name__ == 'run_RBC':
+            HotRod_tol = 1e0  # 1.#4e-1
+            maxiter = 6
         else:
             raise NotImplementedError(
                 'I don\'t have a tolerance for Hot Rod for your problem. Please add one to the\

From 71ac6c7fa8750e7b662c63fbd6bbe94da4d86712 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Wed, 16 Oct 2024 08:51:03 +0200
Subject: [PATCH 03/14] Implemented configurations for faults in RBC

---
 pySDC/helpers/plot_helper.py             |  2 ++
 pySDC/projects/Resilience/RBC.py         | 25 ++++++++++------
 pySDC/projects/Resilience/README.rst     | 21 +++++++++----
 pySDC/projects/Resilience/fault_stats.py |  9 +++++-
 pySDC/projects/Resilience/paper_plots.py | 38 ++++++++++++++++++++----
 pySDC/projects/Resilience/strategies.py  | 27 +++++++++++++----
 6 files changed, 97 insertions(+), 25 deletions(-)

diff --git a/pySDC/helpers/plot_helper.py b/pySDC/helpers/plot_helper.py
index 51e8783b96..8ddea7cd0c 100644
--- a/pySDC/helpers/plot_helper.py
+++ b/pySDC/helpers/plot_helper.py
@@ -43,11 +43,13 @@ def figsize_by_journal(journal, scale, ratio):  # pragma: no cover
         'JSC_beamer': 426.79135,
         'Springer_Numerical_Algorithms': 338.58778,
         'JSC_thesis': 434.26027,
+        'TUHH_thesis': 426.79135,
     }
     # store text height in points here, get this from LaTeX using \the\textheight
     textheights = {
         'JSC_beamer': 214.43411,
         'JSC_thesis': 635.5,
+        'TUHH_thesis': 631.65118,
     }
     assert (
         journal in textwidths.keys()
diff --git a/pySDC/projects/Resilience/RBC.py b/pySDC/projects/Resilience/RBC.py
index a72bd8e5c3..87923b0abe 100644
--- a/pySDC/projects/Resilience/RBC.py
+++ b/pySDC/projects/Resilience/RBC.py
@@ -28,7 +28,7 @@ def u_exact(self, t, u_init=None, t_init=None, recompute=False):
         from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
 
         convergence_controllers = {
-            Adaptivity: {'e_tol': 1e-8, 'dt_rel_min_slope': 0.25},
+            Adaptivity: {'e_tol': 1e-8, 'dt_rel_min_slope': 0.25, 'dt_min': 1e-5},
         }
         desc = {'convergence_controllers': convergence_controllers}
 
@@ -54,6 +54,8 @@ def u_exact(self, t, u_init=None, t_init=None, recompute=False):
 RayleighBenard.u_exact = u_exact
 EstimateExtrapolationErrorNonMPI.get_extrapolated_error = get_extrapolated_error_DAE
 
+PROBLEM_PARAMS = {'Rayleigh': 2e4, 'nx': 256, 'nz': 128}
+
 
 class ReachTendExactly(ConvergenceController):
 
@@ -61,25 +63,27 @@ def setup(self, controller, params, description, **kwargs):
         defaults = {
             "control_order": +50,
             "Tend": None,
+            'min_step_size': 1e-9,
         }
         return {**defaults, **super().setup(controller, params, description, **kwargs)}
 
     def get_new_step_size(self, controller, step, **kwargs):
         L = step.levels[0]
         dt = L.status.dt_new if L.status.dt_new else L.params.dt
-        if self.params.Tend - L.time - L.dt < dt:
-            L.status.dt_new = min([dt, self.params.Tend - L.time - L.dt])
+        time_left = self.params.Tend - L.time - L.dt
+        if time_left < dt:
+            L.status.dt_new = min([dt, max([time_left, self.params.min_step_size])])
 
 
 def run_RBC(
     custom_description=None,
     num_procs=1,
-    Tend=14.0,
+    Tend=21.0,
     hook_class=LogData,
     fault_stuff=None,
     custom_controller_params=None,
     u0=None,
-    t0=11,
+    t0=20.0,
     use_MPI=False,
     **kwargs,
 ):
@@ -112,7 +116,7 @@ def run_RBC(
 
     from mpi4py import MPI
 
-    problem_params = {'comm': MPI.COMM_SELF}
+    problem_params = {'comm': MPI.COMM_SELF, **PROBLEM_PARAMS}
 
     step_params = {}
     step_params['maxiter'] = 5
@@ -176,8 +180,11 @@ def run_RBC(
 
 
 def generate_data_for_fault_stats():
-    prob = RayleighBenard()
-    for t in [11.0, 14.0]:
+    prob = RayleighBenard(**PROBLEM_PARAMS)
+    for t in [
+        20,
+        21,
+    ]:
         prob.u_exact(t)
 
 
@@ -251,5 +258,5 @@ def test_order(t=14, dt=1e-1, steps=6):
 
 if __name__ == '__main__':
     generate_data_for_fault_stats()
-    test_order()
+    # test_order()
     # stats, _, _ = run_RBC()
diff --git a/pySDC/projects/Resilience/README.rst b/pySDC/projects/Resilience/README.rst
index becbd2742f..3092c023b2 100644
--- a/pySDC/projects/Resilience/README.rst
+++ b/pySDC/projects/Resilience/README.rst
@@ -50,10 +50,21 @@ Then, navigate to this directory, `pySDC/projects/Resilience/` and run the follo
 .. code-block:: bash
  
     mpirun -np 4 python work_precision.py
-    mpirun -np 4 python fault_stats.py prob run_vdp
-    mpirun -np 4 python fault_stats.py prob run_quench
-    mpirun -np 4 python fault_stats.py prob run_AC
-    mpirun -np 4 python fault_stats.py prob run_Schroedinger
-    python paper_plots.py
+    python paper_plots.py --target=adaptivity
 
 Possibly, you need to create some directories in this one to store and load things, if path errors occur.
+
+Reproduction of the plots in the resilience paper
+-------------------------------------------------
+To reproduce the plots you need to install pySDC using this project's `environment.yml` file, which is in the same directory as this README.
+
+.. code-block:: bash
+ 
+    mpirun -np 4 python work_precision.py
+    mpirun -np 4 python fault_stats.py prob run_Lorenz
+    mpirun -np 4 python fault_stats.py prob run_Schroedinger
+    mpirun -np 4 python fault_stats.py prob run_AC
+    mpirun -np 4 python fault_stats.py prob run_RBC
+    python paper_plots.py --target=resilience
+
+Please be aware that generating the fault data for Rayleigh-Benard requires generating reference solutions, which may take several hours.
diff --git a/pySDC/projects/Resilience/fault_stats.py b/pySDC/projects/Resilience/fault_stats.py
index f28aa11746..7b22317253 100644
--- a/pySDC/projects/Resilience/fault_stats.py
+++ b/pySDC/projects/Resilience/fault_stats.py
@@ -21,6 +21,7 @@
 from pySDC.projects.Resilience.Schroedinger import run_Schroedinger
 from pySDC.projects.Resilience.quench import run_quench
 from pySDC.projects.Resilience.AC import run_AC
+from pySDC.projects.Resilience.RBC import run_RBC
 
 from pySDC.projects.Resilience.strategies import BaseStrategy, AdaptivityStrategy, IterateStrategy, HotRodStrategy
 import logging
@@ -632,6 +633,8 @@ def get_name(self, strategy=None, faults=True, mode=None):
             prob_name = 'Quench'
         elif self.prob.__name__ == 'run_AC':
             prob_name = 'Allen-Cahn'
+        elif self.prob.__name__ == 'run_RBC':
+            prob_name = 'Rayleigh-Benard'
         else:
             raise NotImplementedError(f'Name not implemented for problem {self.prob}')
 
@@ -1565,6 +1568,10 @@ def parse_args():
                 kwargs['prob'] = run_Schroedinger
             elif sys.argv[i + 1] == 'run_quench':
                 kwargs['prob'] = run_quench
+            elif sys.argv[i + 1] == 'run_AC':
+                kwargs['prob'] = run_AC
+            elif sys.argv[i + 1] == 'run_RBC':
+                kwargs['prob'] = run_RBC
             else:
                 raise NotImplementedError
         elif 'num_procs' in sys.argv[i]:
@@ -1654,7 +1661,7 @@ def compare_adaptivity_modes():
 
 def main():
     kwargs = {
-        'prob': run_AC,
+        'prob': run_RBC,
         'num_procs': 1,
         'mode': 'default',
         'runs': 2000,
diff --git a/pySDC/projects/Resilience/paper_plots.py b/pySDC/projects/Resilience/paper_plots.py
index e3ad6d435e..093ca29dbc 100644
--- a/pySDC/projects/Resilience/paper_plots.py
+++ b/pySDC/projects/Resilience/paper_plots.py
@@ -9,6 +9,7 @@
     run_vdp,
     run_quench,
     run_AC,
+    run_RBC,
     RECOVERY_THRESH_ABS,
 )
 from pySDC.projects.Resilience.strategies import (
@@ -196,7 +197,7 @@ def compare_recovery_rate_problems(**kwargs):  # pragma: no cover
         None
     """
     stats = [
-        get_stats(run_vdp, **kwargs),
+        get_stats(run_RBC, **kwargs),
         get_stats(run_quench, **kwargs),
         get_stats(run_Schroedinger, **kwargs),
         get_stats(run_AC, **kwargs),
@@ -614,8 +615,35 @@ def make_plots_for_notes():  # pragma: no cover
     analyse_resilience(run_quench, format='png')
 
 
+def make_plots_for_thesis():  # pragma: no cover
+    global JOURNAL
+    JOURNAL = 'TUHH_thesis'
+
+    # plot_adaptivity_stuff()
+    # plot_vdp_solution()
+    compare_recovery_rate_problems(num_procs=1, strategy_type='SDC')
+
+
 if __name__ == "__main__":
-    # make_plots_for_notes()
-    # make_plots_for_SIAM_CSE23()
-    # make_plots_for_TIME_X_website()
-    make_plots_for_adaptivity_paper()
+    import argparse
+
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        '--target', choices=['adaptivity', 'resilience', 'thesis', 'notes', 'SIAM_CSE23', 'TIME_X_website'], type=str
+    )
+    args = parser.parse_args()
+
+    if args.target == 'adaptivity':
+        make_plots_for_adaptivity_paper()
+    elif args.target == 'resilience':
+        make_plots_for_resilience_paper()
+    elif args.target == 'thesis':
+        make_plots_for_thesis()
+    elif args.target == 'notes':
+        make_plots_for_notes()
+    elif args.target == 'SIAM_CSE23':
+        make_plots_for_SIAM_CSE23()
+    elif args.target == 'TIME_X_website':
+        make_plots_for_TIME_X_website()
+    else:
+        raise NotImplementedError(f'Don\'t know how to make plots for target {args.target}')
diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py
index f7c2e67fe5..74a056bcf7 100644
--- a/pySDC/projects/Resilience/strategies.py
+++ b/pySDC/projects/Resilience/strategies.py
@@ -131,7 +131,7 @@ def get_fault_args(self, problem, num_procs):
         elif problem.__name__ == "run_AC":
             args['time'] = 1e-2
         elif problem.__name__ == "run_RBC":
-            args['time'] = 11.5
+            args['time'] = 20.09
 
         return args
 
@@ -209,7 +209,7 @@ def get_Tend(cls, problem, num_procs=1):
         elif problem.__name__ == "run_AC":
             return 0.025
         elif problem.__name__ == "run_RBC":
-            return 14
+            return 21
         else:
             raise NotImplementedError('I don\'t have a final time for your problem!')
 
@@ -274,7 +274,7 @@ def get_base_parameters(self, problem, num_procs=1):
             }
             custom_description['level_params'] = {'restol': -1, 'dt': 0.1 * eps**2}
         elif problem.__name__ == 'run_RBC':
-            custom_description['level_params']['dt'] = 5e-3
+            custom_description['level_params']['dt'] = 5e-2
             custom_description['step_params'] = {'maxiter': 5}
 
         custom_description['convergence_controllers'] = {
@@ -381,7 +381,7 @@ def get_custom_description(self, problem, num_procs=1):
             'maxiter': 15,
         }
 
-        if self.newton_inexactness and problem.__name__ not in ['run_Schroedinger', 'run_AC']:
+        if self.newton_inexactness and problem.__name__ not in ['run_Schroedinger', 'run_AC', 'run_RBC']:
             if problem.__name__ == 'run_quench':
                 inexactness_params['ratio'] = 1e-1
                 inexactness_params['min_tol'] = 1e-11
@@ -769,6 +769,8 @@ def get_custom_description(self, problem, num_procs):
             restol = 1e-7
         elif problem.__name__ == "run_AC":
             restol = 1e-11
+        elif problem.__name__ == "run_RBC":
+            restol = 1e-4
         else:
             raise NotImplementedError(
                 'I don\'t have a residual tolerance for your problem. Please add one to the \
@@ -840,6 +842,8 @@ def get_custom_description(self, problem, num_procs, *args, **kwargs):
         elif problem.__name__ == "run_AC":
             desc['level_params']['restol'] = 1e-11
             desc['level_params']['dt'] = 0.4 * desc['problem_params']['eps'] ** 2 / 8.0
+        elif problem.__name__ == "run_RBC":
+            desc['level_params']['restol'] = 1e-4
         return desc
 
     def get_custom_description_for_faults(self, problem, *args, **kwargs):
@@ -942,7 +946,7 @@ def get_custom_description(self, problem, num_procs):
             HotRod_tol = 9.564437e-06
             maxiter = 6
         elif problem.__name__ == 'run_RBC':
-            HotRod_tol = 1e0  # 1.#4e-1
+            HotRod_tol = 6.34e-6
             maxiter = 6
         else:
             raise NotImplementedError(
@@ -1880,7 +1884,10 @@ def get_custom_description(self, problem, num_procs):
         dt_max = np.inf
         restol_rel = 1e-4
         restol_min = 1e-12
+        restol_max = 1e-5
+        dt_slope_min = 0
         dt_min = 0
+        abort_at_growing_residual = True
         level_params = {}
         problem_params = {}
 
@@ -1906,6 +1913,13 @@ def get_custom_description(self, problem, num_procs):
             e_tol = 1.0e-4
             restol_rel = 1e-3
             # dt_max = 0.1 * base_params['problem_params']['eps'] ** 2
+        elif problem.__name__ == "run_RBC":
+            e_tol = 1e-3
+            dt_slope_min = 0.25
+            abort_at_growing_residual = False
+            restol_rel = 1e-4
+            restol_max = 1e-1
+            restol_min = 1e-6
         else:
             raise NotImplementedError(
                 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\
@@ -1917,14 +1931,17 @@ def get_custom_description(self, problem, num_procs):
                 'e_tol': e_tol,
                 'restol_rel': restol_rel if self.use_restol_rel else 1e-11,
                 'restol_min': restol_min if self.use_restol_rel else 1e-12,
+                'restol_max': restol_max if self.use_restol_rel else 1e-5,
                 'restart_at_maxiter': True,
                 'factor_if_not_converged': self.max_slope,
                 'interpolate_between_restarts': self.interpolate_between_restarts,
+                'abort_at_growing_residual': abort_at_growing_residual,
             },
             StepSizeLimiter: {
                 'dt_max': dt_max,
                 'dt_slope_max': self.max_slope,
                 'dt_min': dt_min,
+                'dt_rel_min_slope': dt_slope_min,
             },
         }
         custom_description['level_params'] = level_params

From fccd0450ad8ffb9ec01592d03c4e185061a1b2f4 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Wed, 16 Oct 2024 10:02:02 +0200
Subject: [PATCH 04/14] Added plot of RBC solution

---
 pySDC/projects/Resilience/paper_plots.py | 67 +++++++++++++++++++-----
 1 file changed, 55 insertions(+), 12 deletions(-)

diff --git a/pySDC/projects/Resilience/paper_plots.py b/pySDC/projects/Resilience/paper_plots.py
index 093ca29dbc..f746c9c3d8 100644
--- a/pySDC/projects/Resilience/paper_plots.py
+++ b/pySDC/projects/Resilience/paper_plots.py
@@ -188,21 +188,24 @@ def plot_recovery_rate_recoverable_only(stats_analyser, fig, ax, **kwargs):  # p
         )
 
 
-def compare_recovery_rate_problems(**kwargs):  # pragma: no cover
+def compare_recovery_rate_problems(target='resilience', **kwargs):  # pragma: no cover
     """
-    Compare the recovery rate for vdP, Lorenz and Schroedinger problems.
+    Compare the recovery rate for various problems.
     Only faults that can be recovered are shown.
 
     Returns:
         None
     """
-    stats = [
-        get_stats(run_RBC, **kwargs),
-        get_stats(run_quench, **kwargs),
-        get_stats(run_Schroedinger, **kwargs),
-        get_stats(run_AC, **kwargs),
-    ]
-    titles = ['Van der Pol', 'Quench', r'Schr\"odinger', 'Allen-Cahn']
+    if target == 'resilience':
+        problems = [run_Lorenz, run_Schroedinger, run_AC, run_RBC]
+        titles = ['Lorenz', r'Schr\"odinger', 'Allen-Cahn', 'Rayleigh-Benard']
+    elif target == 'thesis':
+        problems = [run_vdp, run_Lorenz, run_AC, run_RBC]  # TODO: swap in Gray-Scott
+        titles = ['Van der Pol', 'Lorenz', 'Allen-Cahn', 'Rayleigh-Benard']
+    else:
+        raise NotImplementedError()
+
+    stats = [get_stats(problem, **kwargs) for problem in problems]
 
     my_setup_mpl()
     fig, axs = plt.subplots(2, 2, figsize=figsize_by_journal(JOURNAL, 1, 0.8), sharey=True)
@@ -421,6 +424,42 @@ def plot_quench_solution():  # pragma: no cover
     savefig(fig, 'quench_sol')
 
 
+def plot_RBC_solution():  # pragma: no cover
+    """
+    Plot solution of Rayleigh-Benard convection
+    """
+    my_setup_mpl()
+
+    from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+    plt.rcParams['figure.constrained_layout.use'] = True
+    fig, axs = plt.subplots(2, 1, sharex=True, sharey=True, figsize=figsize_by_journal(JOURNAL, 1.0, 0.45))
+    caxs = []
+    divider = make_axes_locatable(axs[0])
+    caxs += [divider.append_axes('right', size='3%', pad=0.03)]
+    divider2 = make_axes_locatable(axs[1])
+    caxs += [divider2.append_axes('right', size='3%', pad=0.03)]
+
+    from pySDC.projects.Resilience.RBC import RayleighBenard, PROBLEM_PARAMS
+
+    prob = RayleighBenard(**PROBLEM_PARAMS)
+
+    def _plot(t, ax, cax):
+        u_hat = prob.u_exact(t)
+        u = prob.itransform(u_hat)
+        im = ax.pcolormesh(prob.X, prob.Z, u[prob.index('T')], rasterized=True)
+        fig.colorbar(im, cax, label=f'$T(t={{{t}}})$')
+
+    _plot(0, axs[0], caxs[0])
+    _plot(21, axs[1], caxs[1])
+
+    axs[1].set_xlabel('$x$')
+    axs[0].set_ylabel('$z$')
+    axs[1].set_ylabel('$z$')
+
+    savefig(fig, 'RBC_sol', tight_layout=False)
+
+
 def plot_Schroedinger_solution():  # pragma: no cover
     from pySDC.implementations.problem_classes.NonlinearSchroedinger_MPIFFT import nonlinearschroedinger_imex
 
@@ -600,7 +639,9 @@ def make_plots_for_resilience_paper():  # pragma: no cover
     plot_recovery_rate(get_stats(run_vdp))
     plot_fault_vdp(0)
     plot_fault_vdp(13)
-    compare_recovery_rate_problems(num_procs=1, strategy_type='SDC')
+    compare_recovery_rate_problems(target='resilience', num_procs=1, strategy_type='SDC')
+
+    plot_RBC_solution()
 
 
 def make_plots_for_notes():  # pragma: no cover
@@ -619,9 +660,11 @@ def make_plots_for_thesis():  # pragma: no cover
     global JOURNAL
     JOURNAL = 'TUHH_thesis'
 
-    # plot_adaptivity_stuff()
+    plot_RBC_solution()
     # plot_vdp_solution()
-    compare_recovery_rate_problems(num_procs=1, strategy_type='SDC')
+
+    # plot_adaptivity_stuff()
+    compare_recovery_rate_problems(target='thesis', num_procs=1, strategy_type='SDC')
 
 
 if __name__ == "__main__":

From 4d5982a73ba7d223aadcb78b9be48b1d035dc1e9 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Wed, 16 Oct 2024 15:15:27 +0200
Subject: [PATCH 05/14] Limit fault targets in space index

---
 pySDC/projects/Resilience/fault_stats.py | 4 ++--
 pySDC/projects/Resilience/strategies.py  | 3 +++
 2 files changed, 5 insertions(+), 2 deletions(-)

diff --git a/pySDC/projects/Resilience/fault_stats.py b/pySDC/projects/Resilience/fault_stats.py
index 7b22317253..1d0b7db909 100644
--- a/pySDC/projects/Resilience/fault_stats.py
+++ b/pySDC/projects/Resilience/fault_stats.py
@@ -1665,7 +1665,7 @@ def main():
         'num_procs': 1,
         'mode': 'default',
         'runs': 2000,
-        'reload': True,
+        'reload': False,
         **parse_args(),
     }
 
@@ -1689,7 +1689,7 @@ def main():
         stats_path='data/stats-jusuf',
         **kwargs,
     )
-    stats_analyser.run_stats_generation(runs=kwargs['runs'], step=12)
+    stats_analyser.run_stats_generation(runs=kwargs['runs'], step=24)
 
     if MPI.COMM_WORLD.rank > 0:  # make sure only one rank accesses the data
         return None
diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py
index 74a056bcf7..1f3cc2e3a8 100644
--- a/pySDC/projects/Resilience/strategies.py
+++ b/pySDC/projects/Resilience/strategies.py
@@ -159,6 +159,9 @@ def get_random_params(self, problem, num_procs):
             rnd_params['iteration'] = 5
         elif problem.__name__ == 'run_Lorenz':
             rnd_params['iteration'] = 5
+        elif problem.__name__ == 'run_RBC':
+            rnd_params['problem_pos'] = [3, 16, 16]
+
         return rnd_params
 
     @property

From 785aac42c96a27b913ae7a286ae25c2a3834e2cc Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Wed, 16 Oct 2024 15:35:11 +0200
Subject: [PATCH 06/14] Fix

---
 pySDC/projects/Resilience/RBC.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/pySDC/projects/Resilience/RBC.py b/pySDC/projects/Resilience/RBC.py
index 87923b0abe..3502242535 100644
--- a/pySDC/projects/Resilience/RBC.py
+++ b/pySDC/projects/Resilience/RBC.py
@@ -52,7 +52,6 @@ def u_exact(self, t, u_init=None, t_init=None, recompute=False):
 
 RayleighBenard._u_exact = RayleighBenard.u_exact
 RayleighBenard.u_exact = u_exact
-EstimateExtrapolationErrorNonMPI.get_extrapolated_error = get_extrapolated_error_DAE
 
 PROBLEM_PARAMS = {'Rayleigh': 2e4, 'nx': 256, 'nz': 128}
 
@@ -104,6 +103,8 @@ def run_RBC(
         controller: The controller
         bool: If the code crashed
     """
+    EstimateExtrapolationErrorNonMPI.get_extrapolated_error = get_extrapolated_error_DAE
+
     level_params = {}
     level_params['dt'] = 1e-3
     level_params['restol'] = -1

From b57255611274d0879d130874c3dce1e79848d943 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Wed, 16 Oct 2024 16:59:29 +0200
Subject: [PATCH 07/14] Fixed fault insertion time for Hot Rod

---
 pySDC/projects/Resilience/strategies.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py
index 1f3cc2e3a8..0dadc4a34f 100644
--- a/pySDC/projects/Resilience/strategies.py
+++ b/pySDC/projects/Resilience/strategies.py
@@ -131,7 +131,7 @@ def get_fault_args(self, problem, num_procs):
         elif problem.__name__ == "run_AC":
             args['time'] = 1e-2
         elif problem.__name__ == "run_RBC":
-            args['time'] = 20.09
+            args['time'] = 20.19
 
         return args
 

From d53717b27d2b54dcf69bdbd185747ca4858517a6 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Thu, 17 Oct 2024 09:47:38 +0200
Subject: [PATCH 08/14] Changed colormap for RBC plots in the paper

---
 pySDC/projects/Resilience/paper_plots.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/pySDC/projects/Resilience/paper_plots.py b/pySDC/projects/Resilience/paper_plots.py
index f746c9c3d8..04bea2d954 100644
--- a/pySDC/projects/Resilience/paper_plots.py
+++ b/pySDC/projects/Resilience/paper_plots.py
@@ -447,7 +447,7 @@ def plot_RBC_solution():  # pragma: no cover
     def _plot(t, ax, cax):
         u_hat = prob.u_exact(t)
         u = prob.itransform(u_hat)
-        im = ax.pcolormesh(prob.X, prob.Z, u[prob.index('T')], rasterized=True)
+        im = ax.pcolormesh(prob.X, prob.Z, u[prob.index('T')], rasterized=True, cmap='plasma')
         fig.colorbar(im, cax, label=f'$T(t={{{t}}})$')
 
     _plot(0, axs[0], caxs[0])
@@ -636,12 +636,11 @@ def make_plots_for_adaptivity_paper():  # pragma: no cover
 
 
 def make_plots_for_resilience_paper():  # pragma: no cover
+    compare_recovery_rate_problems(target='resilience', num_procs=1, strategy_type='SDC')
+    plot_RBC_solution()
     plot_recovery_rate(get_stats(run_vdp))
     plot_fault_vdp(0)
     plot_fault_vdp(13)
-    compare_recovery_rate_problems(target='resilience', num_procs=1, strategy_type='SDC')
-
-    plot_RBC_solution()
 
 
 def make_plots_for_notes():  # pragma: no cover

From 3d2c9948961e946a9b03b52540c8cc85ec50d7a7 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Fri, 18 Oct 2024 13:39:37 +0200
Subject: [PATCH 09/14] Fixed bug in reference solution

---
 .../sweeper_classes/Runge_Kutta.py            |  9 ++-
 pySDC/projects/Resilience/RBC.py              | 56 +++++++++++++------
 2 files changed, 46 insertions(+), 19 deletions(-)

diff --git a/pySDC/implementations/sweeper_classes/Runge_Kutta.py b/pySDC/implementations/sweeper_classes/Runge_Kutta.py
index 801beb9620..a0805f1f2f 100644
--- a/pySDC/implementations/sweeper_classes/Runge_Kutta.py
+++ b/pySDC/implementations/sweeper_classes/Runge_Kutta.py
@@ -441,9 +441,12 @@ def update_nodes(self):
                 rhs += lvl.dt * (self.QI[m + 1, j] * lvl.f[j].impl + self.QE[m + 1, j] * lvl.f[j].expl)
 
             # implicit solve with prefactor stemming from the diagonal of Qd, use previous stage as initial guess
-            lvl.u[m + 1][:] = prob.solve_system(
-                rhs, lvl.dt * self.QI[m + 1, m + 1], lvl.u[m], lvl.time + lvl.dt * self.coll.nodes[m + 1]
-            )
+            if self.QI[m + 1, m + 1] != 0:
+                lvl.u[m + 1][:] = prob.solve_system(
+                    rhs, lvl.dt * self.QI[m + 1, m + 1], lvl.u[m], lvl.time + lvl.dt * self.coll.nodes[m + 1]
+                )
+            else:
+                lvl.u[m + 1][:] = rhs[:]
 
             # update function values (we don't usually need to evaluate the RHS at the solution of the step)
             if m < M - self.coll.num_solution_stages or self.params.eval_rhs_at_right_boundary:
diff --git a/pySDC/projects/Resilience/RBC.py b/pySDC/projects/Resilience/RBC.py
index 3502242535..0f2d12453c 100644
--- a/pySDC/projects/Resilience/RBC.py
+++ b/pySDC/projects/Resilience/RBC.py
@@ -15,7 +15,7 @@
 import numpy as np
 
 
-def u_exact(self, t, u_init=None, t_init=None, recompute=False):
+def u_exact(self, t, u_init=None, t_init=None, recompute=False, _t0=None):
     import pickle
     import os
 
@@ -38,7 +38,10 @@ def u_exact(self, t, u_init=None, t_init=None, recompute=False):
         if t == t0:
             return u0
         else:
-            stats, _, _ = run_RBC(Tend=t, u0=u0, t0=t0, custom_description=desc)
+            u0 = u0 if _t0 is None else self.u_exact(_t0)
+            _t0 = t0 if _t0 is None else _t0
+
+            stats, _, _ = run_RBC(Tend=t, u0=u0, t0=_t0, custom_description=desc)
 
         u = get_sorted(stats, type='u', recomputed=False)[-1]
         data = u[1]
@@ -60,9 +63,9 @@ class ReachTendExactly(ConvergenceController):
 
     def setup(self, controller, params, description, **kwargs):
         defaults = {
-            "control_order": +50,
+            "control_order": +93,
             "Tend": None,
-            'min_step_size': 1e-9,
+            'min_step_size': 1e-10,
         }
         return {**defaults, **super().setup(controller, params, description, **kwargs)}
 
@@ -70,8 +73,16 @@ def get_new_step_size(self, controller, step, **kwargs):
         L = step.levels[0]
         dt = L.status.dt_new if L.status.dt_new else L.params.dt
         time_left = self.params.Tend - L.time - L.dt
-        if time_left < dt:
-            L.status.dt_new = min([dt, max([time_left, self.params.min_step_size])])
+        if time_left <= dt + self.params.min_step_size and not step.status.restart and time_left > 0:
+            L.status.dt_new = (
+                min([dt + self.params.min_step_size, max([time_left, self.params.min_step_size])])
+                + np.finfo('float').eps * 10
+            )
+            if dt != L.status.dt_new:
+                self.log(
+                    f'Reducing step size from {dt:12e} to {L.status.dt_new:.12e} because there is only {time_left:.12e} left.',
+                    step,
+                )
 
 
 def run_RBC(
@@ -182,16 +193,15 @@ def run_RBC(
 
 def generate_data_for_fault_stats():
     prob = RayleighBenard(**PROBLEM_PARAMS)
-    for t in [
-        20,
-        21,
-    ]:
-        prob.u_exact(t)
+    _ts = np.arange(22, dtype=float)
+    for i in range(len(_ts) - 1):
+        prob.u_exact(_ts[i + 1], _t0=_ts[i])
 
 
 def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recompute=False):
     from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun
     from pySDC.implementations.hooks.log_work import LogSDCIterations, LogWork
+    from pySDC.implementations.convergence_controller_classes.crash import StopAtNan
     from pySDC.helpers.stats_helper import get_sorted
     import pickle
     import os
@@ -199,6 +209,7 @@ def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recomp
     sweeper_params = {'num_nodes': num_nodes}
     level_params = {'e_tol': e_tol, 'restol': restol}
     step_params = {'maxiter': 99}
+    convergence_controllers = {StopAtNan: {}}
 
     errors = []
     dts = []
@@ -214,7 +225,12 @@ def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recomp
 
             level_params['dt'] = dts[-1]
 
-            desc = {'sweeper_params': sweeper_params, 'level_params': level_params, 'step_params': step_params}
+            desc = {
+                'sweeper_params': sweeper_params,
+                'level_params': level_params,
+                'step_params': step_params,
+                'convergence_controllers': convergence_controllers,
+            }
 
             stats, _, _ = run_RBC(
                 Tend=t + dt,
@@ -229,7 +245,10 @@ def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recomp
         e = get_sorted(stats, type='e_global_post_run')
         # k = get_sorted(stats, type='k')
 
-        errors += [e[-1][1]]
+        if len(e) > 0:
+            errors += [e[-1][1]]
+        else:
+            errors += [np.nan]
 
     errors = np.array(errors)
     dts = np.array(dts)
@@ -238,7 +257,7 @@ def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recomp
 
     errors = errors[mask]
     dts = dts[mask]
-    ax.loglog(dts, errors, label=f'{num_nodes} nodes')
+    ax.loglog(dts, errors, label=f'{num_nodes} nodes', marker='x')
     ax.loglog(
         dts, [max_error * (me / dts[0]) ** (2 * num_nodes - 1) for me in dts], ls='--', label=f'order {2*num_nodes-1}'
     )
@@ -248,16 +267,21 @@ def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recomp
 
 
 def test_order(t=14, dt=1e-1, steps=6):
+    prob = RayleighBenard(**PROBLEM_PARAMS)
+    _ts = [0, t, t + dt]
+    for i in range(len(_ts) - 1):
+        prob.u_exact(_ts[i + 1], _t0=_ts[i])
 
     import matplotlib.pyplot as plt
 
     fig, ax = plt.subplots(1, 1)
     for num_nodes in [1, 2, 3, 4]:
-        plot_order(t=t, dt=dt, steps=steps, num_nodes=num_nodes, ax=ax)
+        plot_order(t=t, dt=dt, steps=steps, num_nodes=num_nodes, ax=ax, restol=1e-9)
+    ax.set_title(f't={t:.2f}, dt={dt:.2f}')
     plt.show()
 
 
 if __name__ == '__main__':
     generate_data_for_fault_stats()
-    # test_order()
+    # test_order(t=20, dt=1., steps=7)
     # stats, _, _ = run_RBC()

From ccb9d5f92fd481383eefa6d8f7ecc33dac4b1f92 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Sun, 20 Oct 2024 19:04:21 +0200
Subject: [PATCH 10/14] Changed RBC configurations

---
 pySDC/helpers/pysdc_helper.py                 |  4 +-
 .../step_size_limiter.py                      | 64 +++++++++++++++++++
 pySDC/projects/Resilience/RBC.py              | 51 +++++++++++++--
 pySDC/projects/Resilience/fault_stats.py      |  3 +-
 pySDC/projects/Resilience/strategies.py       | 17 +++--
 pySDC/projects/Resilience/work_precision.py   | 24 ++++---
 6 files changed, 143 insertions(+), 20 deletions(-)

diff --git a/pySDC/helpers/pysdc_helper.py b/pySDC/helpers/pysdc_helper.py
index 4acd969932..666d3a4826 100644
--- a/pySDC/helpers/pysdc_helper.py
+++ b/pySDC/helpers/pysdc_helper.py
@@ -15,7 +15,7 @@ class FrozenClass(object):
 
     def __setattr__(self, key, value):
         """
-        Function called when setting arttributes
+        Function called when setting attributes
 
         Args:
             key: the attribute
@@ -35,7 +35,7 @@ def __getattr__(self, key):
         if key in self.attrs:
             return None
         else:
-            super().__getattr__(key)
+            super().__getattribute__(key)
 
     @classmethod
     def add_attr(cls, key, raise_error_if_exists=False):
diff --git a/pySDC/implementations/convergence_controller_classes/step_size_limiter.py b/pySDC/implementations/convergence_controller_classes/step_size_limiter.py
index 7551ad2cfa..c3c7ede145 100644
--- a/pySDC/implementations/convergence_controller_classes/step_size_limiter.py
+++ b/pySDC/implementations/convergence_controller_classes/step_size_limiter.py
@@ -154,3 +154,67 @@ def get_new_step_size(self, controller, S, **kwargs):
                     )
 
         return None
+
+
+class StepSizeRounding(ConvergenceController):
+    """
+    Class to round step size when using adaptive step size selection.
+    """
+
+    def setup(self, controller, params, description, **kwargs):
+        """
+        Define parameters here
+
+        Args:
+            controller (pySDC.Controller): The controller
+            params (dict): The params passed for this specific convergence controller
+            description (dict): The description object used to instantiate the controller
+
+        Returns:
+            (dict): The updated params dictionary
+        """
+        defaults = {
+            "control_order": +93,
+            "digits": 1,
+            "fac": 5,
+        }
+        return {**defaults, **super().setup(controller, params, description, **kwargs)}
+
+    @staticmethod
+    def _round_step_size(dt, fac, digits):
+        dt_rounded = None
+        exponent = np.log10(dt) // 1
+
+        dt_norm = dt / 10 ** (exponent - digits)
+        dt_norm_round = (dt_norm // fac) * fac
+        dt_rounded = dt_norm_round * 10 ** (exponent - digits)
+        return dt_rounded
+
+    def get_new_step_size(self, controller, S, **kwargs):
+        """
+        Enforce an upper and lower limit to the step size here.
+        Be aware that this is only tested when a new step size has been determined. That means if you set an initial
+        value for the step size outside of the limits, and you don't do any further step size control, that value will
+        go through.
+        Also, the final step is adjusted such that we reach Tend as best as possible, which might give step sizes below
+        the lower limit set here.
+
+        Args:
+            controller (pySDC.Controller): The controller
+            S (pySDC.Step): The current step
+
+        Returns:
+            None
+        """
+        for L in S.levels:
+            if L.status.dt_new is not None:
+                dt_rounded = self._round_step_size(L.status.dt_new, self.params.fac, self.params.digits)
+
+                if L.status.dt_new != dt_rounded:
+                    self.log(
+                        f"Step size rounded from {L.status.dt_new:.6e} to {dt_rounded:.6e}",
+                        S,
+                    )
+                    L.status.dt_new = dt_rounded
+
+        return None
diff --git a/pySDC/projects/Resilience/RBC.py b/pySDC/projects/Resilience/RBC.py
index 0f2d12453c..60aa6ee5fc 100644
--- a/pySDC/projects/Resilience/RBC.py
+++ b/pySDC/projects/Resilience/RBC.py
@@ -9,6 +9,7 @@
 from pySDC.implementations.convergence_controller_classes.estimate_extrapolation_error import (
     EstimateExtrapolationErrorNonMPI,
 )
+from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence
 
 from pySDC.core.errors import ConvergenceError
 
@@ -63,7 +64,7 @@ class ReachTendExactly(ConvergenceController):
 
     def setup(self, controller, params, description, **kwargs):
         defaults = {
-            "control_order": +93,
+            "control_order": +94,
             "Tend": None,
             'min_step_size': 1e-10,
         }
@@ -71,14 +72,19 @@ def setup(self, controller, params, description, **kwargs):
 
     def get_new_step_size(self, controller, step, **kwargs):
         L = step.levels[0]
+
+        if not CheckConvergence.check_convergence(step):
+            return None
+
         dt = L.status.dt_new if L.status.dt_new else L.params.dt
         time_left = self.params.Tend - L.time - L.dt
         if time_left <= dt + self.params.min_step_size and not step.status.restart and time_left > 0:
-            L.status.dt_new = (
+            dt_new = (
                 min([dt + self.params.min_step_size, max([time_left, self.params.min_step_size])])
                 + np.finfo('float').eps * 10
             )
-            if dt != L.status.dt_new:
+            if dt_new != L.status.dt_new:
+                L.status.dt_new = dt_new
                 self.log(
                     f'Reducing step size from {dt:12e} to {L.status.dt_new:.12e} because there is only {time_left:.12e} left.',
                     step,
@@ -135,6 +141,9 @@ def run_RBC(
 
     convergence_controllers = {}
     convergence_controllers[ReachTendExactly] = {'Tend': Tend}
+    from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeRounding
+
+    convergence_controllers[StepSizeRounding] = {}
 
     controller_params = {}
     controller_params['logger_level'] = 30
@@ -193,7 +202,7 @@ def run_RBC(
 
 def generate_data_for_fault_stats():
     prob = RayleighBenard(**PROBLEM_PARAMS)
-    _ts = np.arange(22, dtype=float)
+    _ts = np.linspace(0, 22, 220, dtype=float)
     for i in range(len(_ts) - 1):
         prob.u_exact(_ts[i + 1], _t0=_ts[i])
 
@@ -281,7 +290,41 @@ def test_order(t=14, dt=1e-1, steps=6):
     plt.show()
 
 
+def plot_step_size(t0=0, Tend=30, e_tol=1e-3, recompute=False):
+    import matplotlib.pyplot as plt
+    import pickle
+    import os
+    from pySDC.helpers.stats_helper import get_sorted
+
+    path = f'data/stats/RBC-u-{t0:.8f}-{Tend:.8f}-{e_tol:.2e}.pickle'
+    if os.path.exists(path) and not recompute:
+        with open(path, 'rb') as file:
+            stats = pickle.load(file)
+    else:
+        from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
+
+        convergence_controllers = {
+            Adaptivity: {'e_tol': e_tol, 'dt_rel_min_slope': 0.25, 'dt_min': 1e-5},
+        }
+        desc = {'convergence_controllers': convergence_controllers}
+
+        stats, _, _ = run_RBC(Tend=Tend, t0=t0, custom_description=desc)
+
+        with open(path, 'wb') as file:
+            pickle.dump(stats, file)
+
+    fig, ax = plt.subplots(1, 1)
+
+    dt = get_sorted(stats, type='dt', recomputed=False)
+    ax.plot([me[0] for me in dt], [me[1] for me in dt])
+    ax.set_ylabel(r'$\Delta t$')
+    ax.set_xlabel(r'$t$')
+    ax.set_yscale('log')
+    plt.show()
+
+
 if __name__ == '__main__':
+    # plot_step_size(0, 3)
     generate_data_for_fault_stats()
     # test_order(t=20, dt=1., steps=7)
     # stats, _, _ = run_RBC()
diff --git a/pySDC/projects/Resilience/fault_stats.py b/pySDC/projects/Resilience/fault_stats.py
index 1d0b7db909..b5f0dc9629 100644
--- a/pySDC/projects/Resilience/fault_stats.py
+++ b/pySDC/projects/Resilience/fault_stats.py
@@ -1689,7 +1689,8 @@ def main():
         stats_path='data/stats-jusuf',
         **kwargs,
     )
-    stats_analyser.run_stats_generation(runs=kwargs['runs'], step=24)
+    stats_analyser.get_HR_tol(True)
+    stats_analyser.run_stats_generation(runs=kwargs['runs'], step=12)
 
     if MPI.COMM_WORLD.rank > 0:  # make sure only one rank accesses the data
         return None
diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py
index 0dadc4a34f..7194e28f49 100644
--- a/pySDC/projects/Resilience/strategies.py
+++ b/pySDC/projects/Resilience/strategies.py
@@ -558,7 +558,7 @@ def get_custom_description(self, problem, num_procs):
             e_tol = 1e-7
             # dt_max = 0.1 * base_params['problem_params']['eps'] ** 2
         elif problem.__name__ == 'run_RBC':
-            e_tol = 1e-3
+            e_tol = 1e-4
             dt_slope_min = 0.25
 
         else:
@@ -846,7 +846,8 @@ def get_custom_description(self, problem, num_procs, *args, **kwargs):
             desc['level_params']['restol'] = 1e-11
             desc['level_params']['dt'] = 0.4 * desc['problem_params']['eps'] ** 2 / 8.0
         elif problem.__name__ == "run_RBC":
-            desc['level_params']['restol'] = 1e-4
+            desc['level_params']['dt'] = 7e-2
+            desc['level_params']['restol'] = 1e-6
         return desc
 
     def get_custom_description_for_faults(self, problem, *args, **kwargs):
@@ -1325,6 +1326,7 @@ def get_custom_description(self, problem, num_procs):
             The custom descriptions you can supply to the problem when running it
         '''
         from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK, Adaptivity
+        from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeSlopeLimiter
         from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestarting
         from pySDC.implementations.sweeper_classes.Runge_Kutta import ARK548L2SA
 
@@ -1346,6 +1348,9 @@ def get_custom_description(self, problem, num_procs):
             },
         }
 
+        if problem.__name__ == "run_RBC":
+            rk_params['convergence_controllers'][StepSizeSlopeLimiter] = {'dt_rel_min_slope': 0.25}
+
         custom_description = merge_descriptions(adaptivity_description, rk_params)
 
         return custom_description
@@ -1880,6 +1885,7 @@ def get_custom_description(self, problem, num_procs):
         '''
         from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError
         from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter
+        from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence
 
         base_params = super().get_custom_description(problem, num_procs)
         custom_description = {}
@@ -1917,12 +1923,13 @@ def get_custom_description(self, problem, num_procs):
             restol_rel = 1e-3
             # dt_max = 0.1 * base_params['problem_params']['eps'] ** 2
         elif problem.__name__ == "run_RBC":
-            e_tol = 1e-3
+            e_tol = 5e-3
             dt_slope_min = 0.25
             abort_at_growing_residual = False
-            restol_rel = 1e-4
+            restol_rel = 1e-3
             restol_max = 1e-1
-            restol_min = 1e-6
+            restol_min = 5e-7
+            self.max_slope = 4
         else:
             raise NotImplementedError(
                 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\
diff --git a/pySDC/projects/Resilience/work_precision.py b/pySDC/projects/Resilience/work_precision.py
index 0a30ce47e0..0d95e40d9a 100644
--- a/pySDC/projects/Resilience/work_precision.py
+++ b/pySDC/projects/Resilience/work_precision.py
@@ -12,6 +12,7 @@
 from pySDC.projects.Resilience.Schroedinger import run_Schroedinger
 from pySDC.projects.Resilience.quench import run_quench
 from pySDC.projects.Resilience.AC import run_AC
+from pySDC.projects.Resilience.RBC import run_RBC
 
 from pySDC.helpers.stats_helper import get_sorted, filter_stats
 from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal
@@ -38,6 +39,7 @@ def std_log(x):
     'k_linear': ('work_linear', sum, None),
     'k_Newton_no_restart': ('work_newton', sum, False),
     'k_rhs': ('work_rhs', sum, None),
+    'k_factorizations': ('work_factorizations', sum, None),
     'num_steps': ('dt', len, None),
     'restart': ('restart', sum, None),
     'dt_mean': ('dt', np.mean, False),
@@ -274,6 +276,8 @@ def record_work_precision(
                 exponents = [0, 1, 2, 3, 5][::-1]
             else:
                 exponents = [-3, -2, -1, 0, 0.2, 0.8, 1][::-1]
+        if problem.__name__ == 'run_RBC':
+            exponents = [2, 1, 0, -1]
     elif param == 'dt':
         power = 2.0
         exponents = [-1, 0, 1, 2, 3][::-1]
@@ -294,6 +298,9 @@ def record_work_precision(
             param_range = [1e-5, 1e-6, 1e-7, 1e-8, 1e-9]
         elif param == 'dt':
             param_range = [1.25, 2.5, 5.0, 10.0, 20.0][::-1]
+    if problem.__name__ == 'run_RBC':
+        if param == 'dt':
+            param_range = [7e-2, 6e-2, 5e-2, 4e-2]
 
     # run multiple times with different parameters
     for i in range(len(param_range)):
@@ -538,6 +545,7 @@ def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only
         'k_Newton': 'Newton iterations',
         'k_Newton_no_restart': 'Newton iterations (restarts excluded)',
         'k_rhs': 'right hand side evaluations',
+        'k_factorizations': 'matrix factorizations',
         't': 'wall clock time / s',
         'e_global': 'global error',
         'e_global_rel': 'relative global error',
@@ -785,14 +793,14 @@ def get_configs(mode, problem):
             AdaptivityPolynomialError,
         )
 
-        if problem.__name__ in ['run_Schroedinger', 'run_AC']:
+        if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_RBC']:
             from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
         else:
             from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
                 generic_implicit_MPI as parallel_sweeper,
             )
 
-        newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
+        newton_inexactness = False if problem.__name__ in ['run_vdp', 'run_RBC'] else True
 
         desc = {}
         desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': "EE"}
@@ -811,7 +819,7 @@ def get_configs(mode, problem):
         RK_strategies = []
         if problem.__name__ in ['run_Lorenz']:
             RK_strategies.append(ERKStrategy(useMPI=True))
-        if problem.__name__ in ['run_Schroedinger', 'run_AC']:
+        if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_RBC']:
             RK_strategies.append(ARKStrategy(useMPI=True))
         else:
             RK_strategies.append(ESDIRKStrategy(useMPI=True))
@@ -1585,23 +1593,23 @@ def aggregate_parallel_efficiency_plot():  # pragma: no cover
 
     record = comm_world.size > 1
     for mode in [
-        # 'compare_strategies',
+        'compare_strategies',
         # 'RK_comp',
         # 'parallel_efficiency',
     ]:
         params = {
             'mode': mode,
-            'runs': 5,
+            'runs': 1,
             'plotting': comm_world.rank == 0,
         }
         params_single = {
             **params,
-            'problem': run_AC,
+            'problem': run_RBC,
         }
         single_problem(**params_single, work_key='t', precision_key='e_global_rel', record=record)
 
     all_params = {
-        'record': True,
+        'record': False,
         'runs': 5,
         'work_key': 't',
         'precision_key': 'e_global_rel',
@@ -1613,7 +1621,7 @@ def aggregate_parallel_efficiency_plot():  # pragma: no cover
         'parallel_efficiency',
         'compare_strategies',
     ]:
-        all_problems(**all_params, mode=mode)
+        # all_problems(**all_params, mode=mode)
         comm_world.Barrier()
 
     if comm_world.rank == 0:

From d6c92ab2345f43783c093ded8f38ba74c0c85800 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Sun, 20 Oct 2024 19:23:57 +0200
Subject: [PATCH 11/14] A few more changes to RBC config

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

diff --git a/pySDC/projects/Resilience/fault_stats.py b/pySDC/projects/Resilience/fault_stats.py
index b5f0dc9629..6bb22c9827 100644
--- a/pySDC/projects/Resilience/fault_stats.py
+++ b/pySDC/projects/Resilience/fault_stats.py
@@ -1689,7 +1689,6 @@ def main():
         stats_path='data/stats-jusuf',
         **kwargs,
     )
-    stats_analyser.get_HR_tol(True)
     stats_analyser.run_stats_generation(runs=kwargs['runs'], step=12)
 
     if MPI.COMM_WORLD.rank > 0:  # make sure only one rank accesses the data
diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py
index 7194e28f49..95c1723618 100644
--- a/pySDC/projects/Resilience/strategies.py
+++ b/pySDC/projects/Resilience/strategies.py
@@ -559,7 +559,7 @@ def get_custom_description(self, problem, num_procs):
             # dt_max = 0.1 * base_params['problem_params']['eps'] ** 2
         elif problem.__name__ == 'run_RBC':
             e_tol = 1e-4
-            dt_slope_min = 0.25
+            dt_slope_min = 1
 
         else:
             raise NotImplementedError(
@@ -1924,7 +1924,7 @@ def get_custom_description(self, problem, num_procs):
             # dt_max = 0.1 * base_params['problem_params']['eps'] ** 2
         elif problem.__name__ == "run_RBC":
             e_tol = 5e-3
-            dt_slope_min = 0.25
+            dt_slope_min = 1.0
             abort_at_growing_residual = False
             restol_rel = 1e-3
             restol_max = 1e-1

From 63eab1072a051e90b4f32e73471a6172eebdfbc3 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Mon, 21 Oct 2024 13:58:08 +0200
Subject: [PATCH 12/14] Implemented some tests

---
 pySDC/projects/Resilience/strategies.py       |  6 +-
 pySDC/projects/Resilience/tests/test_RBC.py   | 69 +++++++++++++++++++
 pySDC/projects/Resilience/work_precision.py   | 18 ++++-
 .../test_step_size_limiter.py                 | 12 ++++
 4 files changed, 100 insertions(+), 5 deletions(-)
 create mode 100644 pySDC/projects/Resilience/tests/test_RBC.py

diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py
index 95c1723618..09bcbff62d 100644
--- a/pySDC/projects/Resilience/strategies.py
+++ b/pySDC/projects/Resilience/strategies.py
@@ -297,7 +297,7 @@ def get_base_parameters(self, problem, num_procs=1):
             'run_Schroedinger': 150,
             'run_quench': 150,
             'run_AC': 150,
-            'run_RBC': 500,
+            'run_RBC': 1000,
         }
 
         custom_description['convergence_controllers'][StopAtMaxRuntime] = {
@@ -847,7 +847,7 @@ def get_custom_description(self, problem, num_procs, *args, **kwargs):
             desc['level_params']['dt'] = 0.4 * desc['problem_params']['eps'] ** 2 / 8.0
         elif problem.__name__ == "run_RBC":
             desc['level_params']['dt'] = 7e-2
-            desc['level_params']['restol'] = 1e-6
+            desc['level_params']['restol'] = 1e-9
         return desc
 
     def get_custom_description_for_faults(self, problem, *args, **kwargs):
@@ -856,6 +856,8 @@ def get_custom_description_for_faults(self, problem, *args, **kwargs):
             desc['level_params']['dt'] = 5.0
         elif problem.__name__ == 'run_AC':
             desc['level_params']['dt'] = 0.6 * desc['problem_params']['eps'] ** 2
+        elif problem.__name__ == 'run_RBC':
+            desc['level_params']['restol'] = 1e-6
         return desc
 
     def get_reference_value(self, problem, key, op, num_procs=1):
diff --git a/pySDC/projects/Resilience/tests/test_RBC.py b/pySDC/projects/Resilience/tests/test_RBC.py
new file mode 100644
index 0000000000..71dce01bb8
--- /dev/null
+++ b/pySDC/projects/Resilience/tests/test_RBC.py
@@ -0,0 +1,69 @@
+import pytest
+
+
+@pytest.mark.mpi4py
+@pytest.mark.parametrize('Tend', [1.2345, 1, 4.0 / 3.0])
+@pytest.mark.parametrize('dt', [0.1, 1.0 / 3.0, 0.999999])
+def test_ReachTendExactly(Tend, dt, num_procs=1):
+    from pySDC.projects.Resilience.RBC import ReachTendExactly
+    from pySDC.implementations.hooks.log_solution import LogSolution
+    from pySDC.implementations.problem_classes.Van_der_Pol_implicit import vanderpol
+    from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
+    from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
+    from pySDC.helpers.stats_helper import get_sorted
+    import numpy as np
+
+    level_params = {}
+    level_params['dt'] = dt
+
+    sweeper_params = {}
+    sweeper_params['quad_type'] = 'RADAU-RIGHT'
+    sweeper_params['num_nodes'] = 2
+    sweeper_params['QI'] = 'LU'
+
+    problem_params = {
+        'mu': 0.0,
+        'newton_tol': 1e-1,
+        'newton_maxiter': 9,
+        'u0': np.array([2.0, 0.0]),
+        'relative_tolerance': True,
+    }
+
+    step_params = {}
+    step_params['maxiter'] = 2
+
+    convergence_controllers = {ReachTendExactly: {'Tend': Tend}}
+
+    controller_params = {}
+    controller_params['logger_level'] = 15
+    controller_params['hook_class'] = LogSolution
+    controller_params['mssdc_jac'] = False
+
+    description = {}
+    description['problem_class'] = vanderpol
+    description['problem_params'] = problem_params
+    description['sweeper_class'] = generic_implicit
+    description['sweeper_params'] = sweeper_params
+    description['level_params'] = level_params
+    description['step_params'] = step_params
+    description['convergence_controllers'] = convergence_controllers
+
+    t0 = 0.0
+
+    controller = controller_nonMPI(num_procs=num_procs, controller_params=controller_params, description=description)
+
+    P = controller.MS[0].levels[0].prob
+    uinit = P.u_exact(t0)
+
+    uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
+
+    u = get_sorted(stats, type='u')
+    t_last = u[-1][0]
+
+    assert np.isclose(
+        t_last, Tend, atol=1e-10
+    ), f'Expected {Tend=}, but got {t_last}, which is off by {t_last - Tend:.8e}'
+
+
+if __name__ == '__main__':
+    test_ReachTendExactly(1.2345, 0.1)
diff --git a/pySDC/projects/Resilience/work_precision.py b/pySDC/projects/Resilience/work_precision.py
index 0d95e40d9a..14b2dffe67 100644
--- a/pySDC/projects/Resilience/work_precision.py
+++ b/pySDC/projects/Resilience/work_precision.py
@@ -72,6 +72,14 @@ def get_forbidden_combinations(problem, strategy, **kwargs):
     return False
 
 
+Tends = {
+    'run_RBC': 16,
+}
+t0s = {
+    'run_RBC': 0.2,
+}
+
+
 def single_run(
     problem,
     strategy,
@@ -134,12 +142,16 @@ def single_run(
     }
     problem_args = {} if problem_args is None else problem_args
 
+    Tend = Tends.get(problem.__name__, None) if Tend is None else Tend
+    t0 = t0s.get(problem.__name__, None)
+
     stats, controller, crash = problem(
         custom_description=description,
         Tend=strategy.get_Tend(problem, num_procs) if Tend is None else Tend,
         hook_class=[LogData, LogWork, LogGlobalErrorPostRun] + hooks,
         custom_controller_params=controller_params,
         use_MPI=True,
+        t0=t0,
         comm=comm_time,
         **problem_args,
     )
@@ -277,7 +289,7 @@ def record_work_precision(
             else:
                 exponents = [-3, -2, -1, 0, 0.2, 0.8, 1][::-1]
         if problem.__name__ == 'run_RBC':
-            exponents = [2, 1, 0, -1]
+            exponents = [1, 0, -1, -2, -3, -4]
     elif param == 'dt':
         power = 2.0
         exponents = [-1, 0, 1, 2, 3][::-1]
@@ -300,7 +312,7 @@ def record_work_precision(
             param_range = [1.25, 2.5, 5.0, 10.0, 20.0][::-1]
     if problem.__name__ == 'run_RBC':
         if param == 'dt':
-            param_range = [7e-2, 6e-2, 5e-2, 4e-2]
+            param_range = [2e-1, 1e-1, 8e-2, 6e-2]
 
     # run multiple times with different parameters
     for i in range(len(param_range)):
@@ -1621,7 +1633,7 @@ def aggregate_parallel_efficiency_plot():  # pragma: no cover
         'parallel_efficiency',
         'compare_strategies',
     ]:
-        # all_problems(**all_params, mode=mode)
+        all_problems(**all_params, mode=mode)
         comm_world.Barrier()
 
     if comm_world.rank == 0:
diff --git a/pySDC/tests/test_convergence_controllers/test_step_size_limiter.py b/pySDC/tests/test_convergence_controllers/test_step_size_limiter.py
index 3c988976aa..618fbea66e 100644
--- a/pySDC/tests/test_convergence_controllers/test_step_size_limiter.py
+++ b/pySDC/tests/test_convergence_controllers/test_step_size_limiter.py
@@ -108,5 +108,17 @@ def test_step_size_limiter():
     assert L.status.dt_new == 0.5
 
 
+@pytest.mark.base
+@pytest.mark.parametrize('dt', [1 / 3, 2 / 30])
+def test_step_size_rounding(dt):
+    from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeRounding
+
+    expect = {
+        1 / 3: 0.3,
+        2 / 30: 0.065,
+    }
+    assert StepSizeRounding._round_step_size(dt, 5, 1) == expect[dt]
+
+
 if __name__ == '__main__':
     test_step_size_slope_limiter()

From 308c2fc4f9b9fae541377a9356b082e25309bc8e Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Mon, 21 Oct 2024 14:47:59 +0200
Subject: [PATCH 13/14] Restored behaviour of work_precision script

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

diff --git a/pySDC/projects/Resilience/work_precision.py b/pySDC/projects/Resilience/work_precision.py
index 14b2dffe67..86f8600ed8 100644
--- a/pySDC/projects/Resilience/work_precision.py
+++ b/pySDC/projects/Resilience/work_precision.py
@@ -1605,7 +1605,7 @@ def aggregate_parallel_efficiency_plot():  # pragma: no cover
 
     record = comm_world.size > 1
     for mode in [
-        'compare_strategies',
+        # 'compare_strategies',
         # 'RK_comp',
         # 'parallel_efficiency',
     ]:
@@ -1621,7 +1621,7 @@ def aggregate_parallel_efficiency_plot():  # pragma: no cover
         single_problem(**params_single, work_key='t', precision_key='e_global_rel', record=record)
 
     all_params = {
-        'record': False,
+        'record': True,
         'runs': 5,
         'work_key': 't',
         'precision_key': 'e_global_rel',

From 1bc2a6312dc98672634ee09b8783100330865d26 Mon Sep 17 00:00:00 2001
From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com>
Date: Mon, 21 Oct 2024 15:29:57 +0200
Subject: [PATCH 14/14] Small fixes

---
 pySDC/projects/Resilience/RBC.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/pySDC/projects/Resilience/RBC.py b/pySDC/projects/Resilience/RBC.py
index 60aa6ee5fc..46366c5d56 100644
--- a/pySDC/projects/Resilience/RBC.py
+++ b/pySDC/projects/Resilience/RBC.py
@@ -61,6 +61,10 @@ def u_exact(self, t, u_init=None, t_init=None, recompute=False, _t0=None):
 
 
 class ReachTendExactly(ConvergenceController):
+    """
+    This convergence controller will adapt the step size of (hopefully) the last step such that `Tend` is reached very closely.
+    Please pass the same `Tend` that you pass to the controller to the params for this to work.
+    """
 
     def setup(self, controller, params, description, **kwargs):
         defaults = {
@@ -275,7 +279,7 @@ def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recomp
     ax.legend(frameon=False)
 
 
-def test_order(t=14, dt=1e-1, steps=6):
+def check_order(t=14, dt=1e-1, steps=6):
     prob = RayleighBenard(**PROBLEM_PARAMS)
     _ts = [0, t, t + dt]
     for i in range(len(_ts) - 1):
@@ -326,5 +330,5 @@ def plot_step_size(t0=0, Tend=30, e_tol=1e-3, recompute=False):
 if __name__ == '__main__':
     # plot_step_size(0, 3)
     generate_data_for_fault_stats()
-    # test_order(t=20, dt=1., steps=7)
+    # check_order(t=20, dt=1., steps=7)
     # stats, _, _ = run_RBC()