Merge pull request #284 from brownbaerchen/adaptivity_Q

Adaptivity based on error estimate based on interpolating between quadrature rules based on Thibaut

Author: pancetta

pySDC/implementations/convergence_controller_classes/adaptive_collocation.py

@@ -105,7 +105,7 @@ def switch_sweeper(self, S):
             P = L.prob
 
             # store solution of current level which will be interpolated to new level
-            u_old = L.u.copy()
+            u_old = [me.flatten() for me in L.u]
             nodes_old = L.sweep.coll.nodes.copy()
 
             # change sweeper
@@ -119,12 +119,13 @@ def switch_sweeper(self, S):
             # interpolate solution of old collocation problem to new one
             nodes_new = L.sweep.coll.nodes.copy()
             interpolator = LagrangeApproximation(points=np.append(0, nodes_old))
+
             u_inter = interpolator.getInterpolationMatrix(np.append(0, nodes_new)) @ u_old
 
             # assign the interpolated values to the nodes in the level
             for i in range(0, len(u_inter)):
                 me = P.dtype_u(P.init)
-                me[:] = u_inter[i]
+                me[:] = np.reshape(u_inter[i], P.init[0])
                 L.u[i] = me
 
             # reevaluate rhs
@@ -195,3 +196,29 @@ def post_spread_processing(self, controller, S, **kwargs):
             None
         """
         self.switch_sweeper(S)
+
+    def check_parameters(self, controller, params, description, **kwargs):
+        """
+        Check if we allow the scheme to solve the collocation problems to convergence.
+
+        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:
+            bool: Whether the parameters are compatible
+            str: The error message
+        """
+        if description["level_params"].get("restol", -1.0) <= 1e-16:
+            return (
+                False,
+                "Switching the collocation problems requires solving them to some tolerance that can be reached. Please set attainable `restol` in the level params",
+            )
+        if description["step_params"].get("maxiter", -1.0) < 99:
+            return (
+                False,
+                "Switching the collocation problems requires solving them exactly, which may require many iterations please set `maxiter` to at least 99 in the step params",
+            )
+
+        return True, ""

pySDC/implementations/convergence_controller_classes/adaptivity.py

@@ -1,5 +1,5 @@
 import numpy as np
-from pySDC.core.ConvergenceController import ConvergenceController
+from pySDC.core.ConvergenceController import ConvergenceController, Status
 from pySDC.implementations.convergence_controller_classes.step_size_limiter import (
     StepSizeLimiter,
 )
@@ -202,7 +202,7 @@ def check_parameters(self, controller, params, description, **kwargs):
             bool: Whether the parameters are compatible
             str: The error message
         """
