Merge pull request #282 from brownbaerchen/adaptivity_Q

Estimate embedded error using different collocation problems

Author: pancetta

pySDC/implementations/convergence_controller_classes/estimate_embedded_error.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from pySDC.core.ConvergenceController import ConvergenceController, Pars
+from pySDC.core.ConvergenceController import ConvergenceController, Pars, Status
 from pySDC.implementations.convergence_controller_classes.store_uold import StoreUOld
 from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate
 
@@ -203,3 +203,97 @@ def post_iteration_processing(self, controller, S, **kwargs):
                     self.send(comm, dest=S.status.slot + 1, data=temp, blocking=True)
 
         return None
+
+
+class EstimateEmbeddedErrorCollocation(ConvergenceController):
+    """
+    Estimates an embedded error based on changing the underlying quadrature rule. The error estimate is stored as
+    `error_embedded_estimate_collocation` in the status of the level. Note that we only compute the estimate on the
+    finest level. The error is stored as a tuple with the first index denoting to which iteration it belongs. This
+    is useful since the error estimate is not available immediately after, but only when the next collocation problem
+    is converged to make sure the two solutions are of different accuracy.
+
+    Changing the collocation method between iterations happens using the `AdaptiveCollocation` convergence controller.
+    Please refer to that for documentation on how to use this. Just pass the parameters for that convergence controller
+    as `adaptive_coll_params` to the parameters for this one and they will be passed on when the `AdaptiveCollocation`
+    convergence controller is automatically added while loading dependencies.
+    """
+
+    def setup(self, controller, params, description, **kwargs):
+        """
+        Add a default value for control order to the parameters
+
+        Args:
+            controller (pySDC.Controller): The controller
+            params (dict): Parameters for the convergence controller
+            description (dict): The description object used to instantiate the controller
+
+        Returns:
+            dict: Updated parameters
+        """
+        defaults = {
+            "control_order": 210,
+            "adaptive_coll_params": {},
+            **super().setup(controller, params, description, **kwargs),
+        }
+        return defaults
+
+    def dependencies(self, controller, description, **kwargs):
+        """
+        Load the `AdaptiveCollocation` convergence controller to switch between collocation problems between iterations.
+
+        Args:
+            controller (pySDC.Controller): The controller
+            description (dict): The description object used to instantiate the controller
+        """
+        from pySDC.implementations.convergence_controller_classes.adaptive_collocation import AdaptiveCollocation
+
+        controller.add_convergence_controller(
+            AdaptiveCollocation, params=self.params.adaptive_coll_params, description=description
+        )
+
+    def post_iteration_processing(self, controller, step, **kwargs):
+        """
+        Compute the embedded error as the difference between the interpolated and the current solution on the finest
+        level.
+
+        Args:
+            controller (pySDC.Controller.controller): The controller
+            step (pySDC.Step.step): The current step
+        """
+        if step.status.done:
+            lvl = step.levels[0]
+            lvl.sweep.compute_end_point()
+            self.status.u += [lvl.uend]
+            self.status.iter += [step.status.iter]
+
+            if len(self.status.u) > 1:
+                lvl.status.error_embedded_estimate_collocation = (
+                    self.status.iter[-2],
+                    abs(self.status.u[-1] - self.status.u[-2]),
+                )
+
+    def setup_status_variables(self, controller, **kwargs):
+        """
+        Add the embedded error variable to the levels and add a status variable for previous steps.
+
+        Args:
+            controller (pySDC.Controller): The controller
+        """
+        self.status = Status(['u', 'iter'])
+        self.status.u = []  # the solutions of converged collocation problems
+        self.status.iter = []  # the iteration in which the solution converged
+
+        if 'comm' in kwargs.keys():
+            steps = [controller.S]
+        else:
+            if 'active_slots' in kwargs.keys():
+                steps = [controller.MS[i] for i in kwargs['active_slots']]
+            else:
+                steps = controller.MS
+        where = ["levels", "status"]
+        for S in steps:
+            self.add_variable(S, name='error_embedded_estimate_collocation', where=where, init=None)
+
+    def reset_status_variables(self, controller, **kwargs):
+        self.setup_status_variables(controller, **kwargs)

