Verified order and added configurations for dt-adaptivity and base

strategy

Author: brownbaerchen

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

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,16 +130,14 @@ 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
     description['sweeper_class'] = imex_1st_order
     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()

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\