Merge pull request #277 from brownbaerchen/dirk

Added A-stable DIRK34 method

Author: pancetta

pySDC/implementations/sweeper_classes/Runge_Kutta.py

@@ -333,3 +333,24 @@ def __init__(self, params):
         matrix[5, :5] = [1631.0 / 55296.0, 175.0 / 512.0, 575.0 / 13824.0, 44275.0 / 110592.0, 253.0 / 4096.0]
         params['butcher_tableau'] = ButcherTableauEmbedded(weights, nodes, matrix)
         super(Cash_Karp, self).__init__(params)
+
+
+class DIRK34(RungeKutta):
+    """
+    Embedded A-stable diagonally implicit RK pair of order 3 and 4.
+
+    Taken from [here](https://doi.org/10.1007/BF01934920).
+    """
+
+    def __init__(self, params):
+        nodes = np.array([5.0 / 6.0, 10.0 / 39.0, 0, 1.0 / 6.0])
+        weights = np.array(
+            [[32.0 / 75.0, 169.0 / 300.0, 1.0 / 100.0, 0], [61.0 / 150.0, 2197.0 / 2100.0, 19.0 / 100.0, -9.0 / 14.0]]
+        )
+        matrix = np.zeros((4, 4))
+        matrix[0, 0] = 5.0 / 6.0
+        matrix[1, :2] = [-15.0 / 26.0, 5.0 / 6.0]
+        matrix[2, :3] = [215.0 / 54.0, -130.0 / 27.0, 5.0 / 6.0]
+        matrix[3, :] = [4007.0 / 6075.0, -31031.0 / 24300.0, -133.0 / 2700.0, 5.0 / 6.0]
+        params['butcher_tableau'] = ButcherTableauEmbedded(weights, nodes, matrix)
+        super().__init__(params)

pySDC/tests/test_Runge_Kutta_sweeper.py

@@ -9,6 +9,7 @@
     CrankNicholson,
     Cash_Karp,
     Heun_Euler,
+    DIRK34,
 )
 from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK
 from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import EstimateEmbeddedErrorNonMPI
@@ -131,7 +132,7 @@ def plot_stability_single(sweeper, ax=None, description=None, implicit=True, re=
     )
 
     # check if we think the method should be A-stable
-    Astable_methods = [RK1, CrankNicholson, MidpointMethod]  # only the implicit versions are A-stable
+    Astable_methods = [RK1, CrankNicholson, MidpointMethod, DIRK34]  # only the implicit versions are A-stable
     assert (
         implicit and sweeper in Astable_methods
     ) == Astable, f"Unexpected region of stability for {sweeper.__name__} sweeper!"
@@ -162,6 +163,8 @@ def test_all_stability():
             plot_stability_single(sweepers[j][i], implicit=impl[j], ax=axs[j], re=re, im=im, crosshair=crosshair[i])
         axs[j].set_title(titles[j])
 
+    plot_stability_single(DIRK34, re=re, im=im)
+
     fig.tight_layout()
 
 
@@ -216,7 +219,7 @@ def test_advection():
 
 
 @pytest.mark.base
-@pytest.mark.parametrize("sweeper", [Cash_Karp, Heun_Euler])
+@pytest.mark.parametrize("sweeper", [Cash_Karp, Heun_Euler, DIRK34])
 def test_embedded_estimate_order(sweeper):
     """
     Test the order of embedded Runge-Kutta schemes. They are not run with adaptivity here,
@@ -244,12 +247,18 @@ def test_embedded_estimate_order(sweeper):
 
     custom_controller_params = {'logger_level': 40}
 
+    expected_order = {
+        Cash_Karp: [5],
+        Heun_Euler: [2],
+        DIRK34: [4],
+    }
+
     Tend = 7e-2
     dt_list = Tend * 2.0 ** (-np.arange(8))
     prob = run_vdp
     plot_all_errors(
         ax,
-        [5] if sweeper == Cash_Karp else [2],
+        expected_order.get(sweeper, None),
         True,
         Tend_fixed=Tend,
         custom_description=description,