Fixed type and confirmed the order of the embedded estimate

Author: Thomas Baumann

pySDC/implementations/convergence_controller_classes/adaptivity.py

@@ -67,7 +67,7 @@ def get_new_step_size(self, controller, S):
         '''
         raise NotImplementedError('Please implement a rule for updating the step size!')
 
-    def compute_optimatal_step_size(self, beta, dt, e_tol, e_est, order):
+    def compute_optimal_step_size(self, beta, dt, e_tol, e_est, order):
         '''
         Compute the optimal step size for the current step based on the order of the scheme.
         This function can be called from `get_new_step_size` for various implementations of adaptivity, but notably not
@@ -192,8 +192,8 @@ def get_new_step_size(self, controller, S):
             # compute next step size
             order = S.status.iter  # embedded error estimate is same order as time marching
             e_est = self.get_local_error_estimate(controller, S)
-            L.status.dt_new = self.compute_optimatal_step_size(self.params.beta, L.params.dt, self.params.e_tol, e_est,
-                                                               order)
+            L.status.dt_new = self.compute_optimal_step_size(self.params.beta, L.params.dt, self.params.e_tol, e_est,
+                                                             order)
             self.log(f'Adjusting step size from {L.params.dt:.2e} to {L.status.dt_new:.2e}', S)
 
         return None
@@ -254,8 +254,8 @@ def get_new_step_size(self, controller, S):
             # compute next step size
             order = self.params.update_order
             e_est = self.get_local_error_estimate(controller, S)
-            L.status.dt_new = self.compute_optimatal_step_size(self.params.beta, L.params.dt, self.params.e_tol, e_est,
-                                                               order)
+            L.status.dt_new = self.compute_optimal_step_size(self.params.beta, L.params.dt, self.params.e_tol, e_est,
+                                                             order)
             self.log(f'Adjusting step size from {L.params.dt:.2e} to {L.status.dt_new:.2e}', S)
 
         return None

pySDC/implementations/convergence_controller_classes/estimate_embedded_error.py

@@ -16,7 +16,7 @@ class EstimateEmbeddedError(ConvergenceController):
 
     def setup(self, controller, params, description):
         '''
-        Add a default value for control order to the parameters
+        Add a default value for control order to the parameters and check if we are using a Runge-Kutta sweeper
 
         Args:
             controller (pySDC.Controller): The controller
@@ -26,7 +26,8 @@ def setup(self, controller, params, description):
         Returns:
             dict: Updated parameters
         '''
-        return {'control_order': -80, **params}
+        sweeper_type = 'RK' if RungeKutta in description['sweeper_class'].__bases__ else 'SDC'
+        return {'control_order': -80, 'sweeper_type': sweeper_type, **params}
 
     def dependencies(self, controller, description):
         '''
@@ -55,11 +56,15 @@ def estimate_embedded_error_serial(self, L):
         Returns:
             dtype_u: The embedded error estimate
         '''
-        sweeper = L.sweep
-        if RungeKutta in type(sweeper).__bases__:
+        if self.params.sweeper_type == 'RK':
+            # lower order solution is stored in the second to last entry of L.u
             return abs(L.u[-2] - L.u[-1])
-        else:
+        elif self.params.sweeper_type == 'SDC':
+            # order rises by one between sweeps, making this so ridiculously easy
             return abs(L.uold[-1] - L.u[-1])
+        else:
+            raise NotImplementedError(f'Don\'t know how to estimate embedded error for sweeper type \
+{self.params.sweeper_type}')
 
 
 class EstimateEmbeddedErrorNonMPI(EstimateEmbeddedError):
@@ -105,9 +110,8 @@ def post_iteration_processing(self, controller, S):
             raise NotImplementedError('Embedded error estimate only works for serial multi-level or parallel single \
 level')
 
-        if S.status.iter > 1:
+        if S.status.iter > 1 or self.params.sweeper_type == 'RK':
             for L in S.levels:
-                # order rises by one between sweeps, making this so ridiculously easy
                 temp = self.estimate_embedded_error_serial(L)
                 L.status.error_embedded_estimate = max([abs(temp - self.buffers.e_em_last), np.finfo(float).eps])
 

pySDC/playgrounds/Runge_Kutta/Runge-Kutta-Methods.ipynb

