More efficient globally stiffly accurate DIRK (#479)

* Removed "collocation update" for globally stiffly accurate RK methods

* Fixed new reference values

Author: brownbaerchen

pySDC/implementations/sweeper_classes/Runge_Kutta.py

@@ -37,17 +37,24 @@ def __init__(self, weights, nodes, matrix):
         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_solution_stages = 1
-        self.num_nodes = matrix.shape[0] + self.num_solution_stages
+        self.globally_stiffly_accurate = np.allclose(matrix[-1], weights)
+
         self.tleft = 0.0
         self.tright = 1.0
-
-        self.nodes = np.append(np.append([0], nodes), [1])
+        self.num_solution_stages = 0 if self.globally_stiffly_accurate else 1
+        self.num_nodes = matrix.shape[0] + self.num_solution_stages
         self.weights = weights
-        self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1])
-        self.Qmat[1:-1, 1:-1] = matrix
-        self.Qmat[-1, 1:-1] = weights  # this is for computing the solution to the step from the previous stages
+
+        if self.globally_stiffly_accurate:
+            # For globally stiffly accurate methods, the last row of the Butcher tableau is the same as the weights.
+            self.nodes = np.append([0], nodes)
+            self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1])
+            self.Qmat[1:, 1:] = matrix
+        else:
+            self.nodes = np.append(np.append([0], nodes), [1])
+            self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1])
+            self.Qmat[1:-1, 1:-1] = matrix
+            self.Qmat[-1, 1:-1] = weights  # this is for computing the solution to the step from the previous stages
 
         self.left_is_node = True
         self.right_is_node = self.nodes[-1] == self.tright
@@ -277,7 +284,7 @@ def update_nodes(self):
                 rhs += lvl.dt * self.QI[m + 1, j] * self.get_full_f(lvl.f[j])
 
             # implicit solve with prefactor stemming from the diagonal of Qd, use previous stage as initial guess
-            if self.coll.implicit:
+            if self.QI[m + 1, m + 1] != 0:
                 lvl.u[m + 1][:] = prob.solve_system(
                     rhs, lvl.dt * self.QI[m + 1, m + 1], lvl.u[m], lvl.time + lvl.dt * self.coll.nodes[m + 1]
                 )

pySDC/projects/Resilience/strategies.py

@@ -1242,7 +1242,7 @@ def get_reference_value(self, problem, key, op, num_procs=1):
         """
         if problem.__name__ == "run_Lorenz":
             if key == 'work_newton' and op == sum:
-                return 1820
+                return 1456
             elif key == 'e_global_post_run' and op == max:
                 return 0.00013730538358736055
 
@@ -1433,7 +1433,7 @@ def get_reference_value(self, problem, key, op, num_procs=1):
         """
         if problem.__name__ == "run_Lorenz":
             if key == 'work_newton' and op == sum:
-                return 984
+                return 820
             elif key == 'e_global_post_run' and op == max:
                 return 3.148061889390874e-06
 

pySDC/tests/test_sweepers/test_Runge_Kutta_sweeper.py

@@ -357,5 +357,5 @@ def test_RK_sweepers_with_GPU(test_name, sweeper_name):
 
 if __name__ == '__main__':
     # test_rhs_evals('ARK54')
-    test_order('ARK548L2SAERK2')
+    test_order('CrankNicholson')
     # test_order('ARK54')