pySDC/implementations/hooks/log_embedded_error_estimate.py

@@ -6,6 +6,26 @@ class LogEmbeddedErrorEstimate(hooks):
     Store the embedded error estimate at the end of each step as "error_embedded_estimate".
     """
 
+    def log_error(self, step, level_number, appendix=''):
+        L = step.levels[level_number]
+
+        for flavour in ['', '_collocation']:
+            if L.status.get(f'error_embedded_estimate{flavour}'):
+                if flavour == '_collocation':
+                    iter, value = L.status.error_embedded_estimate_collocation
+                else:
+                    iter = step.status.iter
+                    value = L.status.error_embedded_estimate
+                self.add_to_stats(
+                    process=step.status.slot,
+                    time=L.time + L.dt,
+                    level=L.level_index,
+                    iter=iter,
+                    sweep=L.status.sweep,
+                    type=f'error_embedded_estimate{flavour}{appendix}',
+                    value=value,
+                )
+
     def post_step(self, step, level_number, appendix=''):
         """
         Record embedded error estimate
@@ -18,18 +38,7 @@ def post_step(self, step, level_number, appendix=''):
             None
         """
         super().post_step(step, level_number)
-
-        L = step.levels[level_number]
-
-        self.add_to_stats(
-            process=step.status.slot,
-            time=L.time + L.dt,
-            level=L.level_index,
-            iter=step.status.iter,
-            sweep=L.status.sweep,
-            type=f'error_embedded_estimate{appendix}',
-            value=L.status.get('error_embedded_estimate'),
-        )
+        self.log_error(step, level_number, appendix)
 
 
 class LogEmbeddedErrorEstimatePostIter(LogEmbeddedErrorEstimate):
@@ -51,23 +60,8 @@ def post_iteration(self, step, level_number):
         Returns:
             None
         """
-        # check if the estimate is available at all
         super().post_iteration(step, level_number)
-
-        L = step.levels[level_number]
-
-        if not L.status.get('error_embedded_estimate'):
-            return None
-
-        self.add_to_stats(
-            process=step.status.slot,
-            time=L.time + L.dt,
-            level=L.level_index,
-            iter=step.status.iter - 1,
-            sweep=L.status.sweep,
-            type='error_embedded_estimate_post_iteration',
-            value=L.status.get('error_embedded_estimate'),
-        )
+        self.log_error(step, level_number, '_post_iteration')
 
     def post_step(self, step, level_number):
         super().post_step(step, level_number, appendix='_post_iteration')

pySDC/projects/Resilience/collocation_adaptivity.py

@@ -9,8 +9,12 @@
 from pySDC.projects.Resilience.hook import LogData
 from pySDC.projects.Resilience.accuracy_check import get_accuracy_order
 from pySDC.implementations.convergence_controller_classes.adaptive_collocation import AdaptiveCollocation
+from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import (
+    EstimateEmbeddedErrorCollocation,
+)
 from pySDC.core.Hooks import hooks
 from pySDC.implementations.hooks.log_errors import LogLocalErrorPostIter
+from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimatePostIter
 
 
 # define global parameters for running problems and plotting
@@ -43,11 +47,11 @@
 }
 
 special_params = {
-    'inexact': {AdaptiveCollocation: coll_params_inexact},
-    'refinement': {AdaptiveCollocation: coll_params_refinement},
-    'reduce': {AdaptiveCollocation: coll_params_reduce},
+    'inexact': {EstimateEmbeddedErrorCollocation: {'adaptive_coll_params': coll_params_inexact}},
+    'refinement': {EstimateEmbeddedErrorCollocation: {'adaptive_coll_params': coll_params_refinement}},
+    'reduce': {EstimateEmbeddedErrorCollocation: {'adaptive_coll_params': coll_params_reduce}},
     'standard': {},
-    'type': {AdaptiveCollocation: coll_params_type},
+    'type': {EstimateEmbeddedErrorCollocation: {'adaptive_coll_params': coll_params_type}},
 }
 
 