@@ -89,7 +89,7 @@
   {
    "cell_type": "code",
    "execution_count": 2,
-   "id": "90ae9863",
+   "id": "bfc3c486",
    "metadata": {},
    "outputs": [
     {
@@ -111,7 +111,7 @@
   },
   {
    "cell_type": "markdown",
-   "id": "d152c6bb",
+   "id": "a367d047",
    "metadata": {},
    "source": [
     "in the above plot you can see the regions of stability for implicit and explicit methods.\n",
@@ -125,7 +125,7 @@
   {
    "cell_type": "code",
    "execution_count": 3,
-   "id": "34d1a1c2",
+   "id": "d247424b",
    "metadata": {},
    "outputs": [
     {
@@ -147,11 +147,12 @@
   },
   {
    "cell_type": "markdown",
-   "id": "a407381b",
+   "id": "927c533f",
    "metadata": {},
    "source": [
     "We can see that the step size in black is refined by about an order of magnitude to properly resolve the fast time-scale region at around $t=5$.\n",
-    "We have to perfrom quite a few restarts, but we indeed achieve good temporal resolution throughout."
+    "We have to perfrom quite a few restarts, but we indeed achieve good temporal resolution throughout.\n",
+    "We also checked that the order of the embedded estimate is one lower than the order of the scheme, meaning it is of order 5."
    ]
   }
  ],

pySDC/projects/Resilience/accuracy_check.py

@@ -141,16 +141,20 @@ def plot(res, ax, k):
     color = plt.rcParams['axes.prop_cycle'].by_key()['color'][k - 2]
 
     for i in range(len(keys)):
+        if all([me == 0. for me in res[keys[i]]]):
+            continue
         order = get_accuracy_order(res, key=keys[i])
         if keys[i] == 'e_embedded':
             label = rf'$k={{{np.mean(order):.2f}}}$'
-            assert np.isclose(np.mean(order), k, atol=3e-1), f'Expected embedded error estimate to have order {k} \
+            assert np.isclose(np.mean(order), k, atol=4e-1), f'Expected embedded error estimate to have order {k} \
 but got {np.mean(order):.2f}'
 
         elif keys[i] == 'e_extrapolated':
             label = None
             assert np.isclose(np.mean(order), k + 1, rtol=3e-1), f'Expected extrapolation error estimate to have order \
 {k+1} but got {np.mean(order):.2f}'
+        else:
+            label = None
         ax.loglog(res['dt'], res[keys[i]], color=color, ls=ls[i], label=label)
 
     ax.set_xlabel(r'$\Delta t$')
@@ -198,11 +202,12 @@ def plot_orders(ax, ks, serial, Tend_fixed=None, custom_description=None, prob=r
         plot_order(res, ax, k)
 
 
-def plot_all_errors(ax, ks, serial, Tend_fixed=None, custom_description=None, prob=run_piline):
+def plot_all_errors(ax, ks, serial, Tend_fixed=None, custom_description=None, prob=run_piline, dt_list=None,
+                    custom_controller_params=None):
     for i in range(len(ks)):
         k = ks[i]
         res = multiple_runs(k=k, serial=serial, Tend_fixed=Tend_fixed, custom_description=custom_description,
-                            prob=prob)
+                            prob=prob, dt_list=dt_list, custom_controller_params=custom_controller_params)
 
         # visualize results
         plot(res, ax, k)

pySDC/projects/Resilience/test_Runge_Kutta_sweeper.py

@@ -1,14 +1,15 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-from pySDC.projects.Resilience.accuracy_check import plot_orders
+from pySDC.projects.Resilience.accuracy_check import plot_orders, plot_all_errors
 from pySDC.projects.Resilience.dahlquist import run_dahlquist, plot_stability
 
 from pySDC.projects.Resilience.advection import run_advection
 from pySDC.projects.Resilience.vdp import run_vdp, plot_step_sizes
 
 from pySDC.implementations.sweeper_classes.Runge_Kutta import RK1, RK4, MidpointMethod, CrankNicholson, Cash_Karp
 from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK
+from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import EstimateEmbeddedErrorNonMPI
 from pySDC.helpers.stats_helper import get_sorted
 
 
@@ -28,6 +29,7 @@ def plot_order(sweeper, prob, dt_list, description=None, ax=None, Tend_fixed=Non
     description = dict() if description is None else description
     description['sweeper_class'] = sweeper
     description['sweeper_params'] = {'implicit': True}
+    description['step_params'] = {'maxiter': 1}
 
     custom_controller_params = {'logger_level': 40}
 
@@ -63,6 +65,7 @@ def plot_stability_single(sweeper, ax=None, description=None, implicit=True, re=
     description = dict() if description is None else description
     description['sweeper_class'] = sweeper
     description['sweeper_params'] = {'implicit': implicit}
+    description['step_params'] = {'maxiter': 1}
 
     custom_controller_params = {'logger_level': 40}
 
@@ -114,6 +117,28 @@ def test_advection():
     plot_all_orders(run_advection, 1.e-3 * 2.**(-np.arange(8)), None, [RK1, MidpointMethod, CrankNicholson])
 
 
+def test_embedded_estimate_order():
+    sweeper = Cash_Karp
+    fig, ax = plt.subplots(1, 1)
+
+    # change only the things in the description that we need for adaptivity
+    convergence_controllers = dict()
+    convergence_controllers[EstimateEmbeddedErrorNonMPI] = {}
+
+    description = dict()
+    description['convergence_controllers'] = convergence_controllers
+    description['sweeper_class'] = sweeper
+    description['step_params'] = {'maxiter': 1}
+
+    custom_controller_params = {'logger_level': 40}
+
+    Tend = 7e-2
+    dt_list = Tend * 2.**(-np.arange(8))
+    prob = run_vdp
+    plot_all_errors(ax, [5], True, Tend_fixed=Tend, custom_description=description, dt_list=dt_list, prob=prob,
+                    custom_controller_params=custom_controller_params)
+
+
 def test_embedded_method():
     sweeper = Cash_Karp
     fig, ax = plt.subplots(1, 1)
@@ -129,6 +154,7 @@ def test_embedded_method():
     description = dict()
     description['convergence_controllers'] = convergence_controllers
     description['sweeper_class'] = sweeper
+    description['step_params'] = {'maxiter': 1}
 
     custom_controller_params = {'logger_level': 40}
 
@@ -145,6 +171,7 @@ def test_embedded_method():
 
 if __name__ == '__main__':
     test_embedded_method()
+    test_embedded_estimate_order()
     test_vdp()
     test_advection()
     plot_all_stability()