Added Cash-Karp method to do embedded Runge-Kutta 5(4)

Author: Thomas Baumann

pySDC/implementations/convergence_controller_classes/estimate_embedded_error.py

@@ -3,6 +3,8 @@
 from pySDC.core.ConvergenceController import ConvergenceController, Pars
 from pySDC.implementations.convergence_controller_classes.store_uold import StoreUOld
 
+from pySDC.implementations.sweeper_classes.Runge_Kutta import RungeKutta
+
 
 class EstimateEmbeddedError(ConvergenceController):
     '''
@@ -40,6 +42,24 @@ def dependencies(self, controller, description):
         controller.add_convergence_controller(StoreUOld, description=description)
         return None
 
+    def estimate_embedded_error_serial(self, L):
+        '''
+        Estimate the serial embedded error, which may need to be modified for a parallel estimate.
+
+        Depending on the type of sweeper, the lower order solution is stored in a different place.
+
+        Args:
+            L (pySDC.level): The level
+
+        Returns:
+            dtype_u: The embedded error estimate
+        '''
+        sweeper = L.sweep
+        if RungeKutta in type(sweeper).__bases__:
+            return abs(L.u[-2] - L.u[-1])
+        else:
+            return abs(L.uold[-1] - L.u[-1])
+
 
 class EstimateEmbeddedErrorNonMPI(EstimateEmbeddedError):
 
@@ -87,7 +107,7 @@ def post_iteration_processing(self, controller, S):
         if S.status.iter > 1:
             for L in S.levels:
                 # order rises by one between sweeps, making this so ridiculously easy
-                temp = abs(L.uold[-1] - L.u[-1])
+                temp = self.estimate_embedded_error_serial(L)
                 L.status.error_embedded_estimate = max([abs(temp - self.buffers.e_em_last), np.finfo(float).eps])
 
             self.buffers.e_em_last = temp * 1.

pySDC/implementations/sweeper_classes/Runge_Kutta.py

@@ -9,7 +9,7 @@
 class ButcherTableau(object):
     def __init__(self, weights, nodes, matrix):
         """
-        Initialization routine for an collocation object
+        Initialization routine to get a quadrature matrix out of a Butcher tableau
 
         Args:
             weights (numpy.ndarray): Butcher tableau weights
@@ -61,6 +61,64 @@ def __init__(self, weights, nodes, matrix):
         self.implicit = any([matrix[i, i] != 0 for i in range(self.num_nodes - 1)])
 
 
+class ButcherTableauEmbedded(object):
+    def __init__(self, weights, nodes, matrix):
+        """
+        Initialization routine to get a quadrature matrix out of a Butcher tableau for embedded RK methods.
+
+        Be aware that the method that generates the final solution should be in the first row of the weights matrix.
+
+        Args:
+            weights (numpy.ndarray): Butcher tableau weights
+            nodes (numpy.ndarray): Butcher tableau nodes
+            matrix (numpy.ndarray): Butcher tableau entries
+        """
+        # check if the arguments have the correct form
+        if type(matrix) != np.ndarray:
+            raise ParameterError('Runge-Kutta matrix needs to be supplied as  a numpy array!')
+        elif len(np.unique(matrix.shape)) != 1 or len(matrix.shape) != 2:
+            raise ParameterError('Runge-Kutta matrix needs to be a square 2D numpy array!')
+
+        if type(weights) != np.ndarray:
+            raise ParameterError('Weights need to be supplied as a numpy array!')
+        elif len(weights.shape) != 2:
+            raise ParameterError(f'Incompatible dimension of weights! Need 2, got {len(weights.shape)}')
+        elif len(weights[0]) != matrix.shape[0]:
+            raise ParameterError(f'Incompatible number of weights! Need {matrix.shape[0]}, got {len(weights[0])}')
+
+        if type(nodes) != np.ndarray:
+            raise ParameterError('Nodes need to be supplied as a numpy array!')
+        elif len(nodes.shape) != 1:
+            raise ParameterError(f'Incompatible dimension of nodes! Need 1, got {len(nodes.shape)}')
+        elif len(nodes) != matrix.shape[0]:
+            raise ParameterError(f'Incompatible number of nodes! Need {matrix.shape[0]}, got {len(nodes)}')
+
+        # Set number of nodes, left and right interval boundaries
+        self.num_nodes = matrix.shape[0] + 2
+        self.tleft = 0.
+        self.tright = 1.
+
+        self.nodes = np.append(np.append([0], nodes), [1, 1])
+        self.weights = weights
+        self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1])
+        self.Qmat[1:-2, 1:-2] = matrix
+        self.Qmat[-1, 1:-2] = weights[0]  # this is for computing the higher order solution
+        self.Qmat[-2, 1:-2] = weights[1]  # this is for computing the lower order solution
+
+        self.left_is_node = True
+        self.right_is_node = self.nodes[-1] == self.tright
+
+        # compute distances between the nodes
+        if self.num_nodes > 1:
+            self.delta_m = self.nodes[1:] - self.nodes[:-1]
+        else:
+            self.delta_m = np.zeros(1)
+        self.delta_m[0] = self.nodes[0] - self.tleft
+
+        # check if the RK scheme is implicit
+        self.implicit = any([matrix[i, i] != 0 for i in range(self.num_nodes - 2)])
+
+
 class RungeKutta(generic_implicit):
     """
     Runge-Kutta scheme that fits the interface of a sweeper.
@@ -215,3 +273,34 @@ def __init__(self, params):
         matrix[3, 2] = 1.
         params['butcher_tableau'] = ButcherTableau(weights, nodes, matrix)
         super(RK4, self).__init__(params)
+
+
+class Heun_Euler(RungeKutta):
+    '''
+    Second order explicit embedded Runge-Kutta
+    '''
+    def __init__(self, params):
+        nodes = np.array([0, 1])
+        weights = np.array([[0.5, 0.5], [1, 0]])
+        matrix = np.zeros((2, 2))
+        matrix[1, 0] = 1
+        params['butcher_tableau'] = ButcherTableauEmbedded(weights, nodes, matrix)
+        super(Heun_Euler, self).__init__(params)
+
+
+class Cash_Karp(RungeKutta):
+    '''
+    Fifth order explicit embedded Runge-Kutta
+    '''
+    def __init__(self, params):
+        nodes = np.array([0, 0.2, 0.3, 0.6, 1., 7. / 8.])
+        weights = np.array([[37. / 378., 0., 250. / 621., 125. / 594., 0., 512. / 1771.],
+                            [2825. / 27648., 0., 18575. / 48384., 13525. / 55296., 277. / 14336., 1. / 4.]])
+        matrix = np.zeros((6, 6))
+        matrix[1, 0] = 1. / 5.
+        matrix[2, :2] = [3. / 40., 9. / 40.]
+        matrix[3, :3] = [0.3, -0.9, 1.2]
+        matrix[4, :4] = [-11. / 54., 5. / 2., -70. / 27., 35. / 27.]
+        matrix[5, :5] = [1631. / 55296., 175. / 512., 575. / 13824., 44275. / 110592., 253. / 4096.]
+        params['butcher_tableau'] = ButcherTableauEmbedded(weights, nodes, matrix)
+        super(Cash_Karp, self).__init__(params)

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

@@ -5,16 +5,19 @@
 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
+from pySDC.projects.Resilience.vdp import run_vdp, plot_step_sizes
 
-from pySDC.implementations.sweeper_classes.Runge_Kutta import RK1, RK4, MidpointMethod, CrankNicholson
+from pySDC.implementations.sweeper_classes.Runge_Kutta import RK1, RK4, MidpointMethod, CrankNicholson, Cash_Karp
+from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
+from pySDC.helpers.stats_helper import get_sorted
 
 
 colors = {
     RK1: 'blue',
     MidpointMethod: 'red',
     RK4: 'orange',
     CrankNicholson: 'purple',
+    Cash_Karp: 'teal',
 }
 
 
@@ -38,10 +41,11 @@ def plot_order(sweeper, prob, dt_list, description=None, ax=None, Tend_fixed=Non
         MidpointMethod: 3,
         RK4: 5,
         CrankNicholson: 3,
+        Cash_Karp: 6,
     }
     numerical_order = float(ax.get_lines()[-1].get_label()[7:])
     expected_order = orders.get(sweeper, numerical_order)
-    assert np.isclose(numerical_order, expected_order, atol=2.5e-1),\
+    assert np.isclose(numerical_order, expected_order, atol=2.6e-1),\
         f"Expected order {expected_order}, got {numerical_order}!"
 
     # decorate
@@ -80,11 +84,11 @@ def plot_all_stability():
     fig, axs = plt.subplots(1, 2, figsize=(11, 5))
 
     impl = [True, False]
-    sweepers = [[RK1, MidpointMethod, CrankNicholson], [RK1, MidpointMethod, RK4]]
+    sweepers = [[RK1, MidpointMethod, CrankNicholson], [RK1, MidpointMethod, RK4, Cash_Karp]]
     titles = ['implicit', 'explicit']
-    re = np.linspace(-3, 3, 400)
-    im = np.linspace(-3, 3, 400)
-    crosshair = [True, False, False]
+    re = np.linspace(-4, 4, 400)
+    im = np.linspace(-4, 4, 400)
+    crosshair = [True, False, False, False]
 
     for j in range(len(impl)):
         for i in range(len(sweepers[j])):
@@ -102,15 +106,42 @@ def plot_all_orders(prob, dt_list, Tend, sweepers):
 
 
 def test_vdp():
-    Tend = 5e-2
-    plot_all_orders(run_vdp, Tend * 2.**(-np.arange(8)), Tend, [RK1, MidpointMethod, CrankNicholson, RK4])
+    Tend = 7e-2
+    plot_all_orders(run_vdp, Tend * 2.**(-np.arange(8)), Tend, [RK1, MidpointMethod, CrankNicholson, RK4, Cash_Karp])
 
 
 def test_advection():
     plot_all_orders(run_advection, 1.e-3 * 2.**(-np.arange(8)), None, [RK1, MidpointMethod, CrankNicholson])
 
 
+def test_embedded_method():
+    sweeper = Cash_Karp
+    fig, ax = plt.subplots(1, 1)
+
+    # change only the things in the description that we need for adaptivity
+    adaptivity_params = dict()
+    adaptivity_params['e_tol'] = 1e-7
+
+    convergence_controllers = dict()
+    convergence_controllers[Adaptivity] = adaptivity_params
+
+    description = dict()
+    description['convergence_controllers'] = convergence_controllers
+    description['sweeper_class'] = sweeper
+
+    custom_controller_params = {'logger_level': 40}
+
+    stats, _, _ = run_vdp(description, 1, custom_controller_params=custom_controller_params)
+    plot_step_sizes(stats, ax)
+
+    fig.tight_layout()
+
+    dt_last = get_sorted(stats, type='dt')[-1][1]
+    assert np.isclose(dt_last, 0.11091455589374277), "Cash-Karp has computed a different last step size than before!"
+
+
 if __name__ == '__main__':
+    test_embedded_method()
     test_vdp()
     test_advection()
     plot_all_stability()

pySDC/projects/Resilience/test_Runge_Kutta_sweeper.py

@@ -1,7 +1,8 @@
 import pytest
 
-from pySDC.projects.Resilience.test_Runge_Kutta_sweeper import test_vdp, test_advection
+from pySDC.projects.Resilience.test_Runge_Kutta_sweeper import test_vdp, test_advection, test_embedded_method
 
 def test_main():
     test_vdp()
     test_advection()
+    test_embedded_method()