From c2701a0ee4b1811365051c53c5dc1d857c72d7cd Mon Sep 17 00:00:00 2001 From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com> Date: Wed, 11 Sep 2024 08:36:30 +0200 Subject: [PATCH] More efficient globally stiffly accurate DIRK (#479) * Removed "collocation update" for globally stiffly accurate RK methods * Fixed new reference values --- .../sweeper_classes/Runge_Kutta.py | 25 ++++++++++++------- pySDC/projects/Resilience/strategies.py | 4 +-- .../test_sweepers/test_Runge_Kutta_sweeper.py | 2 +- 3 files changed, 19 insertions(+), 12 deletions(-) diff --git a/pySDC/implementations/sweeper_classes/Runge_Kutta.py b/pySDC/implementations/sweeper_classes/Runge_Kutta.py index ac9f9d39a1..bc43e3dfbe 100644 --- a/pySDC/implementations/sweeper_classes/Runge_Kutta.py +++ b/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] ) diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py index e080e90cc3..2d4edaa689 100644 --- a/pySDC/projects/Resilience/strategies.py +++ b/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 diff --git a/pySDC/tests/test_sweepers/test_Runge_Kutta_sweeper.py b/pySDC/tests/test_sweepers/test_Runge_Kutta_sweeper.py index c1a3c65c26..35ed4f7150 100644 --- a/pySDC/tests/test_sweepers/test_Runge_Kutta_sweeper.py +++ b/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')