diff --git a/pySDC/implementations/convergence_controller_classes/estimate_embedded_error.py b/pySDC/implementations/convergence_controller_classes/estimate_embedded_error.py index 905332b9cb..545c2f6efe 100644 --- a/pySDC/implementations/convergence_controller_classes/estimate_embedded_error.py +++ b/pySDC/implementations/convergence_controller_classes/estimate_embedded_error.py @@ -93,10 +93,10 @@ def estimate_embedded_error_serial(self, L): dtype_u: The embedded error estimate """ if self.params.sweeper_type == "RK": - # lower order solution is stored in the second to last entry of L.u - return abs(L.u[-2] - L.u[-1]) + L.sweep.compute_end_point() + return abs(L.uend - L.sweep.u_secondary) elif self.params.sweeper_type == "SDC": - # order rises by one between sweeps, making this so ridiculously easy + # order rises by one between sweeps return abs(L.uold[-1] - L.u[-1]) elif self.params.sweeper_type == 'MPI': comm = L.sweep.comm diff --git a/pySDC/implementations/sweeper_classes/Runge_Kutta.py b/pySDC/implementations/sweeper_classes/Runge_Kutta.py index a0805f1f2f..5fb8dd7063 100644 --- a/pySDC/implementations/sweeper_classes/Runge_Kutta.py +++ b/pySDC/implementations/sweeper_classes/Runge_Kutta.py @@ -17,44 +17,16 @@ def __init__(self, weights, nodes, matrix): 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) != 1: - raise ParameterError(f'Incompatible dimension of weights! Need 1, got {len(weights.shape)}') - elif len(weights) != matrix.shape[0]: - raise ParameterError(f'Incompatible number of weights! Need {matrix.shape[0]}, got {len(weights)}') - - 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)}') - - self.globally_stiffly_accurate = np.allclose(matrix[-1], weights) + self.check_method(weights, nodes, matrix) self.tleft = 0.0 self.tright = 1.0 - self.num_solution_stages = 0 if self.globally_stiffly_accurate else 1 - self.num_nodes = matrix.shape[0] + self.num_solution_stages + self.num_nodes = matrix.shape[0] self.weights = weights - 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.nodes = np.append([0], nodes) + self.Qmat = np.zeros([self.num_nodes + 1, self.num_nodes + 1]) + self.Qmat[1:, 1:] = matrix self.left_is_node = True self.right_is_node = self.nodes[-1] == self.tright @@ -67,34 +39,17 @@ def __init__(self, weights, nodes, matrix): 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 - self.num_solution_stages)) + self.implicit = any(matrix[i, i] != 0 for i in range(self.num_nodes)) - -class ButcherTableauEmbedded(object): - def __init__(self, weights, nodes, matrix): + def check_method(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 that the method is entered in the correct format """ - # 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: @@ -102,31 +57,40 @@ 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 = 2 - self.num_nodes = matrix.shape[0] + self.num_solution_stages - self.tleft = 0.0 - self.tright = 1.0 + self.check_weights(weights, nodes, matrix) - 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 + def check_weights(self, weights, nodes, matrix): + """ + Check that the weights of the method are entered in the correct format + """ + if type(weights) != np.ndarray: + raise ParameterError('Weights need to be supplied as a numpy array!') + elif len(weights.shape) != 1: + raise ParameterError(f'Incompatible dimension of weights! Need 1, got {len(weights.shape)}') + elif len(weights) != matrix.shape[0]: + raise ParameterError(f'Incompatible number of weights! Need {matrix.shape[0]}, got {len(weights)}') - self.left_is_node = True - self.right_is_node = self.nodes[-1] == self.tright + @property + def globally_stiffly_accurate(self): + return np.allclose(self.Qmat[-1, 1:], self.weights) - # 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 - self.num_solution_stages)) +class ButcherTableauEmbedded(ButcherTableau): + + def check_weights(self, weights, nodes, matrix): + """ + Check that the weights of the method are entered in the correct format + """ + 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])}') + + @property + def globally_stiffly_accurate(self): + return np.allclose(self.Qmat[-1, 1:], self.weights[0]) class RungeKutta(Sweeper): @@ -292,8 +256,7 @@ def update_nodes(self): lvl.u[m + 1][:] = rhs[:] # update function values (we don't usually need to evaluate the RHS at the solution of the step) - if m < M - self.coll.num_solution_stages or self.params.eval_rhs_at_right_boundary: - lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], lvl.time + lvl.dt * self.coll.nodes[m + 1]) + lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], lvl.time + lvl.dt * self.coll.nodes[m + 1]) # indicate presence of new values at this level lvl.status.updated = True @@ -304,7 +267,28 @@ def compute_end_point(self): """ In this Runge-Kutta implementation, the solution to the step is always stored in the last node """ - self.level.uend = self.level.u[-1] + lvl = self.level + + if lvl.f[1] is None: + lvl.uend = lvl.prob.dtype_u(lvl.u[0]) + if type(self.coll) == ButcherTableauEmbedded: + self.u_secondary = lvl.prob.dtype_u(lvl.u[0]) + elif self.coll.globally_stiffly_accurate: + lvl.uend = lvl.prob.dtype_u(lvl.u[-1]) + if type(self.coll) == ButcherTableauEmbedded: + self.u_secondary = lvl.prob.dtype_u(lvl.u[0]) + for w2, k in zip(self.coll.weights[1], lvl.f[1:]): + self.u_secondary += lvl.dt * w2 * k + else: + lvl.uend = lvl.prob.dtype_u(lvl.u[0]) + if type(self.coll) == ButcherTableau: + for w, k in zip(self.coll.weights, lvl.f[1:]): + lvl.uend += lvl.dt * w * k + elif type(self.coll) == ButcherTableauEmbedded: + self.u_secondary = lvl.prob.dtype_u(lvl.u[0]) + for w1, w2, k in zip(self.coll.weights[0], self.coll.weights[1], lvl.f[1:]): + lvl.uend += lvl.dt * w1 * k + self.u_secondary += lvl.dt * w2 * k @property def level(self): @@ -356,6 +340,7 @@ class RungeKuttaIMEX(RungeKutta): """ matrix_explicit = None + weights_explicit = None ButcherTableauClass_explicit = ButcherTableau def __init__(self, params): @@ -366,6 +351,7 @@ def __init__(self, params): params: parameters for the sweeper """ super().__init__(params) + type(self).weights_explicit = self.weights if self.weights_explicit is None else self.weights_explicit self.coll_explicit = self.get_Butcher_tableau_explicit() self.QE = self.coll_explicit.Qmat @@ -388,7 +374,7 @@ def predict(self): @classmethod def get_Butcher_tableau_explicit(cls): - return cls.ButcherTableauClass_explicit(cls.weights, cls.nodes, cls.matrix_explicit) + return cls.ButcherTableauClass_explicit(cls.weights_explicit, cls.nodes, cls.matrix_explicit) def integrate(self): """ @@ -448,15 +434,47 @@ def update_nodes(self): else: lvl.u[m + 1][:] = rhs[:] - # update function values (we don't usually need to evaluate the RHS at the solution of the step) - if m < M - self.coll.num_solution_stages or self.params.eval_rhs_at_right_boundary: - lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], lvl.time + lvl.dt * self.coll.nodes[m + 1]) + # update function values + lvl.f[m + 1] = prob.eval_f(lvl.u[m + 1], lvl.time + lvl.dt * self.coll.nodes[m + 1]) # indicate presence of new values at this level lvl.status.updated = True return None + def compute_end_point(self): + """ + In this Runge-Kutta implementation, the solution to the step is always stored in the last node + """ + lvl = self.level + + if lvl.f[1] is None: + lvl.uend = lvl.prob.dtype_u(lvl.u[0]) + if type(self.coll) == ButcherTableauEmbedded: + self.u_secondary = lvl.prob.dtype_u(lvl.u[0]) + elif self.coll.globally_stiffly_accurate and self.coll_explicit.globally_stiffly_accurate: + lvl.uend = lvl.u[-1] + if type(self.coll) == ButcherTableauEmbedded: + self.u_secondary = lvl.prob.dtype_u(lvl.u[0]) + for w2, w2E, k in zip(self.coll.weights[1], self.coll_explicit.weights[1], lvl.f[1:]): + self.u_secondary += lvl.dt * (w2 * k.impl + w2E * k.expl) + else: + lvl.uend = lvl.prob.dtype_u(lvl.u[0]) + if type(self.coll) == ButcherTableau: + for w, wE, k in zip(self.coll.weights, self.coll_explicit.weights, lvl.f[1:]): + lvl.uend += lvl.dt * (w * k.impl + wE * k.expl) + elif type(self.coll) == ButcherTableauEmbedded: + self.u_secondary = lvl.u[0].copy() + for w1, w2, w1E, w2E, k in zip( + self.coll.weights[0], + self.coll.weights[1], + self.coll_explicit.weights[0], + self.coll_explicit.weights[1], + lvl.f[1:], + ): + lvl.uend += lvl.dt * (w1 * k.impl + w1E * k.expl) + self.u_secondary += lvl.dt * (w2 * k.impl + w2E * k.expl) + class ForwardEuler(RungeKutta): """ @@ -480,6 +498,14 @@ class BackwardEuler(RungeKutta): nodes, weights, matrix = generator.genCoeffs() +class IMEXEuler(RungeKuttaIMEX): + nodes = BackwardEuler.nodes + weights = BackwardEuler.weights + + matrix = BackwardEuler.matrix + matrix_explicit = ForwardEuler.matrix + + class CrankNicolson(RungeKutta): """ Implicit Runge-Kutta method of second order, A-stable. @@ -521,8 +547,13 @@ class Heun_Euler(RungeKutta): Second order explicit embedded Runge-Kutta method. """ + ButcherTableauClass = ButcherTableauEmbedded + generator = RK_SCHEMES["HEUN"]() - nodes, weights, matrix = generator.genCoeffs() + nodes, _weights, matrix = generator.genCoeffs() + weights = np.zeros((2, len(_weights))) + weights[0] = _weights + weights[1] = matrix[-1] @classmethod def get_update_order(cls): @@ -697,3 +728,41 @@ class ARK548L2SA(RungeKuttaIMEX): @classmethod def get_update_order(cls): return 5 + + +class ARK324L2SAERK(RungeKutta): + generator = RK_SCHEMES["ARK324L2SAERK"]() + nodes, weights, matrix = generator.genCoeffs(embedded=True) + ButcherTableauClass = ButcherTableauEmbedded + + @classmethod + def get_update_order(cls): + return 3 + + +class ARK324L2SAESDIRK(ARK324L2SAERK): + generator = RK_SCHEMES["ARK324L2SAESDIRK"]() + matrix = generator.Q + + +class ARK32(RungeKuttaIMEX): + ButcherTableauClass = ButcherTableauEmbedded + ButcherTableauClass_explicit = ButcherTableauEmbedded + + nodes = ARK324L2SAESDIRK.nodes + weights = ARK324L2SAESDIRK.weights + + matrix = ARK324L2SAESDIRK.matrix + matrix_explicit = ARK324L2SAERK.matrix + + @classmethod + def get_update_order(cls): + return 3 + + +class ARK2(RungeKuttaIMEX): + generator_IMP = RK_SCHEMES["ARK222EDIRK"]() + generator_EXP = RK_SCHEMES["ARK222ERK"]() + + nodes, weights, matrix = generator_IMP.genCoeffs() + _, weights_explicit, matrix_explicit = generator_EXP.genCoeffs() diff --git a/pySDC/implementations/sweeper_classes/Runge_Kutta_Nystrom.py b/pySDC/implementations/sweeper_classes/Runge_Kutta_Nystrom.py index 7837c1f348..ce3910ef02 100644 --- a/pySDC/implementations/sweeper_classes/Runge_Kutta_Nystrom.py +++ b/pySDC/implementations/sweeper_classes/Runge_Kutta_Nystrom.py @@ -4,11 +4,74 @@ from pySDC.core.sweeper import Sweeper, _Pars from pySDC.core.errors import ParameterError from pySDC.implementations.datatype_classes.particles import particles, fields, acceleration -from pySDC.implementations.sweeper_classes.Runge_Kutta import ButcherTableau -from copy import deepcopy from pySDC.implementations.sweeper_classes.Runge_Kutta import RungeKutta +class ButcherTableauNoCollUpdate(object): + """Version of Butcher Tableau that does not need a collocation update because the weights are put in the last line of Q""" + + def __init__(self, weights, nodes, matrix): + """ + Initialization routine to get a quadrature matrix out of a Butcher tableau + + 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) != 1: + raise ParameterError(f'Incompatible dimension of weights! Need 1, got {len(weights.shape)}') + elif len(weights) != matrix.shape[0]: + raise ParameterError(f'Incompatible number of weights! Need {matrix.shape[0]}, got {len(weights)}') + + 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)}') + + self.globally_stiffly_accurate = np.allclose(matrix[-1], weights) + + self.tleft = 0.0 + self.tright = 1.0 + 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 + + 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 + + # 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 - self.num_solution_stages)) + + class RungeKuttaNystrom(RungeKutta): """ Runge-Kutta scheme that fits the interface of a sweeper. @@ -34,9 +97,13 @@ class RungeKuttaNystrom(RungeKutta): All of these variables are either determined by the RK rule, or are not part of an RK scheme. Attribues: - butcher_tableau (ButcherTableau): Butcher tableau for the Runge-Kutta scheme that you want + butcher_tableau (ButcherTableauNoCollUpdate): Butcher tableau for the Runge-Kutta scheme that you want """ + ButcherTableauClass = ButcherTableauNoCollUpdate + weights_bar = None + matrix_bar = None + def __init__(self, params): """ Initialization routine for the custom sweeper @@ -44,46 +111,17 @@ def __init__(self, params): Args: params: parameters for the sweeper """ - # set up logger - self.logger = logging.getLogger('sweeper') - - essential_keys = ['butcher_tableau'] - for key in essential_keys: - if key not in params: - msg = 'need %s to instantiate step, only got %s' % (key, str(params.keys())) - self.logger.error(msg) - raise ParameterError(msg) - - # check if some parameters are set which only apply to actual sweepers - for key in ['initial_guess', 'collocation_class', 'num_nodes']: - if key in params: - self.logger.warning(f'"{key}" will be ignored by Runge-Kutta sweeper') - - # set parameters to their actual values - params['initial_guess'] = 'zero' - params['collocation_class'] = type(params['butcher_tableau']) - params['num_nodes'] = params['butcher_tableau'].num_nodes - - # disable residual computation by default - params['skip_residual_computation'] = params.get( - 'skip_residual_computation', ('IT_CHECK', 'IT_FINE', 'IT_COARSE', 'IT_UP', 'IT_DOWN') - ) - - self.params = _Pars(params) - - self.coll = params['butcher_tableau'] - self.coll_bar = params['butcher_tableau_bar'] - - # This will be set as soon as the sweeper is instantiated at the level - self.__level = None - - self.parallelizable = False - self.QI = self.coll.Qmat + super().__init__(params) + self.coll_bar = self.get_Butcher_tableau_bar() self.Qx = self.coll_bar.Qmat + @classmethod + def get_Butcher_tableau_bar(cls): + return cls.ButcherTableauClass(cls.weights_bar, cls.nodes, cls.matrix_bar) + def get_full_f(self, f): """ - Test the right hand side funtion is the correct type + Test the right hand side function is the correct type Args: f (dtype_f): Right hand side at a single node @@ -118,7 +156,7 @@ def update_nodes(self): for m in range(0, M): # build rhs, consisting of the known values from above and new values from previous nodes (at k+1) - rhs = deepcopy(L.u[0]) + rhs = P.dtype_u(L.u[0]) rhs.pos += L.dt * self.coll.nodes[m + 1] * L.u[0].vel for j in range(1, m + 1): @@ -147,7 +185,7 @@ def update_nodes(self): if self.coll.implicit: # That is why it only works for the Velocity-Verlet scheme L.f[0] = P.eval_f(L.u[0], L.time) - L.f[m + 1] = deepcopy(L.f[0]) + L.f[m + 1] = P.dtype_f(L.f[0]) else: if m != self.coll.num_nodes - 1: L.f[m + 1] = P.eval_f(L.u[m + 1], L.time + L.dt * self.coll.nodes[m]) @@ -173,23 +211,18 @@ class RKN(RungeKuttaNystrom): Chapter: II.14 Numerical methods for Second order differential equations """ - def __init__(self, params): - nodes = np.array([0.0, 0.5, 0.5, 1]) - weights = np.array([1.0, 2.0, 2.0, 1.0]) / 6.0 - matrix = np.zeros([4, 4]) - matrix[1, 0] = 0.5 - matrix[2, 1] = 0.5 - matrix[3, 2] = 1.0 - - weights_bar = np.array([1.0, 1.0, 1.0, 0]) / 6.0 - matrix_bar = np.zeros([4, 4]) - matrix_bar[1, 0] = 1 / 8 - matrix_bar[2, 0] = 1 / 8 - matrix_bar[3, 2] = 1 / 2 - params['butcher_tableau'] = ButcherTableau(weights, nodes, matrix) - params['butcher_tableau_bar'] = ButcherTableau(weights_bar, nodes, matrix_bar) + nodes = np.array([0.0, 0.5, 0.5, 1]) + weights = np.array([1.0, 2.0, 2.0, 1.0]) / 6.0 + matrix = np.zeros([4, 4]) + matrix[1, 0] = 0.5 + matrix[2, 1] = 0.5 + matrix[3, 2] = 1.0 - super(RKN, self).__init__(params) + weights_bar = np.array([1.0, 1.0, 1.0, 0]) / 6.0 + matrix_bar = np.zeros([4, 4]) + matrix_bar[1, 0] = 1 / 8 + matrix_bar[2, 0] = 1 / 8 + matrix_bar[3, 2] = 1 / 2 class Velocity_Verlet(RungeKuttaNystrom): @@ -198,15 +231,9 @@ class Velocity_Verlet(RungeKuttaNystrom): https://de.wikipedia.org/wiki/Verlet-Algorithmus """ - def __init__(self, params): - nodes = np.array([1.0, 1.0]) - weights = np.array([1 / 2, 0]) - matrix = np.zeros([2, 2]) - matrix[1, 1] = 1 - weights_bar = np.array([1 / 2, 0]) - matrix_bar = np.zeros([2, 2]) - params['butcher_tableau'] = ButcherTableau(weights, nodes, matrix) - params['butcher_tableau_bar'] = ButcherTableau(weights_bar, nodes, matrix_bar) - params['Velocity_verlet'] = True - - super(Velocity_Verlet, self).__init__(params) + nodes = np.array([1.0, 1.0]) + weights = np.array([1 / 2, 0]) + matrix = np.zeros([2, 2]) + matrix[1, 1] = 1 + weights_bar = np.array([1 / 2, 0]) + matrix_bar = np.zeros([2, 2]) diff --git a/pySDC/tests/test_sweepers/test_Runge_Kutta_sweeper.py b/pySDC/tests/test_sweepers/test_Runge_Kutta_sweeper.py index def6d4ca66..7164043dcc 100644 --- a/pySDC/tests/test_sweepers/test_Runge_Kutta_sweeper.py +++ b/pySDC/tests/test_sweepers/test_Runge_Kutta_sweeper.py @@ -18,10 +18,15 @@ 'ARK548L2SAERK', 'ARK548L2SAESDIRK2', 'ARK548L2SAERK2', + 'ARK324L2SAERK', + 'ARK324L2SAESDIRK', ] IMEX_SWEEPERS = [ 'ARK54', 'ARK548L2SA', + 'IMEXEuler', + 'ARK32', + 'ARK2', ] @@ -61,6 +66,7 @@ def single_run(sweeper_name, dt, lambdas, use_RK_sweeper=True, Tend=None, useGPU from pySDC.implementations.hooks.log_embedded_error_estimate import LogEmbeddedErrorEstimate from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import EstimateEmbeddedError from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + from pySDC.implementations.sweeper_classes.Runge_Kutta import ButcherTableauEmbedded level_params = {'dt': dt} @@ -93,7 +99,7 @@ def single_run(sweeper_name, dt, lambdas, use_RK_sweeper=True, Tend=None, useGPU 'problem_class': problem_class, 'sweeper_params': sweeper_params, 'problem_params': problem_params, - 'convergence_controllers': {EstimateEmbeddedError: {}}, + 'convergence_controllers': {}, } controller_params = { @@ -101,6 +107,12 @@ def single_run(sweeper_name, dt, lambdas, use_RK_sweeper=True, Tend=None, useGPU 'hook_class': [LogWork, LogGlobalErrorPostRun, LogSolution, LogEmbeddedErrorEstimate], } + if ( + hasattr(description['sweeper_class'], 'ButcherTableauClass') + and description['sweeper_class'].ButcherTableauClass == ButcherTableauEmbedded + ): + description['convergence_controllers'][EstimateEmbeddedError] = {} + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) if not use_RK_sweeper: @@ -108,11 +120,16 @@ def single_run(sweeper_name, dt, lambdas, use_RK_sweeper=True, Tend=None, useGPU sweeper = controller.MS[0].levels[0].sweep sweeper.QI = rk_sweeper.get_Q_matrix() sweeper.coll = rk_sweeper.get_Butcher_tableau() + _compute_end_point = type(sweeper).compute_end_point + type(sweeper).compute_end_point = rk_sweeper.compute_end_point prob = controller.MS[0].levels[0].prob ic = prob.u_exact(0) u_end, stats = controller.run(ic, 0.0, 5 * dt if Tend is None else Tend) + if not use_RK_sweeper: + type(sweeper).compute_end_point = _compute_end_point + return stats, ic, controller @@ -149,6 +166,11 @@ def test_order(sweeper_name, useGPU=False): 'ARK548L2SAESDIRK2': 6, 'ARK54': 6, 'ARK548L2SA': 6, + 'IMEXEuler': 2, + 'ARK2': 3, + 'ARK32': 4, + 'ARK324L2SAERK': 4, + 'ARK324L2SAESDIRK': 4, } dt_max = { @@ -160,6 +182,7 @@ def test_order(sweeper_name, useGPU=False): 'ARK548L2SAERK2': 1e0, 'ARK54': 5e-2, 'ARK548L2SA': 5e-2, + 'IMEXEuler': 1e-2, } lambdas = [[-1.0e-1 + 0j]] @@ -233,6 +256,8 @@ def test_stability(sweeper_name, useGPU=False): 'ARK548L2SAERK': False, 'ARK548L2SAESDIRK2': True, 'ARK548L2SAERK2': False, + 'ARK324L2SAERK': False, + 'ARK324L2SAESDIRK': True, } re = -np.logspace(-3, 2, 50) @@ -271,7 +296,7 @@ def test_rhs_evals(sweeper_name, useGPU=False): stats, _, controller = single_run(sweeper_name, 1.0, lambdas, Tend=10.0, useGPU=useGPU) sweep = controller.MS[0].levels[0].sweep - num_stages = sweep.coll.num_nodes - sweep.coll.num_solution_stages + num_stages = sweep.coll.num_nodes rhs_evaluations = [me[1] for me in get_sorted(stats, type='work_rhs')] @@ -357,5 +382,7 @@ def test_RK_sweepers_with_GPU(test_name, sweeper_name): if __name__ == '__main__': # test_rhs_evals('ARK54') - test_order('CrankNicolson') + test_order('ARK2') + # test_order('ARK54') + # test_sweeper_equivalence('Cash_Karp') # test_order('ARK54')