Implemented configurations for faults in RBC

Author: brownbaerchen

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()

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,32 +54,36 @@ 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):
 
     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()

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.

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,

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}')

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