@@ -200,7 +204,7 @@ def check_order(prob, coll_name, ax, k_ax):
             Tend=dt_range[i],
             custom_description=custom_description,
             custom_controller_params=custom_controller_parameters,
-            hook_class=[LogData, LogSweeperParams, LogLocalErrorPostIter],
+            hook_class=[LogData, LogSweeperParams, LogLocalErrorPostIter, LogEmbeddedErrorEstimatePostIter],
         )
 
         sweeper_params = get_sorted(stats, type='sweeper_params', sortby='iter')
@@ -212,31 +216,43 @@ def check_order(prob, coll_name, ax, k_ax):
         e_loc = np.array([me[1] for me in get_sorted(stats, type='e_local_post_iteration', sortby='iter')])[
             converged_solution
         ]
+
+        e_em_raw = [
+            me[1] for me in get_sorted(stats, type='error_embedded_estimate_collocation_post_iteration', sortby='iter')
+        ]
+        e_em = np.array((e_em_raw + [None] if coll_name == 'refinement' else [None] + e_em_raw))
         coll_order = np.array([me[1] for me in get_sorted(stats, type='coll_order', sortby='iter')])[converged_solution]
 
-        res += [(dt_range[i], e_loc, idx[1:] - idx[:-1], labels, coll_order)]
-        # res += [(dt_range[i], np.array([me[1] for me in e_loc])[converged_solution], (idx[1:]-idx[:-1])/(idx[:-1]+1)*100, labels)]
+        res += [(dt_range[i], e_loc, idx[1:] - idx[:-1], labels, coll_order, e_em)]
 
     # assemble sth we can compute the order from
     result = {'dt': [me[0] for me in res]}
+    embedded_errors = {'dt': [me[0] for me in res]}
     num_sols = len(res[0][1])
     for i in range(num_sols):
         result[i] = [me[1][i] for me in res]
+        embedded_errors[i] = [me[5][i] for me in res]
+
         label = res[0][3][i]
         expected_order = res[0][4][i] + 1
 
-        order = get_accuracy_order(result, key=i, thresh=1e-9)
-        assert np.isclose(
-            np.mean(order), expected_order, atol=0.3
-        ), f"Expected order: {expected_order}, got {order:.2f}!"
+        ax.scatter(result['dt'], embedded_errors[i], color=CMAP[i])
+
+        for me in [result, embedded_errors]:
+            if None in me[i]:
+                continue
+            order = get_accuracy_order(me, key=i, thresh=1e-9)
+            assert np.isclose(
+                np.mean(order), expected_order, atol=0.3
+            ), f"Expected order: {expected_order}, got {order:.2f}!"
         ax.loglog(result['dt'], result[i], label=f'{label} nodes: order: {np.mean(order):.1f}', color=CMAP[i])
 
         if i > 0:
             extra_iter = [me[2][i - 1] for me in res]
             k_ax.plot(result['dt'], extra_iter, ls='--', color=CMAP[i])
     ax.legend(frameon=False)
     ax.set_xlabel(r'$\Delta t$')
-    ax.set_ylabel(r'$e_\mathrm{local}$')
+    ax.set_ylabel(r'$e_\mathrm{local}$ (lines), $e_\mathrm{embedded}$ (dots)')
     k_ax.set_ylabel(r'extra iterations')
 
 
@@ -247,6 +263,7 @@ def order_stuff(prob):
     for i in range(len(modes)):
         k_axs += [axs.flatten()[i].twinx()]
         check_order(prob, modes[i], axs.flatten()[i], k_axs[-1])
+        axs.flatten()[i].set_title(modes[i])
 
     for i in range(2):
         k_axs[i].set_ylabel('')