-        if description["step_params"].get("restol", -1.0) >= 0:
+        if description["level_params"].get("restol", -1.0) >= 0:
             return (
                 False,
                 "Adaptivity needs constant order in time and hence restol in the step parameters has to be \
@@ -218,8 +218,7 @@ def check_parameters(self, controller, params, description, **kwargs):
         if "e_tol" not in params.keys():
             return (
                 False,
-                "Adaptivity needs a local tolerance! Please set some up in description['convergence_control\
-_params']['e_tol']!",
+                "Adaptivity needs a local tolerance! Please pass `e_tol` to the parameters for this convergence controller!",
             )
 
         return True, ""
@@ -389,7 +388,7 @@ def check_parameters(self, controller, params, description, **kwargs):
             bool: Whether the parameters are compatible
             str: The error message
         """
-        if description["step_params"].get("restol", -1.0) >= 0:
+        if description["level_params"].get("restol", -1.0) >= 0:
             return (
                 False,
                 "Adaptivity needs constant order in time and hence restol in the step parameters has to be \
@@ -446,3 +445,156 @@ def get_local_error_estimate(self, controller, S, **kwargs):
             float: Embedded error estimate
         """
         return S.levels[0].status.residual
+
+
+class AdaptivityCollocation(AdaptivityBase):
+    """
+    Control the step size via a collocation based estimate of the local error.
+    The error estimate works by subtracting two solutions to collocation problems with different order. You can
+    interpolate between collocation methods as much as you want but the adaptive step size selection will always be
+    based on the last switch of quadrature.
+    """
+
+    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 = {
+            "adaptive_coll_params": {},
+            "num_colls": 0,
+            **super().setup(controller, params, description, **kwargs),
+            "control_order": 220,
+        }
+
+        for key in defaults['adaptive_coll_params'].keys():
+            if type(defaults['adaptive_coll_params'][key]) == list:
+                defaults['num_colls'] = max([defaults['num_colls'], len(defaults['adaptive_coll_params'][key])])
+
+        return defaults
+
+    def setup_status_variables(self, controller, **kwargs):
+        self.status = Status(['error', 'order'])
+        self.status.error = []
+        self.status.order = []
+
+    def reset_status_variables(self, controller, **kwargs):
+        self.setup_status_variables(controller, **kwargs)
+
+    def dependencies(self, controller, description, **kwargs):
+        """
+        Load the `EstimateEmbeddedErrorCollocation` convergence controller to estimate the local error by switching
+        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.estimate_embedded_error import (
+            EstimateEmbeddedErrorCollocation,
+        )
+
+        super().dependencies(controller, description)
+
+        params = {'adaptive_coll_params': self.params.adaptive_coll_params}
+        controller.add_convergence_controller(
+            EstimateEmbeddedErrorCollocation,
+            params=params,
+            description=description,
+        )
+
+    def get_local_error_estimate(self, controller, S, **kwargs):
+        """
+        Get the collocation based embedded error estimate.
+
+        Args:
+            controller (pySDC.Controller): The controller
+            S (pySDC.Step): The current step
+
+        Returns:
+            float: Embedded error estimate
+        """
+        if len(self.status.error) > 1:
+            return self.status.error[-1][1]
+        else:
+            return 0.0
+
+    def post_iteration_processing(self, controller, step, **kwargs):
+        """
+        Get the error estimate and its order if available.
+
+        Args:
+            controller (pySDC.Controller.controller): The controller
+            step (pySDC.Step.step): The current step
+        """
+        if step.status.done:
+            lvl = step.levels[0]
+            self.status.error += [lvl.status.error_embedded_estimate_collocation]
+            self.status.order += [lvl.sweep.coll.order]
+
+    def get_new_step_size(self, controller, S, **kwargs):
+        if len(self.status.order) == self.params.num_colls:
+            lvl = S.levels[0]
+
+            # compute next step size
+            order = (
+                min(self.status.order[-2::]) + 1
+            )  # local order of less accurate of the last two collocation problems
+            e_est = self.get_local_error_estimate(controller, S)
+
+            lvl.status.dt_new = self.compute_optimal_step_size(
+                self.params.beta, lvl.params.dt, self.params.e_tol, e_est, order
+            )
+            self.log(f'Adjusting step size from {lvl.params.dt:.2e} to {lvl.status.dt_new:.2e}', S)
+
+    def check_parameters(self, controller, params, description, **kwargs):
+        """
+        Check whether parameters are compatible with whatever assumptions went into the step size functions etc.
+        For adaptivity, we need to know the order of the scheme.
+
+        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:
+            bool: Whether the parameters are compatible
+            str: The error message
+        """
+        if controller.params.mssdc_jac:
+            return (
+                False,
+                "Adaptivity needs the same order on all steps, please activate Gauss-Seidel multistep mode!",
+            )
+
+        if "e_tol" not in params.keys():
+            return (
+                False,
+                "Adaptivity needs a local tolerance! Please pass `e_tol` to the parameters for this convergence controller!",
+            )
+
+        return True, ""
+
+    def determine_restart(self, controller, S, **kwargs):
+        """
+        Check if the step wants to be restarted by comparing the estimate of the local error to a preset tolerance
+
+        Args:
+            controller (pySDC.Controller): The controller
+            S (pySDC.Step): The current step
+
+        Returns:
+            None
+        """
+        if len(self.status.order) == self.params.num_colls:
+            e_est = self.get_local_error_estimate(controller, S)
+            if e_est >= self.params.e_tol:
+                S.status.restart = True
+                self.log(f"Restarting: e={e_est:.2e} >= e_tol={self.params.e_tol:.2e}", S)

pySDC/implementations/convergence_controller_classes/estimate_embedded_error.py

@@ -270,7 +270,7 @@ def post_iteration_processing(self, controller, step, **kwargs):
             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]),
+                    max([np.finfo(float).eps, abs(self.status.u[-1] - self.status.u[-2])]),
                 )
 
     def setup_status_variables(self, controller, **kwargs):

pySDC/projects/Resilience/collocation_adaptivity.py

@@ -274,11 +274,55 @@ def order_stuff(prob):
     fig.tight_layout()
 
 
-def main():
+def adaptivity_collocation(plotting=False):
+    from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityCollocation
+
+    e_tol = 1e-7
+
+    adaptive_coll_params = {
+        'num_nodes': [2, 3],
+    }
+
+    convergence_controllers = {}
+    convergence_controllers[AdaptivityCollocation] = {'adaptive_coll_params': adaptive_coll_params, 'e_tol': e_tol}
+
+    step_params = {}
+    step_params['maxiter'] = 99
+
+    level_params = {}
+    level_params['restol'] = 1e-8
+
+    description = {}
+    description['convergence_controllers'] = convergence_controllers
+    description['step_params'] = step_params
+    description['level_params'] = level_params
+
+    controller_params = {'logger_level': 30}
+
+    stats, controller, _ = run_vdp(custom_description=description, custom_controller_params=controller_params)
+
+    e_em = get_sorted(stats, type='error_embedded_estimate_collocation', recomputed=False)
+    assert (
+        max([me[1] for me in e_em]) <= e_tol
+    ), "Exceeded threshold for local tolerance when using collocation based adaptivity"
+    assert (
+        min([me[1] for me in e_em][1:-1]) >= e_tol / 10
+    ), "Over resolved problem when using collocation based adaptivity"
+
+    if plotting:
+        from pySDC.projects.Resilience.vdp import plot_step_sizes
+
+        fig, ax = plt.subplots()
+
+        plot_step_sizes(stats, ax, 'error_embedded_estimate_collocation')
+
+
+def main(plotting=False):
+    adaptivity_collocation(plotting)
     order_stuff(run_advection)
     compare_adaptive_collocation(run_vdp)
 
 
 if __name__ == "__main__":
-    main()
+    main(True)
     plt.show()

pySDC/projects/Resilience/vdp.py

@@ -11,7 +11,7 @@
 from pySDC.projects.Resilience.hook import LogData, hook_collection
 
 
-def plot_step_sizes(stats, ax):
+def plot_step_sizes(stats, ax, e_em_key='error_embedded_estimate'):
     """
     Plot solution and step sizes to visualize the dynamics in the van der Pol equation.
 
@@ -28,7 +28,7 @@ def plot_step_sizes(stats, ax):
     p = np.array([me[1][1] for me in get_sorted(stats, type='u', recomputed=False, sortby='time')])
     t = np.array([me[0] for me in get_sorted(stats, type='u', recomputed=False, sortby='time')])
 
-    e_em = np.array(get_sorted(stats, type='error_embedded_estimate', recomputed=False, sortby='time'))[:, 1]
+    e_em = np.array(get_sorted(stats, type=e_em_key, recomputed=False, sortby='time'))[:, 1]
     dt = np.array(get_sorted(stats, type='dt', recomputed=False, sortby='time'))
     restart = np.array(get_sorted(stats, type='restart', recomputed=None, sortby='time'))
 

pySDC/tests/test_projects/test_resilience/test_adaptive_collocation.py

@@ -2,7 +2,21 @@
 
 
 @pytest.mark.base
-def test_main():
-    from pySDC.projects.Resilience.collocation_adaptivity import main
+def test_adaptivity_collocation():
+    from pySDC.projects.Resilience.collocation_adaptivity import adaptivity_collocation
 
-    main()
+    adaptivity_collocation(plotting=False)
+
+
+@pytest.mark.base
+def test_error_estimate_order():
+    from pySDC.projects.Resilience.collocation_adaptivity import order_stuff, run_advection
+
+    order_stuff(run_advection)
+
+
+@pytest.mark.base
+def test_adaptive_collocation():
+    from pySDC.projects.Resilience.collocation_adaptivity import compare_adaptive_collocation, run_vdp
+
+    compare_adaptive_collocation(run_vdp)