Resilience experiments with Rayleigh-Benard convection (#495)

* Added basic infrastructure for Rayleigh-Benard in resilience project

* Verified order and added configurations for dt-adaptivity and base
strategy

* Implemented configurations for faults in RBC

* Added plot of RBC solution

* Limit fault targets in space index

* Fix

* Fixed fault insertion time for Hot Rod

* Changed colormap for RBC plots in the paper

* Fixed bug in reference solution

* Changed RBC configurations

* A few more changes to RBC config

* Implemented some tests

* Restored behaviour of work_precision script

* Small fixes

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/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):

pySDC/implementations/convergence_controller_classes/step_size_limiter.py

@@ -146,11 +146,75 @@ 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}",
                         S,
                     )
 
         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

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/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: