Updated resilience plots with Lorenz attractor (#497)

* Added plots for solution of Lorenz with and without faults

* Bugfix

* Fixed resilience parameters for Allen-Cahn

* Changed parameters for Lorenz

* Fixed rest of reference values for test

* Updated reference values

* Fix tests

* Fixes

* Hopefully sped up tests somewhat

* Split the tests

* Removed some lengthy tests

* Removed test that is apparently flaky based on hardware

pySDC/projects/Resilience/fault_stats.py

@@ -271,8 +271,10 @@ def generate_stats(self, strategy=None, runs=1000, reload=True, faults=True, com
                 dat['bit'][i] = faults_run[0][1][4]
                 dat['target'][i] = faults_run[0][1][5]
                 dat['rank'][i] = faults_run[0][1][6]
+
             if crash:
                 print('Code crashed!')
+                dat['error'][i] = np.inf
                 continue
 
             # record the rest of the data
@@ -876,7 +878,7 @@ def plot_things_per_things(
         args=None,
         strategies=None,
         name=None,
-        store=True,
+        store=False,
         ax=None,
         fig=None,
         plotting_args=None,
@@ -933,7 +935,7 @@ def plot_things_per_things(
         return None
 
     def plot_recovery_thresholds(
-        self, strategies=None, thresh_range=None, ax=None, mask=None, **kwargs
+        self, strategies=None, thresh_range=None, ax=None, recoverable_only=False, **kwargs
     ):  # pragma: no cover
         '''
         Plot the recovery rate for a range of thresholds
@@ -942,7 +944,6 @@ def plot_recovery_thresholds(
             strategies (list): List of the strategies you want to plot, if None, all will be plotted
             thresh_range (list): thresholds for deciding whether to accept as recovered
             ax (Matplotlib.axes): Somewhere to plot
-            mask (Numpy.ndarray of shape (n)): The mask you want to know about
 
         Returns:
             None
@@ -961,16 +962,21 @@ def plot_recovery_thresholds(
             fault_free = self.load(strategy=strategy, faults=False)
             with_faults = self.load(strategy=strategy, faults=True)
 
+            if recoverable_only:
+                recoverable_mask = self.get_fixable_faults_only(strategy)
+            else:
+                recoverable_mask = self.get_mask()
+
             for thresh_idx in range(len(thresh_range)):
                 rec_mask = self.get_mask(
                     strategy=strategy,
                     key='error',
                     val=(thresh_range[thresh_idx] * fault_free['error'].mean()),
                     op='gt',
-                    old_mask=mask,
+                    old_mask=recoverable_mask,
                 )
                 rec_rates[strategy_idx][thresh_idx] = 1.0 - len(with_faults['error'][rec_mask]) / len(
-                    with_faults['error']
+                    with_faults['error'][recoverable_mask]
                 )
 
             ax.plot(
@@ -1689,7 +1695,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=25)
 
     if MPI.COMM_WORLD.rank > 0:  # make sure only one rank accesses the data
         return None

pySDC/projects/Resilience/paper_plots.py

@@ -20,6 +20,7 @@
     DIRKStrategy,
     ERKStrategy,
     AdaptivityPolynomialError,
+    cmap,
 )
 from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal
 from pySDC.helpers.stats_helper import get_sorted
@@ -392,6 +393,102 @@ def plot_fault_vdp(bit=0):  # pragma: no cover
     savefig(fig, f'fault_bit_{bit}')
 
 
+def plot_fault_Lorenz(bit=0):  # pragma: no cover
+    """
+    Make a plot showing the impact of a fault on the Lorenz attractor without any resilience.
+    The faults are inserted in the last iteration in the last node in x such that you can best see the impact.
+
+    Args:
+        bit (int): The bit that you want to flip
+
+    Returns:
+        None
+    """
+    from pySDC.projects.Resilience.fault_stats import (
+        FaultStats,
+        BaseStrategy,
+    )
+    from pySDC.projects.Resilience.hook import LogData
+
+    stats_analyser = FaultStats(
+        prob=run_Lorenz,
+        strategies=[BaseStrategy()],
+        faults=[False, True],
+        reload=True,
+        recovery_thresh=1.1,
+        num_procs=1,
+        mode='combination',
+    )
+
+    strategy = BaseStrategy()
+
+    my_setup_mpl()
+    fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.8, 0.5))
+    colors = ['grey', strategy.color, 'magenta']
+    ls = ['--', '-']
+    markers = [None, strategy.marker]
+    do_faults = [False, True]
+    superscripts = ['*', '']
+    labels = ['x', 'x']
+
+    run = 19 + 20 * bit
+
+    for i in range(len(do_faults)):
+        stats, controller, Tend = stats_analyser.single_run(
+            strategy=BaseStrategy(),
+            run=run,
+            faults=do_faults[i],
+            hook_class=[LogData],
+        )
+        u = get_sorted(stats, type='u')
+        faults = get_sorted(stats, type='bitflip')
+        ax.plot(
+            [me[0] for me in u],
+            [me[1][0] for me in u],
+            ls=ls[i],
+            color=colors[i],
+            label=rf'${{{labels[i]}}}^{{{superscripts[i]}}}$',
+            marker=markers[i],
+            markevery=500,
+        )
+        for idx in range(len(faults)):
+            ax.axvline(faults[idx][0], color='black', label='Fault', ls=':')
+            print(
+                f'Fault at t={faults[idx][0]:.2e}, iter={faults[idx][1][1]}, node={faults[idx][1][2]}, space={faults[idx][1][3]}, bit={faults[idx][1][4]}'
+            )
+            ax.set_title(f'Fault in bit {faults[idx][1][4]}')
+
+    ax.legend(frameon=True, loc='lower left')
+    ax.set_xlabel(r'$t$')
+    savefig(fig, f'fault_bit_{bit}')
+
+
+def plot_Lorenz_solution():  # pragma: no cover
+    my_setup_mpl()
+
+    fig, axs = plt.subplots(1, 2, figsize=figsize_by_journal(JOURNAL, 1, 0.4), sharex=True)
+
+    strategy = BaseStrategy()
+    desc = strategy.get_custom_description(run_Lorenz, num_procs=1)
+    stats, controller, _ = run_Lorenz(custom_description=desc, Tend=strategy.get_Tend(run_Lorenz))
+
+    u = get_sorted(stats, recomputed=False, type='u')
+
+    axs[0].plot([me[1][0] for me in u], [me[1][2] for me in u])
+    axs[0].set_ylabel('$z$')
+    axs[0].set_xlabel('$x$')
+
+    axs[1].plot([me[1][0] for me in u], [me[1][1] for me in u])
+    axs[1].set_ylabel('$y$')
+    axs[1].set_xlabel('$x$')
+
+    for ax in axs:
+        ax.set_box_aspect(1.0)
+
+    path = 'data/paper/Lorenz_sol.pdf'
+    fig.savefig(path, bbox_inches='tight', transparent=True, dpi=200)
+
+
 def plot_quench_solution():  # pragma: no cover
     """
     Plot the solution of Quench problem over time
@@ -586,6 +683,12 @@ def work_precision():  # pragma: no cover
     all_problems(**{**all_params, 'work_key': 'param'}, mode='compare_strategies')
 
 
+def plot_recovery_rate_per_acceptance_threshold(problem):  # pragma no cover
+    stats_analyser = get_stats(problem)
+
+    stats_analyser.plot_recovery_thresholds(thresh_range=np.linspace(0.5, 1.5, 1000), recoverable_only=True)
+
+
 def make_plots_for_TIME_X_website():  # pragma: no cover
     global JOURNAL, BASE_PATH
     JOURNAL = 'JSC_beamer'
@@ -636,11 +739,14 @@ 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_Lorenz_solution()
+    plot_fault_Lorenz(0)
+    plot_fault_Lorenz(20)
     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_recovery_rate(get_stats(run_Lorenz))
+    # plot_recovery_rate_per_acceptance_threshold(run_Lorenz)
+    plt.show()
 
 
 def make_plots_for_notes():  # pragma: no cover