From f05aed239548f631ab4e41f327dc1f3e34f211fc Mon Sep 17 00:00:00 2001 From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com> Date: Mon, 20 Nov 2023 09:33:34 +0100 Subject: [PATCH] Remaining changes for current state of adaptivity paper (#378) * Remaining changes for dt-k adaptivity * Fix * Made plots a tiny bit more pretty --- .../adaptive_collocation.py | 5 - .../adaptivity.py | 90 +- .../estimate_polynomial_error.py | 19 +- pySDC/projects/Resilience/Lorenz.py | 13 +- pySDC/projects/Resilience/Schroedinger.py | 25 +- pySDC/projects/Resilience/fault_stats.py | 110 ++- pySDC/projects/Resilience/paper_plots.py | 170 ++-- pySDC/projects/Resilience/quench.py | 53 +- pySDC/projects/Resilience/strategies.py | 640 +++++++++++-- pySDC/projects/Resilience/vdp.py | 8 +- pySDC/projects/Resilience/work_precision.py | 873 +++++++++++------- .../test_resilience/test_strategies.py | 57 +- 12 files changed, 1406 insertions(+), 657 deletions(-) diff --git a/pySDC/implementations/convergence_controller_classes/adaptive_collocation.py b/pySDC/implementations/convergence_controller_classes/adaptive_collocation.py index e914273ee0..e4169170a5 100644 --- a/pySDC/implementations/convergence_controller_classes/adaptive_collocation.py +++ b/pySDC/implementations/convergence_controller_classes/adaptive_collocation.py @@ -248,10 +248,5 @@ def check_parameters(self, controller, params, description, **kwargs): False, "Switching the collocation problems requires solving them to some tolerance that can be reached. Please set attainable `restol` in the level params", ) - if description["step_params"].get("maxiter", -1.0) < 99: - return ( - False, - "Switching the collocation problems requires solving them exactly, which may require many iterations please set `maxiter` to at least 99 in the step params", - ) return True, "" diff --git a/pySDC/implementations/convergence_controller_classes/adaptivity.py b/pySDC/implementations/convergence_controller_classes/adaptivity.py index 2ce3b1ec4e..e9f74f8ec1 100644 --- a/pySDC/implementations/convergence_controller_classes/adaptivity.py +++ b/pySDC/implementations/convergence_controller_classes/adaptivity.py @@ -35,6 +35,10 @@ def setup(self, controller, params, description, **kwargs): controller.add_hook(LogStepSize) + from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence + + self.communicate_convergence = CheckConvergence.communicate_convergence + return {**defaults, **super().setup(controller, params, description, **kwargs)} def dependencies(self, controller, description, **kwargs): @@ -149,6 +153,32 @@ def determine_restart(self, controller, S, **kwargs): class AdaptivityForConvergedCollocationProblems(AdaptivityBase): + def dependencies(self, controller, description, **kwargs): + """ + Load interpolation between restarts. + + Args: + controller (pySDC.Controller): The controller + description (dict): The description object used to instantiate the controller + + Returns: + None + """ + super().dependencies(controller, description, **kwargs) + + if self.params.interpolate_between_restarts: + from pySDC.implementations.convergence_controller_classes.interpolate_between_restarts import ( + InterpolateBetweenRestarts, + ) + + controller.add_convergence_controller( + InterpolateBetweenRestarts, + description=description, + params={}, + ) + self.interpolator = controller.convergence_controllers[-1] + return None + def get_convergence(self, controller, S, **kwargs): raise NotImplementedError("Please implement a way to check if the collocation problem is converged!") @@ -167,10 +197,13 @@ def setup(self, controller, params, description, **kwargs): defaults = { 'restol_rel': None, 'e_tol_rel': None, - 'restart_at_maxiter': False, + 'restart_at_maxiter': True, 'restol_min': 1e-12, 'restol_max': 1e-5, + 'factor_if_not_converged': 4.0, + 'residual_max_tol': 1e9, 'maxiter': description['sweeper_params'].get('maxiter', 99), + 'interpolate_between_restarts': True, **super().setup(controller, params, description, **kwargs), } if defaults['restol_rel']: @@ -182,20 +215,40 @@ def setup(self, controller, params, description, **kwargs): if defaults['restart_at_maxiter']: defaults['maxiter'] = description['step_params'].get('maxiter', 99) + + self.res_last_iter = np.inf + return defaults def determine_restart(self, controller, S, **kwargs): if self.get_convergence(controller, S, **kwargs): - if self.get_local_error_estimate(controller, S, **kwargs) > self.params.e_tol: - S.status.restart = True - elif S.status.iter >= self.params.maxiter and self.params.restart_at_maxiter: + self.res_last_iter = np.inf + + if self.params.restart_at_maxiter and S.levels[0].status.residual > S.levels[0].params.restol: + self.trigger_restart_upon_nonconvergence(S) + elif self.get_local_error_estimate(controller, S, **kwargs) > self.params.e_tol: S.status.restart = True - for L in S.levels: - L.status.dt_new = L.params.dt / 2.0 - self.log( - f'Collocation problem not converged after max. number of iterations, halving step size to {L.status.dt_new:.2e}', - S, - ) + elif S.status.time_size == 1 and self.res_last_iter < S.levels[0].status.residual and S.status.iter > 0: + self.trigger_restart_upon_nonconvergence(S) + elif S.levels[0].status.residual > self.params.residual_max_tol: + self.trigger_restart_upon_nonconvergence(S) + + if self.params.useMPI: + self.communicate_convergence(self, controller, S, **kwargs) + + self.res_last_iter = S.levels[0].status.residual * 1.0 + + def trigger_restart_upon_nonconvergence(self, S): + S.status.restart = True + S.status.force_done = True + for L in S.levels: + L.status.dt_new = L.params.dt / self.params.factor_if_not_converged + self.log( + f'Collocation problem not converged. Reducing step size to {L.status.dt_new:.2e}', + S, + ) + if self.params.interpolate_between_restarts: + self.interpolator.status.skip_interpolation = True class Adaptivity(AdaptivityBase): @@ -251,7 +304,7 @@ def dependencies(self, controller, description, **kwargs): """ from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import EstimateEmbeddedError - super().dependencies(controller, description) + super().dependencies(controller, description, **kwargs) controller.add_convergence_controller( EstimateEmbeddedError.get_implementation(self.params.embedded_error_flavor, self.params.useMPI), @@ -398,6 +451,7 @@ def setup(self, controller, params, description, **kwargs): - control_order (int): The order relative to other convergence controllers - e_tol_low (float): Lower absolute threshold for the residual - e_tol (float): Upper absolute threshold for the residual + - use_restol (bool): Restart if the residual tolerance was not reached - max_restarts: Override maximum number of restarts Args: @@ -412,6 +466,7 @@ def setup(self, controller, params, description, **kwargs): "control_order": -45, "e_tol_low": 0, "e_tol": np.inf, + "use_restol": False, "max_restarts": 99 if "e_tol_low" in params else None, "allowed_modifications": ['increase', 'decrease'], # what we are allowed to do with the step size } @@ -481,7 +536,9 @@ def get_new_step_size(self, controller, S, **kwargs): dt_planned = L.status.dt_new if L.status.dt_new is not None else L.params.dt - if res > self.params.e_tol and 'decrease' in self.params.allowed_modifications: + if ( + res > self.params.e_tol or (res > L.params.restol and self.params.use_restol) + ) and 'decrease' in self.params.allowed_modifications: L.status.dt_new = min([dt_planned, L.params.dt / 2.0]) self.log(f'Adjusting step size from {L.params.dt:.2e} to {L.status.dt_new:.2e}', S) elif res < self.params.e_tol_low and 'increase' in self.params.allowed_modifications: @@ -561,7 +618,7 @@ def dependencies(self, controller, description, **kwargs): EstimateEmbeddedErrorCollocation, ) - super().dependencies(controller, description) + super().dependencies(controller, description, **kwargs) params = {'adaptive_coll_params': self.params.adaptive_coll_params} controller.add_convergence_controller( @@ -694,7 +751,7 @@ def dependencies(self, controller, description, **kwargs): EstimateExtrapolationErrorWithinQ, ) - super().dependencies(controller, description) + super().dependencies(controller, description, **kwargs) controller.add_convergence_controller( EstimateExtrapolationErrorWithinQ, @@ -761,11 +818,12 @@ def setup(self, controller, params, description, **kwargs): defaults = { 'control_order': -50, + **super().setup(controller, params, description, **kwargs), **params, } self.check_convergence = CheckConvergence.check_convergence - return {**defaults, **super().setup(controller, params, description, **kwargs)} + return defaults def get_convergence(self, controller, S, **kwargs): return self.check_convergence(S) @@ -785,7 +843,7 @@ def dependencies(self, controller, description, **kwargs): EstimatePolynomialError, ) - super().dependencies(controller, description) + super().dependencies(controller, description, **kwargs) controller.add_convergence_controller( EstimatePolynomialError, diff --git a/pySDC/implementations/convergence_controller_classes/estimate_polynomial_error.py b/pySDC/implementations/convergence_controller_classes/estimate_polynomial_error.py index 996c40b2b3..7eb9989d96 100644 --- a/pySDC/implementations/convergence_controller_classes/estimate_polynomial_error.py +++ b/pySDC/implementations/convergence_controller_classes/estimate_polynomial_error.py @@ -58,6 +58,8 @@ def setup(self, controller, params, description, **kwargs): 'You cannot interpolate with lower accuracy to the end point if the end point is a node!' ) + self.interpolation_matrix = None + return defaults def reset_status_variables(self, controller, **kwargs): @@ -131,16 +133,18 @@ def post_iteration_processing(self, controller, S, **kwargs): nodes = np.append(np.append(0, coll.nodes), 1.0) estimate_on_node = self.params.estimate_on_node - interpolator = LagrangeApproximation( - points=[nodes[i] for i in range(coll.num_nodes + 1) if i != estimate_on_node] - ) - interpolation_matrix = interpolator.getInterpolationMatrix([nodes[estimate_on_node]]) + if self.interpolation_matrix is None: + interpolator = LagrangeApproximation( + points=[nodes[i] for i in range(coll.num_nodes + 1) if i != estimate_on_node] + ) + self.interpolation_matrix = interpolator.getInterpolationMatrix([nodes[estimate_on_node]]) + u = [ L.u[i].flatten() if L.u[i] is not None else L.u[i] for i in range(coll.num_nodes + 1) if i != estimate_on_node ] - u_inter = self.matmul(interpolation_matrix, u)[0].reshape(L.prob.init[0]) + u_inter = self.matmul(self.interpolation_matrix, u)[0].reshape(L.prob.init[0]) # compute end point if needed if estimate_on_node == len(nodes) - 1: @@ -161,6 +165,11 @@ def post_iteration_processing(self, controller, S, **kwargs): else: L.status.error_embedded_estimate = abs(u_inter - high_order_sol) + self.debug( + f'Obtained error estimate: {L.status.error_embedded_estimate:.2e} of order {L.status.order_embedded_estimate}', + S, + ) + def check_parameters(self, controller, params, description, **kwargs): """ Check if we allow the scheme to solve the collocation problems to convergence. diff --git a/pySDC/projects/Resilience/Lorenz.py b/pySDC/projects/Resilience/Lorenz.py index 26d79be21f..1f69eee662 100644 --- a/pySDC/projects/Resilience/Lorenz.py +++ b/pySDC/projects/Resilience/Lorenz.py @@ -4,12 +4,12 @@ from pySDC.helpers.stats_helper import get_sorted from pySDC.implementations.problem_classes.Lorenz import LorenzAttractor -from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity from pySDC.core.Errors import ConvergenceError from pySDC.projects.Resilience.hook import LogData, hook_collection from pySDC.projects.Resilience.strategies import merge_descriptions +from pySDC.projects.Resilience.sweepers import generic_implicit_efficient def run_Lorenz( @@ -37,7 +37,7 @@ def run_Lorenz( Returns: dict: The stats object controller: The controller - Tend: The time that was supposed to be integrated to + bool: Whether the code crashed """ # initialize level parameters @@ -72,7 +72,7 @@ def run_Lorenz( description = dict() description['problem_class'] = LorenzAttractor description['problem_params'] = problem_params - description['sweeper_class'] = generic_implicit + description['sweeper_class'] = generic_implicit_efficient description['sweeper_params'] = sweeper_params description['level_params'] = level_params description['step_params'] = step_params @@ -105,12 +105,15 @@ def run_Lorenz( prepare_controller_for_faults(controller, fault_stuff) # call main function to get things done... + crash = False try: uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - except ConvergenceError: + except (ConvergenceError, ZeroDivisionError) as e: + print(f'Warning: Premature termination: {e}') stats = controller.return_stats() + crash = True - return stats, controller, Tend + return stats, controller, crash def plot_solution(stats): # pragma: no cover diff --git a/pySDC/projects/Resilience/Schroedinger.py b/pySDC/projects/Resilience/Schroedinger.py index 0c88b8aa92..56a593521a 100644 --- a/pySDC/projects/Resilience/Schroedinger.py +++ b/pySDC/projects/Resilience/Schroedinger.py @@ -1,18 +1,14 @@ -from mpi4py import MPI import numpy as np -from pySDC.helpers.stats_helper import get_sorted - from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI -from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order -from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit from pySDC.implementations.problem_classes.NonlinearSchroedinger_MPIFFT import ( nonlinearschroedinger_imex, nonlinearschroedinger_fully_implicit, ) -from pySDC.implementations.transfer_classes.TransferMesh_MPIFFT import fft_to_fft from pySDC.projects.Resilience.hook import LogData, hook_collection from pySDC.projects.Resilience.strategies import merge_descriptions +from pySDC.projects.Resilience.sweepers import imex_1st_order_efficient, generic_implicit_efficient +from pySDC.core.Errors import ConvergenceError from pySDC.core.Hooks import hooks @@ -95,7 +91,7 @@ def run_Schroedinger( Returns: dict: The stats object controller: The controller - Tend: The time that was supposed to be integrated to + bool: If the code crashed """ if custom_description is not None: problem_params = custom_description.get('problem_params', {}) @@ -119,6 +115,7 @@ def run_Schroedinger( sweeper_params['quad_type'] = 'RADAU-RIGHT' sweeper_params['num_nodes'] = 3 sweeper_params['QI'] = 'IE' + sweeper_params['QE'] = 'PIC' sweeper_params['initial_guess'] = 'spread' # initialize problem parameters @@ -148,7 +145,7 @@ def run_Schroedinger( description = dict() description['problem_params'] = problem_params description['problem_class'] = nonlinearschroedinger_imex if imex else nonlinearschroedinger_fully_implicit - description['sweeper_class'] = imex_1st_order if imex else generic_implicit + description['sweeper_class'] = imex_1st_order_efficient if imex else generic_implicit_efficient description['sweeper_params'] = sweeper_params description['level_params'] = level_params description['step_params'] = step_params @@ -187,12 +184,20 @@ def run_Schroedinger( prepare_controller_for_faults(controller, fault_stuff, rnd_args) # call main function to get things done... - uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + crash = False + try: + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + except (ConvergenceError, OverflowError) as e: + print(f'Warning: Premature termination: {e}') + stats = controller.return_stats() + crash = True - return stats, controller, Tend + return stats, controller, crash def main(): + from mpi4py import MPI + stats, _, _ = run_Schroedinger(space_comm=MPI.COMM_WORLD, hook_class=live_plotting, imex=False) plt.show() diff --git a/pySDC/projects/Resilience/fault_stats.py b/pySDC/projects/Resilience/fault_stats.py index a7411b3d7b..4558002ea1 100644 --- a/pySDC/projects/Resilience/fault_stats.py +++ b/pySDC/projects/Resilience/fault_stats.py @@ -29,7 +29,7 @@ LOGGER_LEVEL = 40 RECOVERY_THRESH_ABS = { - run_quench: 5e-3, + # run_quench: 5e-3, # run_Schroedinger: 1.7e-6, } @@ -142,7 +142,7 @@ def run_stats_generation( return None comm = MPI.COMM_WORLD if comm is None else comm - step = (runs if step is None else step) if comm.size == 1 else comm.size + step = (runs if comm.size == 1 else comm.size) if step is None else step _runs_partial = step if _runs_partial == 0 else _runs_partial reload = self.reload or _reload @@ -175,7 +175,7 @@ def run_stats_generation( for j in range(0, len(strategies), comm.size): self.generate_stats( - strategy=strategies[j + (comm.rank % len(strategies) % (len(strategies)) - j)], + strategy=strategies[j + (comm.rank % len(strategies) % (len(strategies) - j))], runs=min(_runs_partial, max_runs), faults=faults, reload=reload, @@ -255,28 +255,30 @@ def generate_stats(self, strategy=None, runs=1000, reload=True, faults=True, com continue # perform a single experiment with the correct random seed - stats, controller, Tend = self.single_run(strategy=strategy, run=i, faults=faults, space_comm=space_comm) + stats, controller, crash = self.single_run(strategy=strategy, run=i, faults=faults, space_comm=space_comm) # get the data from the stats faults_run = get_sorted(stats, type='bitflip') - t, u = get_sorted(stats, type='u', recomputed=False)[-1] - error = self.get_error(u, t, controller, strategy, Tend) + if faults: + assert len(faults_run) > 0, f"Did not record a fault in run {i} of {strategy.name}!" + dat['level'][i] = faults_run[0][1][0] + dat['iteration'][i] = faults_run[0][1][1] + dat['node'][i] = faults_run[0][1][2] + dat['problem_pos'] += [faults_run[0][1][3]] + dat['bit'][i] = faults_run[0][1][4] + dat['target'][i] = faults_run[0][1][5] + dat['rank'][i] = faults_run[0][1][6] + if crash: + print('Code crashed!') + continue + + # record the rest of the data + t, u = get_sorted(stats, type='u', recomputed=False)[-1] + error = self.get_error(u, t, controller, strategy, self.get_Tend()) total_iteration = sum([k[1] for k in get_sorted(stats, type='k')]) total_newton_iteration = sum([k[1] for k in get_sorted(stats, type='work_newton')]) - # record the new data point - if faults: - if len(faults_run) > 0: - dat['level'][i] = faults_run[0][1][0] - dat['iteration'][i] = faults_run[0][1][1] - dat['node'][i] = faults_run[0][1][2] - dat['problem_pos'] += [faults_run[0][1][3]] - dat['bit'][i] = faults_run[0][1][4] - dat['target'][i] = faults_run[0][1][5] - dat['rank'][i] = faults_run[0][1][6] - else: - assert self.mode == 'regular', f'No faults where recorded in run {i} of strategy {strategy.name}!' dat['error'][i] = error dat['total_iteration'][i] = total_iteration dat['total_newton_iteration'][i] = total_newton_iteration @@ -320,7 +322,8 @@ def get_error(self, u, t, controller, strategy, Tend): Tend_reached = t + lvl.params.dt_initial >= Tend if Tend_reached: - return abs(u - lvl.prob.u_exact(t=t)) + u_star = lvl.prob.u_exact(t=t) + return abs(u - u_star) / abs(u_star) else: print(f'Final time was not reached! Code crashed at t={t:.2f} instead of reaching Tend={Tend:.2f}') return np.inf @@ -357,7 +360,7 @@ def single_run( force_params = {} if force_params is None else force_params # build the custom description - custom_description = strategy.get_custom_description(self.prob, self.num_procs) + custom_description = strategy.get_custom_description_for_faults(self.prob, self.num_procs) for k in force_params.keys(): custom_description[k] = {**custom_description.get(k, {}), **force_params[k]} @@ -515,7 +518,7 @@ def scrutinize(self, strategy, logger_level=15, **kwargs): u_all = get_sorted(stats, type='u', recomputed=False) t, u = get_sorted(stats, type='u', recomputed=False)[np.argmax([me[0] for me in u_all])] - k = get_sorted(stats, type='k') + k = [me[1] for me in get_sorted(stats, type='k')] fault_hook = [me for me in controller.hooks if type(me) == FaultInjector] @@ -546,8 +549,7 @@ def scrutinize(self, strategy, logger_level=15, **kwargs): recovery_thresh = self.get_thresh(strategy) print( - f'e={error:.2e}, e^*={e_star:.2e}, thresh: {recovery_thresh:.2e} -> recovered: \ -{error < recovery_thresh}' + f'e={error:.2e}, e^*={e_star:.2e}, thresh: {recovery_thresh:.2e} -> recovered: {error < recovery_thresh}, e_rel = {error/abs(u):.2e}' ) print(f'k: sum: {np.sum(k)}, min: {np.min(k)}, max: {np.max(k)}, mean: {np.mean(k):.2f},') @@ -560,7 +562,7 @@ def scrutinize(self, strategy, logger_level=15, **kwargs): # checkout the step size dt = [me[1] for me in get_sorted(stats, type='dt')] - print(f'dt: min: {np.min(dt):.2e}, max: {np.max(dt):.2e}, mean: {np.mean(dt):.2e}') + print(f't: {t:.2f}, dt: min: {np.min(dt):.2e}, max: {np.max(dt):.2e}, mean: {np.mean(dt):.2e}') # restarts restarts = [me[1] for me in get_sorted(stats, type='restart')] @@ -789,7 +791,15 @@ def extra_mean(self, dat, no_faults, thingA, mask): return None def plot_thingA_per_thingB( - self, strategy, thingA, thingB, ax=None, mask=None, recovered=False, op=None + self, + strategy, + thingA, + thingB, + ax=None, + mask=None, + recovered=False, + op=None, + **kwargs, ): # pragma: no cover ''' Plot thingA vs. thingB for a single strategy @@ -832,10 +842,17 @@ def plot_thingA_per_thingB( marker=strategy.marker, ls='--', linewidth=3, + **kwargs, ) ax.plot( - admissable_thingB, me, label=f'{strategy.label}', color=strategy.color, marker=strategy.marker, linewidth=2 + admissable_thingB, + me, + label=f'{strategy.label}', + color=strategy.color, + marker=strategy.marker, + linewidth=2, + **kwargs, ) ax.legend(frameon=False) @@ -856,6 +873,7 @@ def plot_things_per_things( store=True, ax=None, fig=None, + plotting_args=None, ): # pragma: no cover ''' Plot thingA vs thingB for multiple strategies @@ -887,7 +905,16 @@ def plot_things_per_things( # execute the plots for all strategies for s in strategies: - self.plot_thingA_per_thingB(s, thingA=thingA, thingB=thingB, recovered=recovered, ax=ax, mask=mask, op=op) + self.plot_thingA_per_thingB( + s, + thingA=thingA, + thingB=thingB, + recovered=recovered, + ax=ax, + mask=mask, + op=op, + **(plotting_args if plotting_args else {}), + ) # set the parameters [plt.setp(ax, k, v) for k, v in args.items()] @@ -1627,22 +1654,32 @@ def main(): 'prob': run_quench, 'num_procs': 1, 'mode': 'default', - 'runs': 5000, - 'reload': True, + 'runs': 2000, + 'reload': False, **parse_args(), } - from pySDC.projects.Resilience.strategies import DIRKStrategy, ERKStrategy + + strategy_args = { + 'stop_at_nan': True, + } + + from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError, kAdaptivityStrategy stats_analyser = FaultStats( - # strategies=[BaseStrategy(), AdaptivityStrategy(), IterateStrategy(), HotRodStrategy()], - strategies=[DIRKStrategy(), ERKStrategy()], + strategies=[ + BaseStrategy(**strategy_args), + AdaptivityStrategy(**strategy_args), + kAdaptivityStrategy(**strategy_args), + HotRodStrategy(**strategy_args), + AdaptivityPolynomialError(**strategy_args), + ], faults=[False, True], - recovery_thresh=1.5, + recovery_thresh=1.15, recovery_thresh_abs=RECOVERY_THRESH_ABS.get(kwargs.get('prob', None), 0), stats_path='data/stats-jusuf', **kwargs, ) - stats_analyser.run_stats_generation(runs=kwargs['runs']) + stats_analyser.run_stats_generation(runs=kwargs['runs'], step=12) if MPI.COMM_WORLD.rank > 0: # make sure only one rank accesses the data return None @@ -1682,7 +1719,7 @@ def main(): fig, ax = plt.subplots(1, 1, figsize=(13, 4)) stats_analyser.plot_recovery_thresholds(ax=ax, thresh_range=np.logspace(-1, 1, 1000)) - ax.axvline(stats_analyser.get_thresh(BaseStrategy()), color='grey', ls=':', label='recovery threshold') + # ax.axvline(stats_analyser.get_thresh(BaseStrategy()), color='grey', ls=':', label='recovery threshold') ax.set_xscale('log') ax.legend(frameon=False) fig.tight_layout() @@ -1709,7 +1746,8 @@ def main(): 'error', 'bit', True, op=stats_analyser.mean, mask=mask, args={'yscale': 'log'} ) - stats_analyser.plot_recovery_thresholds() + stats_analyser.plot_recovery_thresholds(thresh_range=np.linspace(0.9, 1.3, 100)) + plt.show() if __name__ == "__main__": diff --git a/pySDC/projects/Resilience/paper_plots.py b/pySDC/projects/Resilience/paper_plots.py index 29426c2d45..90e73020f1 100644 --- a/pySDC/projects/Resilience/paper_plots.py +++ b/pySDC/projects/Resilience/paper_plots.py @@ -17,6 +17,7 @@ HotRodStrategy, DIRKStrategy, ERKStrategy, + AdaptivityPolynomialError, ) from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal from pySDC.helpers.stats_helper import get_sorted @@ -43,7 +44,7 @@ def get_stats(problem, path='data/stats-jusuf', num_procs=1, strategy_type='SDC' if strategy_type == 'SDC': strategies = [BaseStrategy(), AdaptivityStrategy(), IterateStrategy()] if JOURNAL not in ['JSC_beamer']: - strategies += [HotRodStrategy()] + strategies += [HotRodStrategy(), AdaptivityPolynomialError()] elif strategy_type == 'RK': strategies = [DIRKStrategy()] if problem.__name__ in ['run_Lorenz', 'run_vdp']: @@ -54,7 +55,7 @@ def get_stats(problem, path='data/stats-jusuf', num_procs=1, strategy_type='SDC' strategies=strategies, faults=[False, True], reload=True, - recovery_thresh=1.5, + recovery_thresh=1.1, recovery_thresh_abs=RECOVERY_THRESH_ABS.get(problem, 0), mode='default', stats_path=path, @@ -107,12 +108,6 @@ def analyse_resilience(problem, path='data/stats', **kwargs): # pragma: no cove not_overflow = stats_analyser.get_mask(strategy=strategy, key='bit', val=1, op='uneq', old_mask=not_fixed) stats_analyser.print_faults(not_overflow) - # special = stats_analyser.get_mask(strategy=strategy, key='bit', val=10, op='eq') - # stats_analyser.print_faults(special) - - # Adaptivity: 19, ... - # stats_analyser.scrutinize(strategy, run=19, faults=True) - compare_strategies(stats_analyser, **kwargs) plot_recovery_rate(stats_analyser, **kwargs) @@ -146,12 +141,19 @@ def plot_recovery_rate(stats_analyser, **kwargs): # pragma: no cover my_setup_mpl() fig, axs = plt.subplots(1, 2, figsize=(TEXTWIDTH, 5 * cm), sharex=True, sharey=True) stats_analyser.plot_things_per_things( - 'recovered', 'bit', False, op=stats_analyser.rec_rate, args={'ylabel': 'recovery rate'}, ax=axs[0] + 'recovered', + 'bit', + False, + op=stats_analyser.rec_rate, + args={'ylabel': 'recovery rate'}, + plotting_args={'markevery': 5}, + ax=axs[0], ) plot_recovery_rate_recoverable_only(stats_analyser, fig, axs[1], ylabel='') - axs[1].get_legend().remove() + axs[0].get_legend().remove() axs[0].set_title('All faults') axs[1].set_title('Only recoverable faults') + axs[0].set_ylim((-0.05, 1.05)) savefig(fig, 'recovery_rate_compared', **kwargs) @@ -180,6 +182,7 @@ def plot_recovery_rate_recoverable_only(stats_analyser, fig, ax, **kwargs): # p ax=ax, fig=fig, strategies=[stats_analyser.strategies[i]], + plotting_args={'markevery': 5}, ) @@ -209,12 +212,11 @@ def compare_recovery_rate_problems(**kwargs): # pragma: no cover for ax in axs.flatten(): ax.get_legend().remove() - # axs[0, 0].set_ylim((0, 1)) - if kwargs.get('strategy_type', 'SDC') == 'SDC': axs[1, 1].legend(frameon=False) else: axs[0, 1].legend(frameon=False) + axs[0, 0].set_ylim((-0.05, 1.05)) axs[1, 0].set_ylabel('recovery rate') axs[0, 0].set_ylabel('recovery rate') @@ -225,27 +227,6 @@ def compare_recovery_rate_problems(**kwargs): # pragma: no cover savefig(fig, f'compare_equations{name}.pdf') -def plot_recovery_thresholds(problem, num_procs_list=None): - num_procs_list = [1, 4] if num_procs_list is None else num_procs_list - stats_dict = {n: get_stats(problem, num_procs=n) for n in num_procs_list} - strategies = [AdaptivityStrategy(useMPI=True)] - thresh_range = np.linspace(0.9, 1.4, 100) - - fig, ax = plt.subplots() - [ - stats.plot_recovery_thresholds( - strategies=strategies, - thresh_range=thresh_range, - ax=ax, - mask=stats.get_fixable_faults_only(strategies[0]), - label=f'{n} procs', - ls='-' if n == 1 else '--', - ) - for n, stats in stats_dict.items() - ] - savefig(fig, 'threshold', tight_layout=False) - - def plot_adaptivity_stuff(): # pragma: no cover """ Plot the solution for a van der Pol problem as well as the local error and cost associated with the base scheme and @@ -289,23 +270,40 @@ def plot_error(stats, ax, iter_ax, strategy, **kwargs): ax.set_ylabel('local error') iter_ax.set_ylabel(r'Newton iterations') - force_params = {'convergence_controllers': {EstimateEmbeddedError: {}}} - # force_params = {'convergence_controllers': {EstimateEmbeddedError: {}}, 'step_params': {'maxiter': 5}, 'level_params': {'dt': 4e-2}} - for strategy in [BaseStrategy, AdaptivityStrategy, IterateStrategy]: + force_params = {} + for strategy in [BaseStrategy, AdaptivityStrategy, IterateStrategy, AdaptivityPolynomialError]: + if strategy == AdaptivityPolynomialError: + from pySDC.implementations.convergence_controller_classes.adaptivity import ( + AdaptivityPolynomialError as adaptivity, + ) + + force_params = {'sweeper_params': {'num_nodes': 2}} + force_params['convergence_controllers'] = { + adaptivity: { + 'e_tol': 7e-5, + 'restol_rel': 1e-4, + 'restol_min': 1e-10, + 'restart_at_maxiter': True, + 'factor_if_not_converged': 4.0, + }, + } + else: + force_params = {} stats, _, _ = stats_analyser.single_run( - strategy=strategy(), force_params=force_params, hook_class=[LogLocalErrorPostStep, LogData, LogWork] + strategy=strategy(useMPI=False), + force_params=force_params, + hook_class=[LogLocalErrorPostStep, LogData, LogWork], ) plot_error(stats, axs[1], axs[2], strategy()) if strategy == BaseStrategy: u = get_sorted(stats, type='u', recomputed=False) axs[0].plot([me[0] for me in u], [me[1][0] for me in u], color='black', label=r'$u$') - # axs[0].plot([me[0] for me in u], [me[1][1] for me in u], color='black', ls='--', label=r'$u_t$') - # axs[0].legend(frameon=False) axs[2].set_xlabel(r'$t$') axs[0].set_ylabel('solution') axs[2].legend(frameon=JOURNAL == 'JSC_beamer') + axs[1].legend(frameon=True) savefig(fig, 'adaptivity') @@ -394,13 +392,13 @@ def plot_quench_solution(): # pragma: no cover strategy = BaseStrategy() - custom_description = strategy.get_custom_description(run_quench) + custom_description = strategy.get_custom_description(run_quench, num_procs=1) stats, controller, _ = run_quench(custom_description=custom_description, Tend=strategy.get_Tend(run_quench)) prob = controller.MS[0].levels[0].prob - u = get_sorted(stats, type='u') + u = get_sorted(stats, type='u', recomputed=False) ax.plot([me[0] for me in u], [max(me[1]) for me in u], color='black', label='$T$') ax.axhline(prob.u_thresh, label='$T_\mathrm{thresh}$', ls='--', color='grey', zorder=-1) @@ -461,40 +459,43 @@ def work_precision(): # pragma: no cover for mode in ['compare_strategies', 'parallel_efficiency', 'RK_comp']: all_problems(**all_params, mode=mode) - # Quench stuff - fig, axs = get_fig(x=3, y=1, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1, 0.47)) - quench_params = { - **all_params, - 'problem': run_quench, - 'decorate': True, - 'configurations': get_configs('step_size_limiting', run_quench), - 'num_procs': 1, - 'runs': 1, - 'comm_world': MPI.COMM_WORLD, - } - quench_params.pop('base_path', None) - execute_configurations(**{**quench_params, 'work_key': 'k_SDC', 'precision_key': 'k_Newton'}, ax=axs[2]) - execute_configurations(**{**quench_params, 'work_key': 'param', 'precision_key': 'restart'}, ax=axs[1]) - execute_configurations(**{**quench_params, 'work_key': 't', 'precision_key': 'e_global_rel'}, ax=axs[0]) - axs[1].set_yscale('linear') - axs[2].set_yscale('linear') - axs[2].set_xscale('linear') - axs[1].set_xlabel(r'$e_\mathrm{tol}$') - axs[0].set_xticks([1e0, 3e0], minor=True) - - for ax in axs: - ax.set_title(ax.get_ylabel()) - ax.set_ylabel('') - fig.suptitle('Quench') - - save_fig( - fig=fig, - name=f'{run_quench.__name__}', - work_key='step-size', - precision_key='limiting', - legend=True, - base_path=all_params["base_path"], - ) + # # Quench stuff + # fig, axs = get_fig(x=3, y=1, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1, 0.47)) + # quench_params = { + # **all_params, + # 'problem': run_quench, + # 'decorate': True, + # 'configurations': get_configs('step_size_limiting', run_quench), + # 'num_procs': 1, + # 'runs': 1, + # 'comm_world': MPI.COMM_WORLD, + # 'mode': 'step_size_limiting', + # } + # quench_params.pop('base_path', None) + # execute_configurations(**{**quench_params, 'work_key': 'k_SDC', 'precision_key': 'k_Newton'}, ax=axs[2]) + # execute_configurations(**{**quench_params, 'work_key': 'param', 'precision_key': 'restart'}, ax=axs[1]) + # execute_configurations(**{**quench_params, 'work_key': 't', 'precision_key': 'e_global_rel'}, ax=axs[0]) + # axs[1].set_yscale('linear') + # # axs[2].set_yscale('linear') + # axs[2].set_xscale('linear') + # axs[1].set_xlabel(r'$e_\mathrm{tol}$') + # # axs[0].set_xticks([1e0, 3e0], [r'$10^{0}$', r'$3\times 10^{0}$'], minor=False) + + # for ax in axs: + # ax.set_title(ax.get_ylabel()) + # ax.set_ylabel('') + # fig.suptitle('Quench') + + # axs[1].set_yticks([4.0, 6.0, 8.0, 10.0, 12.0], minor=False) + + # save_fig( + # fig=fig, + # name=f'{run_quench.__name__}', + # work_key='step-size', + # precision_key='limiting', + # legend=True, + # base_path=all_params["base_path"], + # ) # End Quench stuff # vdp_stiffness_plot(base_path='data/paper') @@ -539,17 +540,17 @@ def make_plots_for_paper(): # pragma: no cover JOURNAL = 'Springer_Numerical_Algorithms' BASE_PATH = 'data/paper' - # work_precision() - # plot_vdp_solution() - # plot_quench_solution() - # plot_recovery_rate(get_stats(run_vdp)) - # plot_fault_vdp(0) - # plot_fault_vdp(13) plot_adaptivity_stuff() - compare_recovery_rate_problems(num_procs=1, strategy_type='RK') - for i in [1, 4]: - compare_recovery_rate_problems(num_procs=i, strategy_type='SDC') + work_precision() + + plot_vdp_solution() + plot_quench_solution() + + plot_recovery_rate(get_stats(run_vdp)) + plot_fault_vdp(0) + plot_fault_vdp(13) + compare_recovery_rate_problems(num_procs=1, strategy_type='SDC') def make_plots_for_notes(): # pragma: no cover @@ -568,5 +569,4 @@ def make_plots_for_notes(): # pragma: no cover # make_plots_for_notes() # make_plots_for_SIAM_CSE23() # make_plots_for_TIME_X_website() - # plot_recovery_thresholds(run_Schroedinger) make_plots_for_paper() diff --git a/pySDC/projects/Resilience/quench.py b/pySDC/projects/Resilience/quench.py index 7efa261005..9cc1767b5b 100644 --- a/pySDC/projects/Resilience/quench.py +++ b/pySDC/projects/Resilience/quench.py @@ -1,12 +1,11 @@ # script to run a quench problem from pySDC.implementations.problem_classes.Quench import Quench, QuenchIMEX -from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit -from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI from pySDC.core.Hooks import hooks from pySDC.helpers.stats_helper import get_sorted from pySDC.projects.Resilience.hook import hook_collection, LogData from pySDC.projects.Resilience.strategies import merge_descriptions +from pySDC.projects.Resilience.sweepers import imex_1st_order_efficient, generic_implicit_efficient import numpy as np import matplotlib.pyplot as plt @@ -98,7 +97,7 @@ def run_quench( Returns: dict: The stats object controller: The controller - Tend: The time that was supposed to be integrated to + bool: If the code crashed """ if custom_description is not None: problem_params = custom_description.get('problem_params', {}) @@ -120,6 +119,8 @@ def run_quench( problem_params = { 'newton_tol': 1e-9, 'direct_solver': False, + 'order': 6, + 'nvars': 2**7, } # initialize step parameters @@ -139,7 +140,7 @@ def run_quench( description = {} description['problem_class'] = QuenchIMEX if imex else Quench description['problem_params'] = problem_params - description['sweeper_class'] = imex_1st_order if imex else generic_implicit + description['sweeper_class'] = imex_1st_order_efficient if imex else generic_implicit_efficient description['sweeper_params'] = sweeper_params description['level_params'] = level_params description['step_params'] = step_params @@ -155,6 +156,7 @@ def run_quench( 'controller_params': controller_params, 'description': description, } + if use_MPI: from mpi4py import MPI from pySDC.implementations.controller_classes.controller_MPI import controller_MPI @@ -175,12 +177,14 @@ def run_quench( prepare_controller_for_faults(controller, fault_stuff) # call main function to get things done... + crash = False try: uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - except ConvergenceError: - print('Warning: Premature termination!') + except ConvergenceError as e: + print(f'Warning: Premature termination! {e}') stats = controller.return_stats() - return stats, controller, Tend + crash = True + return stats, controller, crash def faults(seed=0): # pragma: no cover @@ -377,8 +381,8 @@ def compare_imex_full(plotting=False, leak_type='linear'): assert error[True] < 1.2e-4, f'Expected error of IMEX version to be less than 1.2e-4, but got e={error[True]:.2e}!' assert ( - error[False] < 7.7e-5 - ), f'Expected error of fully implicit version to be less than 7.7e-5, but got e={error[False]:.2e}!' + error[False] < 8e-5 + ), f'Expected error of fully implicit version to be less than 8e-5, but got e={error[False]:.2e}!' def compare_reference_solutions_single(): @@ -455,7 +459,6 @@ def compare_reference_solutions(): fig, ax = plt.subplots() Tend = 500 dt_list = [Tend / 2.0**me for me in [2, 3, 4, 5, 6, 7, 8, 9, 10]] - # dt_list = [Tend / 2.**me for me in [2, 3, 4, 5, 6, 7]] for j in range(len(types)): errors = [None] * len(dt_list) @@ -495,13 +498,14 @@ def check_order(reference_sol_type='scipy'): Tend = 500 maxiter_list = [1, 2, 3, 4, 5] - dt_list = [Tend / 2.0**me for me in [4, 5, 6, 7, 8, 9]] - # dt_list = [Tend / 2.**me for me in [6, 7, 8]] + dt_list = [Tend / 2.0**me for me in [4, 5, 6, 7, 8]] fig, ax = plt.subplots() from pySDC.implementations.sweeper_classes.Runge_Kutta import DIRK43 + controller_params = {'logger_level': 30} + colors = ['black', 'teal', 'magenta', 'orange', 'red'] for j in range(len(maxiter_list)): errors = [None] * len(dt_list) @@ -514,8 +518,11 @@ def check_order(reference_sol_type='scipy'): description['problem_params'] = { 'leak_type': 'linear', 'leak_transition': 'step', - 'nvars': 2**10, + 'nvars': 2**7, 'reference_sol_type': reference_sol_type, + 'newton_tol': 1e-10, + 'lintol': 1e-10, + 'direct_solver': True, } description['convergence_controllers'] = {EstimateEmbeddedError: {}} @@ -523,11 +530,18 @@ def check_order(reference_sol_type='scipy'): # description['sweeper_class'] = DIRK43 # description['sweeper_params'] = {'maxiter': 1} + hook_class = [ + # LogGlobalErrorPostRun, + ] stats, controller, _ = run_quench( - custom_description=description, hook_class=[LogGlobalErrorPostRun], Tend=Tend, imex=False + custom_description=description, + hook_class=hook_class, + Tend=Tend, + imex=False, + custom_controller_params=controller_params, ) - # errors[i] = max([me[1] for me in get_sorted(stats, type='error_embedded_estimate')]) - errors[i] = get_sorted(stats, type='e_global_post_run')[-1][1] + errors[i] = max([me[1] for me in get_sorted(stats, type='error_embedded_estimate')]) + # errors[i] = get_sorted(stats, type='e_global_post_run')[-1][1] print(errors) ax.loglog(dt_list, errors, color=colors[j], label=f'{maxiter_list[j]} iterations') ax.loglog( @@ -550,10 +564,11 @@ def check_order(reference_sol_type='scipy'): if __name__ == '__main__': - compare_reference_solutions_single() - # for reference_sol_type in ['DIRK', 'SDC', 'scipy']: - # check_order(reference_sol_type=reference_sol_type) + # compare_reference_solutions_single() + for reference_sol_type in ['scipy']: + check_order(reference_sol_type=reference_sol_type) # faults(19) # get_crossing_time() # compare_imex_full(plotting=True) + # iteration_counts() plt.show() diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py index 654d1492d5..8b235251a8 100644 --- a/pySDC/projects/Resilience/strategies.py +++ b/pySDC/projects/Resilience/strategies.py @@ -31,12 +31,12 @@ class Strategy: Abstract class for resilience strategies ''' - def __init__(self, useMPI=False, skip_residual_computation='none'): + def __init__(self, useMPI=False, skip_residual_computation='none', stop_at_nan=True, **kwargs): ''' Initialization routine ''' self.useMPI = useMPI - self.max_steps = 1e4 + self.max_steps = 1e5 # set default values for plotting self.linestyle = '-' @@ -50,12 +50,22 @@ def __init__(self, useMPI=False, skip_residual_computation='none'): self.skip_residual_computation = ('IT_CHECK', 'IT_DOWN', 'IT_UP', 'IT_FINE', 'IT_COARSE') elif skip_residual_computation == 'most': self.skip_residual_computation = ('IT_DOWN', 'IT_UP', 'IT_FINE', 'IT_COARSE') - else: + elif skip_residual_computation == 'none': self.skip_residual_computation = () + else: + raise NotImplementedError( + f'Don\'t know when to skip residual computation with rule \"{skip_residual_computation}\"' + ) + + self.stop_at_nan = stop_at_nan # setup custom descriptions self.custom_description = {} self.custom_description['sweeper_params'] = {'skip_residual_computation': self.skip_residual_computation} + self.custom_description['level_params'] = {} + self.custom_description['problem_params'] = {} + self.custom_description['step_params'] = {} + self.custom_description['convergence_controllers'] = {} # prepare parameters for masks to identify faults that cannot be fixed by this strategy self.fixable = [] @@ -80,6 +90,12 @@ def __init__(self, useMPI=False, skip_residual_computation='none'): def __str__(self): return self.name + def get_controller_params(self, **kwargs): + return {'all_to_done': False} + + def get_description_for_tolerance(self, problem, param, **kwargs): + return {} + def get_fixable_params(self, **kwargs): """ Return a list containing dictionaries which can be passed to `FaultStats.get_mask` as keyword arguments to @@ -109,7 +125,7 @@ def get_fault_args(self, problem, num_procs): elif problem.__name__ == "run_Schroedinger": args['time'] = 0.3 elif problem.__name__ == "run_quench": - args['time'] = 31.0 + args['time'] = 41.0 elif problem.__name__ == "run_Lorenz": args['time'] = 0.3 @@ -136,7 +152,9 @@ def get_random_params(self, problem, num_procs): rnd_params['min_node'] = 1 if problem.__name__ == "run_quench": - rnd_params['iteration'] = 1 + rnd_params['iteration'] = 5 + elif problem.__name__ == 'run_Lorenz': + rnd_params['iteration'] = 5 return rnd_params @property @@ -176,7 +194,6 @@ def get_Tend(cls, problem, num_procs=1): ''' if problem.__name__ == "run_vdp": return 11.5 - # return 2.3752559741400825 # old stuff elif problem.__name__ == "run_piline": return 20.0 elif problem.__name__ == "run_Lorenz": @@ -185,6 +202,8 @@ def get_Tend(cls, problem, num_procs=1): return 1.0 elif problem.__name__ == "run_quench": return 500.0 + elif problem.__name__ == "run_AC": + return 1e-2 else: raise NotImplementedError('I don\'t have a final time for your problem!') @@ -202,29 +221,65 @@ def get_base_parameters(self, problem, num_procs=1): from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter custom_description = {} + custom_description['step_params'] = {} + custom_description['level_params'] = {} + custom_description['problem_params'] = {} + if problem.__name__ == "run_vdp": custom_description['step_params'] = {'maxiter': 3} custom_description['problem_params'] = { 'u0': np.array([2, 0], dtype=np.float64), - # 'u0': np.array([0.99995, -0.00999985], dtype=np.float64), # old stuff 'crash_at_maxiter': False, 'newton_tol': 1e-11, + 'stop_at_nan': False, } custom_description['level_params'] = {'dt': 1e-2} elif problem.__name__ == "run_Lorenz": custom_description['step_params'] = {'maxiter': 5} custom_description['level_params'] = {'dt': 1e-2} + custom_description['problem_params'] = {'stop_at_nan': False} elif problem.__name__ == "run_Schroedinger": custom_description['step_params'] = {'maxiter': 5} custom_description['level_params'] = {'dt': 1e-2, 'restol': -1} + custom_description['problem_params']['nvars'] = (256, 256) elif problem.__name__ == "run_quench": custom_description['level_params'] = {'restol': -1, 'dt': 8.0} custom_description['step_params'] = {'maxiter': 5} - custom_description['problem_params'] = {'newton_maxiter': 99, 'newton_tol': 1e-11, 'nvars': 2**7} + custom_description['problem_params'] = { + 'newton_maxiter': 29, + 'newton_tol': 1e-7, + 'nvars': 2**6, + 'direct_solver': False, + 'lintol': 1e-8, + 'liniter': 29, + 'order': 6, + } + elif problem.__name__ == "run_AC": + custom_description['level_params'] = {'restol': -1, 'dt': 1e-4} + custom_description['step_params'] = {'maxiter': 5} + custom_description['problem_params'] = {'newton_maxiter': 99, 'newton_tol': 1e-11} custom_description['convergence_controllers'] = { - StepSizeLimiter: {'dt_min': self.get_Tend(problem=problem, num_procs=num_procs) / self.max_steps} + # StepSizeLimiter: {'dt_min': self.get_Tend(problem=problem, num_procs=num_procs) / self.max_steps} + } + + if self.stop_at_nan: + from pySDC.implementations.convergence_controller_classes.crash import StopAtNan + + custom_description['convergence_controllers'][StopAtNan] = {'thresh': 1e20} + + from pySDC.implementations.convergence_controller_classes.crash import StopAtMaxRuntime + + max_runtime = { + 'run_vdp': 60, + 'run_Lorenz': 60, + 'run_Schroedinger': 150, + 'run_quench': 150, + } + + custom_description['convergence_controllers'][StopAtMaxRuntime] = { + 'max_runtime': max_runtime.get(problem.__name__, 100) } return custom_description @@ -242,6 +297,15 @@ def get_custom_description(self, problem, num_procs=1): custom_description = self.get_base_parameters(problem, num_procs) return merge_descriptions(custom_description, self.custom_description) + def get_custom_description_for_faults(self, *args, **kwargs): + ''' + Get a custom description based on the problem to run the fault stuff + + Returns: + dict: Custom description + ''' + return self.get_custom_description(*args, **kwargs) + def get_reference_value(self, problem, key, op, num_procs=1): """ Get a reference value for a given problem for testing in CI. @@ -258,16 +322,80 @@ def get_reference_value(self, problem, key, op, num_procs=1): raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') +class InexactBaseStrategy(Strategy): + """ + Base class for inexact strategies. + """ + + def __init__( + self, double_adaptivity=False, newton_inexactness=True, linear_inexactness=True, SDC_maxiter=16, **kwargs + ): + kwargs = {**kwargs, 'skip_residual_computation': 'most'} + super().__init__(**kwargs) + self.double_adaptivity = double_adaptivity + self.newton_inexactness = newton_inexactness + self.linear_inexactness = linear_inexactness + self.SDC_maxiter = SDC_maxiter + + def get_controller_params(self, **kwargs): + return {'all_to_done': True} + + def get_custom_description(self, problem, num_procs=1): + from pySDC.implementations.convergence_controller_classes.inexactness import NewtonInexactness + + preconditioner = 'MIN-SR-NS' if problem.__name__ in ['run_vdp', 'run_Lorenz'] else 'MIN-SR-S' + + desc = {} + desc['sweeper_params'] = {'QI': preconditioner} + desc['step_params'] = {'maxiter': self.SDC_maxiter} + desc['problem_params'] = {} + desc['level_params'] = {'restol': 1e-8, 'residual_type': 'last_abs'} + desc['convergence_controllers'] = {} + + inexactness_params = { + 'min_tol': 1e-12, + 'ratio': 1e-2, + 'max_tol': 1e-4, + 'use_e_tol': False, + 'maxiter': 15, + } + + if self.newton_inexactness and problem.__name__ not in ['run_Schroedinger']: + if problem.__name__ == 'run_quench': + inexactness_params['ratio'] = 1e-1 + inexactness_params['min_tol'] = 1e-11 + inexactness_params['maxiter'] = 5 + desc['convergence_controllers'][NewtonInexactness] = inexactness_params + + if problem.__name__ in ['run_vdp']: + desc['problem_params']['stop_at_nan'] = False + + if self.linear_inexactness and problem.__name__ in ['run_AC', 'run_quench']: + desc['problem_params']['inexact_linear_ratio'] = 1e-1 + if problem.__name__ in ['run_quench']: + desc['problem_params']['direct_solver'] = False + desc['problem_params']['liniter'] = 9 + desc['problem_params']['min_lintol'] = 1e-11 + + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestarting + + desc['convergence_controllers'][BasicRestarting.get_implementation(useMPI=self.useMPI)] = { + 'max_restarts': 29, + 'crash_after_max_restarts': True, + } + return merge_descriptions(super().get_custom_description(problem, num_procs), desc) + + class BaseStrategy(Strategy): ''' Do a fixed iteration count ''' - def __init__(self, useMPI=False, skip_residual_computation='all'): + def __init__(self, skip_residual_computation='all', **kwargs): ''' Initialization routine ''' - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + super().__init__(skip_residual_computation=skip_residual_computation, **kwargs) self.color = list(cmap.values())[0] self.marker = 'o' self.name = 'base' @@ -279,6 +407,12 @@ def __init__(self, useMPI=False, skip_residual_computation='all'): def label(self): return r'fixed' + def get_custom_description_for_faults(self, problem, *args, **kwargs): + desc = self.get_custom_description(problem, *args, **kwargs) + if problem.__name__ == "run_quench": + desc['level_params']['dt'] = 5.0 + return desc + def get_reference_value(self, problem, key, op, num_procs=1): """ Get a reference value for a given problem for testing in CI. @@ -306,13 +440,14 @@ class AdaptivityStrategy(Strategy): Adaptivity as a resilience strategy ''' - def __init__(self, useMPI=False, skip_residual_computation='all'): + def __init__(self, **kwargs): ''' Initialization routine ''' from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + kwargs['skip_residual_computation'] = 'all' + super().__init__(**kwargs) self.color = list(cmap.values())[1] self.marker = '*' self.name = 'adaptivity' @@ -355,37 +490,36 @@ def get_custom_description(self, problem, num_procs): The custom descriptions you can supply to the problem when running it ''' from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter custom_description = {} custom_description['convergence_controllers'] = {} dt_max = np.inf - dt_min = 1e-5 dt_slope_max = np.inf if problem.__name__ == "run_piline": e_tol = 1e-7 - dt_min = 1e-2 elif problem.__name__ == "run_vdp": e_tol = 2e-5 - dt_min = 1e-3 elif problem.__name__ == "run_Lorenz": e_tol = 2e-5 - dt_min = 1e-3 elif problem.__name__ == "run_Schroedinger": e_tol = 4e-6 - dt_min = 1e-3 elif problem.__name__ == "run_quench": - e_tol = 1e-5 - dt_min = 1e-3 - # dt_max = 25 - # dt_slope_max = 4. + e_tol = 1e-8 + custom_description['problem_params'] = { + 'newton_tol': 1e-9, + 'lintol': 1e-10, + } from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestarting custom_description['convergence_controllers'][BasicRestarting.get_implementation(useMPI=self.useMPI)] = { - 'max_restarts': 15, + 'max_restarts': 99, } + elif problem.__name__ == "run_AC": + e_tol = 1e-4 else: raise NotImplementedError( 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\ @@ -394,10 +528,11 @@ def get_custom_description(self, problem, num_procs): custom_description['convergence_controllers'][Adaptivity] = { 'e_tol': e_tol, - 'dt_min': dt_min, - 'dt_max': dt_max, 'dt_slope_max': dt_slope_max, } + custom_description['convergence_controllers'][StepSizeLimiter] = { + 'dt_max': dt_max, + } return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) def get_reference_value(self, problem, key, op, num_procs=1): @@ -421,10 +556,20 @@ def get_reference_value(self, problem, key, op, num_procs=1): raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') + def get_custom_description_for_faults(self, problem, num_procs, *args, **kwargs): + from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + + desc = self.get_custom_description(problem, num_procs, *args, **kwargs) + if problem.__name__ == "run_quench": + desc['level_params']['dt'] = 1.1e1 + desc['convergence_controllers'][Adaptivity]['e_tol'] = 1e-6 + return desc + class AdaptivityRestartFirstStep(AdaptivityStrategy): - def __init__(self, useMPI=False): - super().__init__(useMPI=useMPI) + def __init__(self, **kwargs): + super().__init__(**kwargs) self.color = 'teal' self.name = 'adaptivityRestartFirstStep' @@ -458,13 +603,14 @@ class AdaptiveHotRodStrategy(Strategy): Adaptivity + Hot Rod as a resilience strategy ''' - def __init__(self, useMPI=False, skip_residual_computation='all'): + def __init__(self, **kwargs): ''' Initialization routine ''' from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + kwargs['skip_residual_computation'] = 'all' + super().__init__(**kwargs) self.color = list(cmap.values())[4] self.marker = '.' self.name = 'adaptive Hot Rod' @@ -536,11 +682,12 @@ class IterateStrategy(Strategy): Iterate for as much as you want ''' - def __init__(self, useMPI=False, skip_residual_computation='most'): + def __init__(self, **kwargs): ''' Initialization routine ''' - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + kwargs['skip_residual_computation'] = 'most' + super().__init__(**kwargs) self.color = list(cmap.values())[2] self.marker = 'v' self.name = 'iterate' @@ -588,10 +735,27 @@ def get_custom_description(self, problem, num_procs): } if problem.__name__ == "run_quench": - custom_description['level_params']['dt'] = 1 + custom_description['level_params']['dt'] = 1.0 return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + def get_random_params(self, problem, num_procs): + ''' + Routine to get parameters for the randomization of faults + + Args: + problem: A function that runs a pySDC problem, see imports for available problems + num_procs (int): Number of processes you intend to run with + + Returns: + dict: Randomization parameters + ''' + + rnd_params = super().get_random_params(problem, num_procs) + if problem.__name__ == "run_quench": + rnd_params['iteration'] = 1 + return rnd_params + def get_reference_value(self, problem, key, op, num_procs=1): """ Get a reference value for a given problem for testing in CI. @@ -614,16 +778,60 @@ def get_reference_value(self, problem, key, op, num_procs=1): raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') +class kAdaptivityStrategy(IterateStrategy): + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.precision_parameter = 'dt' + self.precision_parameter_loc = ['level_params', 'dt'] + + def get_custom_description(self, problem, num_procs, *args, **kwargs): + desc = super().get_custom_description(problem, num_procs, *args, **kwargs) + desc['level_params']['restol'] = 1e-9 + if problem.__name__ == "run_quench": + desc['problem_params']['newton_tol'] = 1e-9 + desc['problem_params']['lintol'] = 1e-9 + desc['level_params']['dt'] = 2.5 + return desc + + def get_custom_description_for_faults(self, problem, *args, **kwargs): + desc = self.get_custom_description(problem, *args, **kwargs) + if problem.__name__ == 'run_quench': + desc['level_params']['dt'] = 5.0 + return desc + + def get_reference_value(self, problem, key, op, num_procs=1): + """ + Get a reference value for a given problem for testing in CI. + + Args: + problem: A function that runs a pySDC problem, see imports for available problems + key (str): The name of the variable you want to compare + op (function): The operation you want to apply to the data + num_procs (int): Number of processes + + Returns: + The reference value + """ + if problem.__name__ == "run_vdp": + if key == 'work_newton' and op == sum: + return 13241 + elif key == 'e_global_post_run' and op == max: + return 1.0505165626284452e-08 + + raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') + + class HotRodStrategy(Strategy): ''' Hot Rod as a resilience strategy ''' - def __init__(self, useMPI=False, skip_residual_computation='all'): + def __init__(self, **kwargs): ''' Initialization routine ''' - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + kwargs['skip_residual_computation'] = 'all' + super().__init__(**kwargs) self.color = list(cmap.values())[3] self.marker = '^' self.name = 'Hot Rod' @@ -651,7 +859,7 @@ def get_custom_description(self, problem, num_procs): elif num_procs == 5: HotRod_tol = 9.329361e-05 else: # 1 process - HotRod_tol = 1.347949e-06 # 5e-7 + HotRod_tol = 1.347949e-06 HotRod_tol = 7e-6 if num_procs > 1 else 5e-7 maxiter = 4 elif problem.__name__ == "run_Lorenz": @@ -660,7 +868,7 @@ def get_custom_description(self, problem, num_procs): elif num_procs == 4: HotRod_tol = 3.201e-6 else: - HotRod_tol = 7.720589e-07 # 4e-7 + HotRod_tol = 7.720589e-07 maxiter = 6 elif problem.__name__ == "run_Schroedinger": if num_procs == 5: @@ -700,6 +908,12 @@ def get_custom_description(self, problem, num_procs): return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + def get_custom_description_for_faults(self, problem, *args, **kwargs): + desc = self.get_custom_description(problem, *args, **kwargs) + if problem.__name__ == "run_quench": + desc['level_params']['dt'] = 5.0 + return desc + def get_reference_value(self, problem, key, op, num_procs=1): """ Get a reference value for a given problem for testing in CI. @@ -722,18 +936,24 @@ def get_reference_value(self, problem, key, op, num_procs=1): raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') -class AdaptivityCollocationStrategy(Strategy): +class AdaptivityCollocationStrategy(InexactBaseStrategy): ''' Adaptivity based on collocation as a resilience strategy ''' - def __init__(self, useMPI=False, skip_residual_computation='most'): + def __init__(self, **kwargs): ''' Initialization routine ''' + kwargs = { + 'skip_residual_computation': 'most', + **kwargs, + } + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityCollocation - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + self.restol = None + super().__init__(**kwargs) self.color = list(cmap.values())[1] self.marker = '*' self.name = 'adaptivity_coll' @@ -741,8 +961,6 @@ def __init__(self, useMPI=False, skip_residual_computation='most'): self.precision_parameter = 'e_tol' self.adaptive_coll_params = {} self.precision_parameter_loc = ['convergence_controllers', AdaptivityCollocation, 'e_tol'] - self.restol = None - self.maxiter = 99 def get_custom_description(self, problem, num_procs): ''' @@ -758,7 +976,6 @@ def get_custom_description(self, problem, num_procs): from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityCollocation custom_description = {} - custom_description['step_params'] = {'maxiter': self.maxiter} dt_max = np.inf dt_min = 1e-5 @@ -779,27 +996,29 @@ def get_custom_description(self, problem, num_procs): e_tol = 1e-5 dt_min = 1e-3 dt_max = 1e2 + elif problem.__name__ == "run_AC": + e_tol = 1e-4 else: raise NotImplementedError( 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\ strategy' ) - custom_description['level_params'] = {'restol': e_tol / 10 if self.restol is None else self.restol} custom_description['convergence_controllers'] = { AdaptivityCollocation: { 'e_tol': e_tol, 'dt_min': dt_min, 'dt_max': dt_max, 'adaptive_coll_params': self.adaptive_coll_params, + 'restol_rel': 1e-2, } } return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) class AdaptivityCollocationTypeStrategy(AdaptivityCollocationStrategy): - def __init__(self, useMPI=False, skip_residual_computation='most'): - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + def __init__(self, **kwargs): + super().__init__(**kwargs) self.color = list(cmap.values())[4] self.marker = '.' self.adaptive_coll_params = { @@ -826,16 +1045,16 @@ def get_reference_value(self, problem, key, op, num_procs=1): """ if problem.__name__ == "run_vdp": if key == 'work_newton' and op == sum: - return 2694 + return 1954 elif key == 'e_global_post_run' and op == max: - return 2.1707816100224875e-06 + return 1.1341277847964903e-06 raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') class AdaptivityCollocationRefinementStrategy(AdaptivityCollocationStrategy): - def __init__(self, useMPI=False, skip_residual_computation='most'): - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + def __init__(self, **kwargs): + super().__init__(**kwargs) self.color = list(cmap.values())[5] self.marker = '^' self.adaptive_coll_params = { @@ -863,16 +1082,16 @@ def get_reference_value(self, problem, key, op, num_procs=1): """ if problem.__name__ == "run_vdp": if key == 'work_newton' and op == sum: - return 1881 + return 1785 elif key == 'e_global_post_run' and op == max: - return 3.3428689164005654e-06 + return 1.2448609458259874e-06 raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') class AdaptivityCollocationDerefinementStrategy(AdaptivityCollocationStrategy): - def __init__(self, useMPI=False, skip_residual_computation='most'): - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + def __init__(self, **kwargs): + super().__init__(**kwargs) self.color = list(cmap.values())[6] self.marker = '^' self.adaptive_coll_params = {'num_nodes': [4, 3]} @@ -896,9 +1115,9 @@ def get_reference_value(self, problem, key, op, num_procs=1): """ if problem.__name__ == "run_vdp": if key == 'work_newton' and op == sum: - return 3421 + return 2508 elif key == 'e_global_post_run' and op == max: - return 2.1130961994131336e-05 + return 3.452255036812124e-05 raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') @@ -908,13 +1127,13 @@ class DIRKStrategy(AdaptivityStrategy): DIRK4(3) ''' - def __init__(self, useMPI=False, skip_residual_computation='all'): + def __init__(self, **kwargs): ''' Initialization routine ''' from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + super().__init__(**kwargs) self.color = list(cmap.values())[7] self.marker = '^' self.name = 'DIRK' @@ -1003,18 +1222,101 @@ def get_random_params(self, problem, num_procs): return rnd_params +class ARKStrategy(AdaptivityStrategy): + ''' + ARK5(4) + ''' + + def __init__(self, **kwargs): + ''' + Initialization routine + ''' + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK + + super().__init__(**kwargs) + self.color = list(cmap.values())[7] + self.marker = 'P' + self.name = 'ARK' + self.bar_plot_x_label = 'ARK5(4)' + self.precision_parameter = 'e_tol' + self.precision_parameter_loc = ['convergence_controllers', AdaptivityRK, 'e_tol'] + self.max_steps = 1e5 + + @property + def label(self): + return 'ARK5(4)' + + def get_custom_description(self, problem, num_procs): + ''' + Routine to get a custom description that adds adaptivity + + Args: + problem: A function that runs a pySDC problem, see imports for available problems + num_procs (int): Number of processes you intend to run with + + Returns: + The custom descriptions you can supply to the problem when running it + ''' + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK, Adaptivity + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestarting + from pySDC.implementations.sweeper_classes.Runge_Kutta import ARK548L2SA + + adaptivity_description = super().get_custom_description(problem, num_procs) + + e_tol = adaptivity_description['convergence_controllers'][Adaptivity]['e_tol'] / 20.0 + adaptivity_description['convergence_controllers'].pop(Adaptivity, None) + adaptivity_description.pop('sweeper_params', None) + + rk_params = { + 'step_params': {'maxiter': 1}, + 'sweeper_class': ARK548L2SA, + 'convergence_controllers': { + AdaptivityRK: {'e_tol': e_tol}, + BasicRestarting.get_implementation(useMPI=self.useMPI): { + 'max_restarts': 49, + 'crash_after_max_restarts': False, + }, + }, + } + + custom_description = merge_descriptions(adaptivity_description, rk_params) + + return custom_description + + def get_reference_value(self, problem, key, op, num_procs=1): + """ + Get a reference value for a given problem for testing in CI. + + Args: + problem: A function that runs a pySDC problem, see imports for available problems + key (str): The name of the variable you want to compare + op (function): The operation you want to apply to the data + num_procs (int): Number of processes + + Returns: + The reference value + """ + if problem.__name__ == "run_Schroedinger": + if key == 'work_newton' and op == sum: + return 0 + elif key == 'e_global_post_run' and op == max: + return 2.444605657738645e-07 + + raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') + + class ESDIRKStrategy(AdaptivityStrategy): ''' ESDIRK5(3) ''' - def __init__(self, useMPI=False, skip_residual_computation='all'): + def __init__(self, **kwargs): ''' Initialization routine ''' from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + super().__init__(**kwargs) self.color = 'violet' self.marker = '^' self.name = 'ESDIRK' @@ -1027,6 +1329,12 @@ def __init__(self, useMPI=False, skip_residual_computation='all'): def label(self): return 'ESDIRK5(3)' + def get_description_for_tolerance(self, problem, param, **kwargs): + desc = {} + if problem.__name__ == 'run_Schroedinger': + desc['problem_params'] = {'lintol': param} + return desc + def get_custom_description(self, problem, num_procs): ''' Routine to get a custom description that adds adaptivity @@ -1048,11 +1356,13 @@ def get_custom_description(self, problem, num_procs): adaptivity_description['convergence_controllers'].pop(Adaptivity, None) adaptivity_description.pop('sweeper_params', None) + mod = 1e1 if problem.__name__ == 'run_quench' else 1.0 + rk_params = { 'step_params': {'maxiter': 1}, 'sweeper_class': ESDIRK53, 'convergence_controllers': { - AdaptivityRK: {'e_tol': e_tol}, + AdaptivityRK: {'e_tol': e_tol * mod}, BasicRestarting.get_implementation(useMPI=self.useMPI): { 'max_restarts': 49, 'crash_after_max_restarts': False, @@ -1108,16 +1418,22 @@ class ERKStrategy(DIRKStrategy): Explicit embedded RK using Cash-Karp's method """ - def __init__(self, useMPI=False, skip_residual_computation='all'): + def __init__(self, **kwargs): ''' Initialization routine ''' - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) - self.color = list(cmap.values())[9] + super().__init__(**kwargs) + self.color = list(cmap.values())[8] self.marker = 'x' self.name = 'ERK' self.bar_plot_x_label = 'ERK5(4)' + def get_description_for_tolerance(self, problem, param, **kwargs): + desc = {} + if problem.__name__ == 'run_Schroedinger': + desc['problem_params'] = {'lintol': param} + return desc + @property def label(self): return 'CP5(4)' @@ -1172,13 +1488,14 @@ class DoubleAdaptivityStrategy(AdaptivityStrategy): Adaptivity based both on embedded estimate and on residual ''' - def __init__(self, useMPI=False, skip_residual_computation='all'): + def __init__(self, **kwargs): ''' Initialization routine ''' from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity - super().__init__(useMPI=useMPI, skip_residual_computation=skip_residual_computation) + kwargs['skip_residual_computation'] = 'all' + super().__init__(**kwargs) self.color = list(cmap.values())[7] self.marker = '^' self.name = 'double_adaptivity' @@ -1359,18 +1676,19 @@ def get_reference_value(self, problem, key, op, num_procs=1): raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') -class AdaptivityExtrapolationWithinQStrategy(Strategy): +class AdaptivityExtrapolationWithinQStrategy(InexactBaseStrategy): ''' Adaptivity based on extrapolation between collocation nodes as a resilience strategy ''' - def __init__(self, useMPI=False): + def __init__(self, **kwargs): ''' Initialization routine ''' from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityExtrapolationWithinQ - super().__init__(useMPI=useMPI) + self.restol = None + super().__init__(**kwargs) self.color = list(cmap.values())[8] self.marker = '*' self.name = 'adaptivity_extraQ' @@ -1378,8 +1696,6 @@ def __init__(self, useMPI=False): self.precision_parameter = 'e_tol' self.adaptive_coll_params = {} self.precision_parameter_loc = ['convergence_controllers', AdaptivityExtrapolationWithinQ, 'e_tol'] - self.restol = None - self.maxiter = 99 def get_custom_description(self, problem, num_procs): ''' @@ -1395,7 +1711,6 @@ def get_custom_description(self, problem, num_procs): from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityExtrapolationWithinQ custom_description = {} - custom_description['step_params'] = {'maxiter': self.maxiter} dt_max = np.inf dt_min = 1e-5 @@ -1403,43 +1718,36 @@ def get_custom_description(self, problem, num_procs): if problem.__name__ == "run_vdp": e_tol = 2e-5 dt_min = 1e-3 - # elif problem.__name__ == "run_piline": - # e_tol = 1e-7 - # dt_min = 1e-2 - # elif problem.__name__ == "run_Lorenz": - # e_tol = 2e-5 - # dt_min = 1e-3 - # elif problem.__name__ == "run_Schroedinger": - # e_tol = 4e-6 - # dt_min = 1e-3 - # elif problem.__name__ == "run_quench": - # e_tol = 1e-5 - # dt_min = 1e-3 - # dt_max = 1e2 + elif problem.__name__ == "run_piline": + e_tol = 1e-7 + dt_min = 1e-2 + elif problem.__name__ == "run_Lorenz": + e_tol = 2e-5 + dt_min = 1e-3 + elif problem.__name__ == "run_Schroedinger": + e_tol = 4e-6 + dt_min = 1e-3 + elif problem.__name__ == "run_quench": + e_tol = 1e-5 + dt_min = 1e-3 + dt_max = 1e2 + elif problem.__name__ == "run_AC": + e_tol = 1e-4 else: raise NotImplementedError( 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\ strategy' ) - if problem.__name__ in ['run_Schroedinger', 'run_quench']: - from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order - - sweeper_class = imex_1st_order - else: - from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit - - sweeper_class = generic_implicit - - custom_description['level_params'] = {'restol': e_tol / 10 if self.restol is None else self.restol} custom_description['convergence_controllers'] = { AdaptivityExtrapolationWithinQ: { 'e_tol': e_tol, 'dt_min': dt_min, 'dt_max': dt_max, + 'restol_rel': 1e-2, + 'restart_at_maxiter': True, } } - custom_description['sweeper_class'] = sweeper_class return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) def get_reference_value(self, problem, key, op, num_procs=1): @@ -1457,8 +1765,146 @@ def get_reference_value(self, problem, key, op, num_procs=1): """ if problem.__name__ == "run_vdp": if key == 'work_newton' and op == sum: - return 3443 + return 2426 elif key == 'e_global_post_run' and op == max: - return 4.929282266752377e-06 + return 1.7293825160247245e-06 raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') + + +class AdaptivityPolynomialError(InexactBaseStrategy): + ''' + Adaptivity based on extrapolation between collocation nodes as a resilience strategy + ''' + + def __init__(self, interpolate_between_restarts=False, use_restol_rel=True, **kwargs): + ''' + Initialization routine + ''' + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError + + self.restol = None + super().__init__(**kwargs) + self.color = list(cmap.values())[9] + self.marker = '+' + self.name = 'adaptivity-inter' + self.bar_plot_x_label = 'adaptivity Q' + self.precision_parameter = 'e_tol' + self.adaptive_coll_params = {} + self.precision_parameter_loc = ['convergence_controllers', AdaptivityPolynomialError, 'e_tol'] + self.interpolate_between_restarts = interpolate_between_restarts + self.use_restol_rel = use_restol_rel + + def get_custom_description(self, problem, num_procs): + ''' + Routine to get a custom description that adds adaptivity + + Args: + problem: A function that runs a pySDC problem, see imports for available problems + num_procs (int): Number of processes you intend to run with + + Returns: + The custom descriptions you can supply to the problem when running it + ''' + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError + from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter + + custom_description = {} + + dt_max = np.inf + max_slope = 4.0 + restol_rel = 1e-4 + restol_min = 1e-12 + level_params = {} + + if problem.__name__ == "run_vdp": + e_tol = 6e-4 + level_params['dt'] = 0.1 + elif problem.__name__ == "run_piline": + e_tol = 1e-7 + elif problem.__name__ == "run_Lorenz": + e_tol = 7e-4 + elif problem.__name__ == "run_Schroedinger": + e_tol = 3e-4 + elif problem.__name__ == "run_quench": + e_tol = 1e-7 + level_params['dt'] = 50.0 + restol_min = 1e-11 + restol_rel = 1e-1 + elif problem.__name__ == "run_AC": + e_tol = 1e-4 + else: + raise NotImplementedError( + 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\ + strategy' + ) + + custom_description['convergence_controllers'] = { + AdaptivityPolynomialError: { + 'e_tol': e_tol, + 'restol_rel': restol_rel if self.use_restol_rel else 1e-11, + 'restol_min': restol_min if self.use_restol_rel else 1e-12, + 'restart_at_maxiter': True, + 'factor_if_not_converged': max_slope, + 'interpolate_between_restarts': self.interpolate_between_restarts, + }, + StepSizeLimiter: { + 'dt_max': dt_max, + 'dt_slope_max': max_slope, + }, + } + custom_description['level_params'] = level_params + return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + + def get_custom_description_for_faults(self, problem, *args, **kwargs): + desc = self.get_custom_description(problem, *args, **kwargs) + if problem.__name__ == "run_quench": + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError + + desc['convergence_controllers'][AdaptivityPolynomialError]['e_tol'] = 1e-7 * 11 + desc['level_params']['dt'] = 4.0 + return desc + + def get_random_params(self, problem, num_procs): + ''' + Routine to get parameters for the randomization of faults + + Args: + problem: A function that runs a pySDC problem, see imports for available problems + num_procs (int): Number of processes you intend to run with + + Returns: + dict: Randomization parameters + ''' + + rnd_params = super().get_random_params(problem, num_procs) + if problem.__name__ == "run_quench": + rnd_params['iteration'] = 1 + elif problem.__name__ == 'run_Lorenz': + rnd_params['iteration'] = 4 + return rnd_params + + def get_reference_value(self, problem, key, op, num_procs=1): + """ + Get a reference value for a given problem for testing in CI. + + Args: + problem: A function that runs a pySDC problem, see imports for available problems + key (str): The name of the variable you want to compare + op (function): The operation you want to apply to the data + num_procs (int): Number of processes + + Returns: + The reference value + """ + if problem.__name__ == "run_vdp": + if key == 'work_newton' and op == sum: + return 2819 + elif key == 'e_global_post_run' and op == max: + return 1.8147451097405565e-07 + + raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') + + @property + def label(self): + return r'$\Delta t$-$k~\mathrm{adaptivity}$' diff --git a/pySDC/projects/Resilience/vdp.py b/pySDC/projects/Resilience/vdp.py index af2d8bbfe5..cae4e72ad3 100644 --- a/pySDC/projects/Resilience/vdp.py +++ b/pySDC/projects/Resilience/vdp.py @@ -107,7 +107,7 @@ def run_vdp( Returns: dict: The stats object controller: The controller - Tend: The time that was supposed to be integrated to + bool: If the code crashed """ # initialize level parameters @@ -182,13 +182,15 @@ def run_vdp( prepare_controller_for_faults(controller, fault_stuff, {}, {}) # call main function to get things done... + crash = False try: uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - except (ProblemError, ConvergenceError) as e: + except (ProblemError, ConvergenceError, ZeroDivisionError) as e: + crash = True print(f'Warning: Premature termination: {e}') stats = controller.return_stats() - return stats, controller, Tend + return stats, controller, crash def fetch_test_data(stats, comm=None, use_MPI=False): diff --git a/pySDC/projects/Resilience/work_precision.py b/pySDC/projects/Resilience/work_precision.py index 5b9681e1f1..7bde3d8463 100644 --- a/pySDC/projects/Resilience/work_precision.py +++ b/pySDC/projects/Resilience/work_precision.py @@ -17,6 +17,7 @@ setup_mpl(reset=True) LOGGER_LEVEL = 25 +LOG_TO_FILE = False logging.getLogger('matplotlib.texmanager').setLevel(90) @@ -28,11 +29,13 @@ 'k_SDC': ('k', sum, None), 'k_SDC_no_restart': ('k', sum, False), 'k_Newton': ('work_newton', sum, None), + 'k_linear': ('work_linear', sum, None), 'k_Newton_no_restart': ('work_newton', sum, False), 'k_rhs': ('work_rhs', sum, None), 'restart': ('restart', sum, None), 'dt_mean': ('dt', np.mean, False), 'dt_max': ('dt', max, False), + 'dt_min': ('dt', min, False), 'e_embedded_max': ('error_embedded_estimate', max, False), 'u0_increment_max': ('u0_increment', max, None), 'u0_increment_mean': ('u0_increment', np.mean, None), @@ -64,7 +67,16 @@ def get_forbidden_combinations(problem, strategy, **kwargs): def single_run( - problem, strategy, data, custom_description, num_procs=1, comm_world=None, problem_args=None, hooks=None, Tend=None + problem, + strategy, + data, + custom_description, + num_procs=1, + comm_world=None, + problem_args=None, + hooks=None, + Tend=None, + num_procs_sweeper=1, ): """ Make a single run of a particular problem with a certain strategy. @@ -77,6 +89,7 @@ def single_run( num_procs (int): Number of processes for the time communicator comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script hooks (list): List of additional hooks + num_procs_sweeper (int): Number of processes for the sweeper Returns: dict: Stats generated by the run @@ -89,24 +102,36 @@ def single_run( t_last = perf_counter() - comm = comm_world.Split(comm_world.rank < num_procs) - if comm_world.rank >= num_procs: + num_procs_tot = num_procs * num_procs_sweeper + comm = comm_world.Split(comm_world.rank < num_procs_tot) + if comm_world.rank >= num_procs_tot: comm.Free() return None + # make communicators for time and sweepers + comm_time = comm.Split(comm.rank // num_procs) + comm_sweep = comm.Split(comm_time.rank) + strategy_description = strategy.get_custom_description(problem, num_procs) description = merge_descriptions(strategy_description, custom_description) - - controller_params = {'logger_level': LOGGER_LEVEL} + if comm_sweep.size > 1: + description['sweeper_params']['comm'] = comm_sweep + + controller_params = { + 'logger_level': LOGGER_LEVEL, + 'log_to_file': LOG_TO_FILE, + 'fname': 'out.txt', + **strategy.get_controller_params(), + } problem_args = {} if problem_args is None else problem_args - stats, controller, _ = problem( + stats, controller, crash = problem( custom_description=description, Tend=strategy.get_Tend(problem, num_procs) if Tend is None else Tend, hook_class=[LogData, LogWork, LogGlobalErrorPostRun] + hooks, custom_controller_params=controller_params, use_MPI=True, - comm=comm, + comm=comm_time, **problem_args, ) @@ -115,8 +140,14 @@ def single_run( t_last = perf_counter() # record all the metrics + stats_all = filter_stats(stats, comm=comm_sweep) + comm_sweep.Free() + for key, mapping in MAPPINGS.items(): - me = get_sorted(stats, comm=comm, type=mapping[0], recomputed=mapping[2]) + if crash: + data[key] += [np.nan] + continue + me = get_sorted(stats_all, comm=comm_time, type=mapping[0], recomputed=mapping[2]) if len(me) == 0: data[key] += [np.nan] else: @@ -126,6 +157,7 @@ def single_run( logger.debug(f'Recorded all data after {t_now - t_last:.2e} s') t_last = perf_counter() + comm_time.Free() comm.Free() return stats @@ -165,7 +197,7 @@ def set_parameter(dictionary, where, parameter): set_parameter(dictionary[where[0]], where[1:], parameter) -def get_path(problem, strategy, num_procs, handle='', base_path='data/work_precision'): +def get_path(problem, strategy, num_procs, handle='', base_path='data/work_precision', num_procs_sweeper=1, mode=''): """ Get the path to a certain data. @@ -175,11 +207,13 @@ def get_path(problem, strategy, num_procs, handle='', base_path='data/work_preci num_procs (int): Number of processes for the time communicator handle (str): The name of the configuration base_path (str): Some path where all the files are stored + num_procs_sweeper (int): Number of processes for the sweeper + mode (str): The mode this was generated for Returns: str: The path to the data you are looking for """ - return f'{base_path}/{problem.__name__}-{strategy.__class__.__name__}-{handle}{"-wp" if handle else "wp"}-{num_procs}procs.pickle' + return f'{base_path}/{problem.__name__}-{strategy.__class__.__name__}-{handle}{"-wp" if handle else "wp"}-{num_procs}-{num_procs_sweeper}procs-{mode}.pickle' def record_work_precision( @@ -194,6 +228,8 @@ def record_work_precision( param_range=None, Tend=None, hooks=None, + num_procs_sweeper=1, + mode='', ): """ Run problem with strategy and record the cost parameters. @@ -206,6 +242,7 @@ def record_work_precision( handle (str): The name of the configuration runs (int): Number of runs you want to do comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script + num_procs_sweeper (int): Number of processes for the sweeper Returns: None @@ -224,7 +261,7 @@ def record_work_precision( if param == 'e_tol': power = 10.0 set_parameter(description, strategy.precision_parameter_loc[:-1] + ['dt_min'], 0) - exponents = [-3, -2, -1, 0, 1, 2, 3] + exponents = [-3, -2, -1, 0, 1, 2, 3][::-1] if problem.__name__ == 'run_vdp': exponents = [-4, -3, -2, -1, 0, 1, 2] elif param == 'dt': @@ -245,27 +282,30 @@ def record_work_precision( if problem.__name__ == 'run_quench': if param == 'restol': param_range = [1e-5, 1e-6, 1e-7, 1e-8, 1e-9] - elif param == 'e_tol': - param_range = [1e-2 / 2.0**me for me in [4, 5, 6, 7, 8, 9, 10]] elif param == 'dt': - param_range = [500 / 2.0**me for me in [5, 6, 7, 8]] + param_range = [1.25, 2.5, 5.0, 10.0, 20.0][::-1] + + elif problem.__name__ == 'run_AC': + if param == 'e_tol': + param_range = [1e-4, 1e-5, 1e-6, 1e-7, 1e-8][::-1] # run multiple times with different parameters for i in range(len(param_range)): - set_parameter(description, where, param_range[i]) - - if strategy.name == 'adaptivity_coll': - set_parameter(description, ['level_params', 'restol'], param_range[i] / 10.0) - - if problem.__name__ == 'run_Schroedinger' and strategy.name in ['ESDIRK', 'ERK']: - problem_params = custom_description['problem_params'] - if not problem_params.get('imex', True): - description["problem_params"]['lintol'] = param_range[i] + set_parameter(description, where, param_range[i] * strategy.precision_range_fac) data[param_range[i]] = {key: [] for key in MAPPINGS.keys()} data[param_range[i]]['param'] = [param_range[i]] data[param_range[i]][param] = [param_range[i]] + + description = merge_descriptions( + descA=description, descB=strategy.get_description_for_tolerance(problem=problem, param=param_range[i]) + ) for _j in range(runs): + if comm_world.rank == 0: + logger.log( + 24, + f'Starting: {problem.__name__}: {strategy} {handle} {num_procs}-{num_procs_sweeper} procs, {param}={param_range[i]:.2e}', + ) single_run( problem, strategy, @@ -276,14 +316,21 @@ def record_work_precision( num_procs=num_procs, hooks=hooks, Tend=Tend, + num_procs_sweeper=num_procs_sweeper, ) comm_world.Barrier() if comm_world.rank == 0: + if np.isfinite(data[param_range[i]]["k_linear"][-1]): + k_type = "k_linear" + elif np.isfinite(data[param_range[i]]["k_Newton"][-1]): + k_type = 'k_Newton' + else: + k_type = "k_SDC" logger.log( 25, - f'{problem.__name__}: {strategy} {handle} {num_procs} procs, {param}={param_range[i]:.2e}: e={data[param_range[i]]["e_global"][-1]}, t={data[param_range[i]]["t"][-1]}, k={data[param_range[i]]["k_SDC"][-1]}', + f'{problem.__name__}: {strategy} {handle} {num_procs}-{num_procs_sweeper} procs, {param}={param_range[i]:.2e}: e={data[param_range[i]]["e_global"][-1]}, t={data[param_range[i]]["t"][-1]}, {k_type}={data[param_range[i]][k_type][-1]}', ) if comm_world.rank == 0: @@ -295,9 +342,11 @@ def record_work_precision( 'time': time.time, 'runs': runs, } - with open(get_path(problem, strategy, num_procs, handle), 'wb') as f: - logger.debug(f'Dumping file \"{get_path(problem, strategy, num_procs, handle=handle)}\"') + path = get_path(problem, strategy, num_procs, handle, num_procs_sweeper=num_procs_sweeper, mode=mode) + with open(path, 'wb') as f: + logger.debug(f'Dumping file \"{path}\"') pickle.dump(data, f) + return data def load(**kwargs): @@ -352,6 +401,8 @@ def plot_work_precision( handle='', plotting_params=None, comm_world=None, + num_procs_sweeper=1, + mode='', ): # pragma: no cover """ Plot data from running a problem with a strategy. @@ -366,6 +417,8 @@ def plot_work_precision( handle (str): The name of the configuration plotting_params (dict): Will be passed when plotting comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script + num_procs_sweeper (int): Number of processes for the sweeper + mode (str): The of the configurations you want to retrieve Returns: None @@ -373,7 +426,14 @@ def plot_work_precision( if comm_world.rank > 0 or get_forbidden_combinations(problem, strategy): return None - data = load(problem=problem, strategy=strategy, num_procs=num_procs, handle=handle) + data = load( + problem=problem, + strategy=strategy, + num_procs=num_procs, + handle=handle, + num_procs_sweeper=num_procs_sweeper, + mode=mode, + ) work, precision = extract_data(data, work_key, precision_key) keys = [key for key in data.keys() if key not in ['meta']] @@ -382,7 +442,7 @@ def plot_work_precision( rel_variance = [np.std(data[me][key]) / max([np.nanmean(data[me][key]), 1.0]) for me in keys] if not all(me < 1e-1 or not np.isfinite(me) for me in rel_variance): logger.warning( - f"Variance in \"{key}\" for {get_path(problem, strategy, num_procs, handle)} too large! Got {rel_variance}" + f"Variance in \"{key}\" for {get_path(problem, strategy, num_procs, handle, num_procs_sweeper=num_procs_sweeper, mode=mode)} too large! Got {rel_variance}" ) style = merge_descriptions( @@ -410,6 +470,46 @@ def plot_work_precision( ax.text(0.1, 0.2, "No sampling!", transform=ax.transAxes) +def plot_parallel_efficiency_diagonalSDC( + ax, work_key, precision_key, num_procs_sweeper, num_procs=1, **kwargs +): # pragma: no cover + serial_data = load( + num_procs=num_procs, + num_procs_sweeper=1, + **kwargs, + ) + parallel_data = load( + num_procs=num_procs, + num_procs_sweeper=num_procs_sweeper, + **kwargs, + ) + serial_work, serial_precision = extract_data(serial_data, work_key, precision_key) + parallel_work, parallel_precision = extract_data(parallel_data, work_key, precision_key) + # assert np.allclose(serial_precision, parallel_precision) + + serial_work = np.asarray(serial_work) + parallel_work = np.asarray(parallel_work) + + # ax.loglog(serial_work, serial_precision) + # ax.loglog(parallel_work, parallel_precision) + + speedup = serial_work / parallel_work + parallel_efficiency = np.median(speedup) / num_procs_sweeper + ax.plot(serial_precision, speedup) + ax.set_xscale('log') + ax.set_ylabel('speedup') + + if 't' in [work_key, precision_key]: + meta = parallel_data.get('meta', {}) + + if meta.get('hostname', None) in ['thomas-work']: + ax.text(0.1, 0.1, "Laptop timings!", transform=ax.transAxes) + if meta.get('runs', None) == 1: + ax.text(0.1, 0.2, "No sampling!", transform=ax.transAxes) + + return np.median(speedup), parallel_efficiency + + def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only=False): # pragma: no cover """ Decorate a plot @@ -443,6 +543,10 @@ def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only 'u0_increment_mean': r'$\bar{\Delta u_0}$', 'u0_increment_max_no_restart': r'$\| \Delta u_0 \|_{\infty} $ (restarts excluded)', 'u0_increment_mean_no_restart': r'$\bar{\Delta u_0}$ (restarts excluded)', + 'k_linear': 'Linear solver iterations', + 'speedup': 'Speedup', + 'nprocs': r'$N_\mathrm{procs}$', + '': '', } if not title_only: @@ -472,6 +576,8 @@ def execute_configurations( comm_world, plotting, Tend=None, + num_procs_sweeper=1, + mode='', ): """ Run for multiple configurations. @@ -488,6 +594,8 @@ def execute_configurations( runs (int): Number of runs you want to do comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script plotting (bool): Whether to plot something + num_procs_sweeper (int): Number of processes for the sweeper + mode (str): What you want to look at Returns: None @@ -499,6 +607,7 @@ def execute_configurations( 'strategy': strategy, 'handle': config.get('handle', ''), 'num_procs': config.get('num_procs', num_procs), + 'num_procs_sweeper': config.get('num_procs_sweeper', num_procs_sweeper), } if record: logger.debug('Recording work precision') @@ -511,6 +620,7 @@ def execute_configurations( param_range=config.get('param_range', None), hooks=config.get('hooks', None), Tend=Tend, + mode=mode, ) if plotting and comm_world.rank == 0: logger.debug('Plotting') @@ -521,6 +631,7 @@ def execute_configurations( ax=ax, plotting_params=config.get('plotting_params', {}), comm_world=comm_world, + mode=mode, ) if comm_world.rank == 0: @@ -558,7 +669,10 @@ def get_configs(mode, problem): from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter from pySDC.projects.Resilience.strategies import AdaptivityStrategy, ESDIRKStrategy - limits = [25.0, 50.0] + limits = [ + 25.0, + 50.0, + ] colors = ['teal', 'magenta'] markers = ['v', 'x'] markersize = 9 @@ -571,7 +685,7 @@ def get_configs(mode, problem): 'color': colors[i], 'marker': markers[i], 'label': rf'$\Delta t \leq {{{limits[i]:.0f}}}$', - 'ls': '', + # 'ls': '', 'markersize': markersize, }, 'num_procs': 1, @@ -581,31 +695,22 @@ def get_configs(mode, problem): 'handle': 'no limits', 'plotting_params': { 'label': 'no limiter', - 'ls': '', + # 'ls': '', 'markersize': markersize, }, 'strategies': [AdaptivityStrategy(useMPI=True)], 'num_procs': 1, } elif mode == 'dynamic_restarts': + """ + Compare Block Gauss-Seidel SDC with restarting the first step in the block or the first step that exceeded the error threshold. + """ from pySDC.projects.Resilience.strategies import AdaptivityStrategy, AdaptivityRestartFirstStep desc = {} desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'} desc['step_params'] = {'maxiter': 5} - # configurations[1] = { - # 'custom_description': {}, - # 'handle': 'dynamic restarts', - # 'plotting_params': {'label': 'adaptivity (dynamic restarts)', 'ls': '-.'}, - # 'strategies': [AdaptivityStrategy(useMPI=True)], - # } - # configurations[2] = { - # 'custom_description': {}, - # 'handle': 'regular', - # 'plotting_params': {'label': 'adaptivity'}, - # 'strategies': [AdaptivityRestartFirstStep(useMPI=True)], - # } ls = { 1: '-', 2: '--', @@ -629,116 +734,139 @@ def get_configs(mode, problem): } elif mode == 'compare_strategies': - from pySDC.projects.Resilience.strategies import AdaptivityStrategy, BaseStrategy, IterateStrategy - - description_high_order = {'step_params': {'maxiter': 5}} - description_low_order = {'step_params': {'maxiter': 3}} - dashed = {'ls': '--'} + """ + Compare the different SDC strategies. + """ + from pySDC.projects.Resilience.strategies import ( + AdaptivityStrategy, + kAdaptivityStrategy, + AdaptivityPolynomialError, + BaseStrategy, + ) - configurations[0] = { - 'custom_description': description_high_order, - 'handle': r'high order', - 'strategies': [AdaptivityStrategy(useMPI=True), BaseStrategy(useMPI=True)], - } configurations[1] = { - 'custom_description': description_low_order, - 'handle': r'low order', - 'strategies': [AdaptivityStrategy(useMPI=True), BaseStrategy(useMPI=True)], - 'plotting_params': dashed, + 'strategies': [AdaptivityPolynomialError(useMPI=True)], } - - description_large_step = {'level_params': {'dt': 5.0 if problem.__name__ == 'run_quench' else 3e-2}} - description_small_step = {'level_params': {'dt': 1.0 if problem.__name__ == 'run_quench' else 1e-2}} - configurations[2] = { - 'custom_description': description_large_step, - 'handle': r'large step', - 'strategies': [IterateStrategy(useMPI=True)], - 'plotting_params': dashed, + 'strategies': [kAdaptivityStrategy(useMPI=True)], } - configurations[3] = { - 'custom_description': description_small_step, - 'handle': r'small step', - 'strategies': [IterateStrategy(useMPI=True)], - } - elif mode == 'RK': - from pySDC.projects.Resilience.strategies import AdaptivityStrategy, DIRKStrategy, ERKStrategy - - # from pySDC.implementations.sweeper_classes.explicit import explicit - # configurations[3] = { - # 'custom_description': { - # 'step_params': {'maxiter': 5}, - # 'sweeper_params': {'QE': 'EE'}, - # 'sweeper_class': explicit, - # }, - # 'handle': 'explicit order 4', - # 'strategies': [AdaptivityStrategy(useMPI=True)], - # 'plotting_params': {'ls': ':', 'label': 'explicit SDC5(4)'}, - # } + configurations[0] = { - 'strategies': [ERKStrategy(useMPI=True), DIRKStrategy(useMPI=True)], - 'num_procs': 1, - } - configurations[1] = { - 'custom_description': {'step_params': {'maxiter': 5}}, - 'handle': 'order 5', - 'strategies': [AdaptivityStrategy(useMPI=True)], - 'plotting_params': {'label': 'SDC5(4)'}, - } - configurations[2] = { - 'custom_description': {'step_params': {'maxiter': 4}}, - 'handle': 'order 4', - 'strategies': [AdaptivityStrategy(useMPI=True)], - 'plotting_params': {'ls': '--', 'label': 'SDC4(3)'}, + 'custom_description': { + 'step_params': {'maxiter': 5}, + 'sweeper_params': {'num_nodes': 3, 'quad_type': 'RADAU-RIGHT'}, + }, + 'strategies': [AdaptivityStrategy(useMPI=True), BaseStrategy(useMPI=True)], } + + elif mode == 'interpolate_between_restarts': + """ + Compare adaptivity with interpolation between restarts and without + """ + from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError + + i = 0 + for interpolate_between_restarts, handle, ls in zip( + [True, False], ['Interpolation bewteen restarts', 'regular'], ['--', '-'] + ): + configurations[i] = { + 'strategies': [ + AdaptivityPolynomialError(interpolate_between_restarts=interpolate_between_restarts, useMPI=True) + ], + 'plotting_params': {'ls': ls}, + 'handle': handle, + } + i += 1 + elif mode == 'diagonal_SDC': + """ + Run diagonal SDC with different number of nodes and ranks. You can use this to compute a speedup, but it's not strong scaling! + """ + from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError + + if problem.__name__ in ['run_Schroedinger']: + from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper + else: + from pySDC.implementations.sweeper_classes.generic_implicit_MPI import ( + generic_implicit_MPI as parallel_sweeper, + ) + + for parallel in [False, True]: + desc = {'sweeper_class': parallel_sweeper} if parallel else {} + for num_nodes, ls in zip([3, 4, 2], ['-', '--', ':', '-.']): + configurations[num_nodes + (99 if parallel else 0)] = { + 'custom_description': {**desc, 'sweeper_params': {'num_nodes': num_nodes}}, + 'strategies': [ + AdaptivityPolynomialError(useMPI=True, newton_inexactness=True, linear_inexactness=True) + ], + 'num_procs_sweeper': num_nodes if parallel else 1, + 'num_procs': 1, + 'handle': f'{num_nodes} nodes', + 'plotting_params': { + 'ls': ls, + 'label': f'{num_nodes} procs', + # **{'color': 'grey' if parallel else None}, + }, + } + elif mode == 'parallel_efficiency': - from pySDC.projects.Resilience.strategies import ( - AdaptivityStrategy, - BaseStrategy, - IterateStrategy, - ) + """ + Compare parallel runs of the step size adaptive SDC + """ + from pySDC.projects.Resilience.strategies import AdaptivityStrategy, AdaptivityPolynomialError + + if problem.__name__ in ['run_Schroedinger']: + from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper + else: + from pySDC.implementations.sweeper_classes.generic_implicit_MPI import ( + generic_implicit_MPI as parallel_sweeper, + ) desc = {} - desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'} + desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': 'EE'} desc['step_params'] = {'maxiter': 5} - descIterate = {} - descIterate['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'} - ls = { 1: '-', 2: '--', 3: '-.', - 4: ':', + 4: '--', 5: ':', + 12: ':', } - for num_procs in [4, 2, 1]: - plotting_params = {'ls': ls[num_procs], 'label': fr'$\Delta t$ adaptivity {num_procs} procs'} + for num_procs in [4, 1]: + plotting_params = ( + {'ls': ls[num_procs], 'label': fr'$\Delta t$ adaptivity $N$={num_procs}x1'} if num_procs > 1 else {} + ) configurations[num_procs] = { - 'strategies': [AdaptivityStrategy(True)], - 'custom_description': desc, + 'strategies': [AdaptivityStrategy(useMPI=True)], + 'custom_description': desc.copy(), 'num_procs': num_procs, - 'plotting_params': plotting_params, + 'plotting_params': plotting_params.copy(), } - plotting_params = {'ls': ls[num_procs], 'label': fr'$k$ adaptivity {num_procs} procs'} - configurations[num_procs + 100] = { - 'strategies': [IterateStrategy(True)], - 'custom_description': descIterate, + configurations[num_procs * 100 + 79] = { + 'custom_description': {'sweeper_class': parallel_sweeper}, + 'strategies': [ + AdaptivityPolynomialError(useMPI=True, newton_inexactness=True, linear_inexactness=True) + ], + 'num_procs_sweeper': 3, 'num_procs': num_procs, - 'plotting_params': plotting_params, + 'plotting_params': { + 'ls': ls.get(num_procs * 3, '-'), + 'label': rf'$\Delta t$-$k$ adaptivity $N$={num_procs}x3', + }, } - # plotting_params = {'ls': ls[num_procs], 'label': fr'fixed {num_procs} procs'} - # configurations[num_procs + 999] = { - # 'strategies': [BaseStrategy(True)], - # 'custom_description': desc, - # 'num_procs': num_procs, - # 'plotting_params': plotting_params, - # } + configurations[num_procs * 200 + 79] = { + 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=True, linear_inexactness=True)], + 'num_procs': 1, + } elif mode[:13] == 'vdp_stiffness': - from pySDC.projects.Resilience.strategies import AdaptivityStrategy, ERKStrategy, DIRKStrategy, ESDIRKStrategy + """ + Run van der Pol with different parameter for the nonlinear term, which controls the stiffness. + """ + from pySDC.projects.Resilience.strategies import AdaptivityStrategy, ERKStrategy, ESDIRKStrategy mu = float(mode[14:]) @@ -781,15 +909,7 @@ def get_configs(mode, problem): 'handle': mode, 'plotting_params': {'label': 'CP5(4)'}, 'custom_description': problem_desc, - #'param_range': [1e-2], } - # configurations[3] = { - # 'strategies': [DIRKStrategy(useMPI=True)], - # 'num_procs': 1, - # 'handle': mode, - # 'plotting_params': {'label': 'DIRK4(3)'}, - # 'custom_description': problem_desc, - # } configurations[4] = { 'strategies': [ESDIRKStrategy(useMPI=True)], 'num_procs': 1, @@ -797,130 +917,179 @@ def get_configs(mode, problem): 'plotting_params': {'label': 'ESDIRK5(3)'}, 'custom_description': problem_desc, } - - elif mode == 'compare_adaptivity': - # TODO: configurations not final! + elif mode == 'inexactness': + """ + Compare inexact SDC to exact SDC + """ from pySDC.projects.Resilience.strategies import ( - AdaptivityCollocationTypeStrategy, - AdaptivityCollocationRefinementStrategy, - AdaptivityStrategy, - AdaptivityExtrapolationWithinQStrategy, + AdaptivityPolynomialError, ) + if problem.__name__ in ['run_Schroedinger']: + from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper + else: + from pySDC.implementations.sweeper_classes.generic_implicit_MPI import ( + generic_implicit_MPI as parallel_sweeper, + ) + strategies = [ - AdaptivityCollocationTypeStrategy(useMPI=True), - AdaptivityCollocationRefinementStrategy(useMPI=True), + AdaptivityPolynomialError, ] - restol = None - for strategy in strategies: - strategy.restol = restol + inexactness = { + 'newton_inexactness': True, + 'linear_inexactness': True, + } + no_inexactness = { + 'newton_inexactness': False, + 'linear_inexactness': False, + 'SDC_maxiter': 99, + 'use_restol_rel': False, + } configurations[1] = { - 'custom_description': {'step_params': {'maxiter': 99}, 'level_params': {'restol': 1e-11}}, - 'strategies': [AdaptivityExtrapolationWithinQStrategy(useMPI=True)], + 'custom_description': {'sweeper_class': parallel_sweeper}, + 'strategies': [me(useMPI=True, **no_inexactness) for me in strategies], + 'num_procs_sweeper': 3, + 'handle': 'exact', + 'plotting_params': {'ls': '--'}, } - configurations[2] = {'strategies': strategies} - configurations[3] = { - 'custom_description': {'step_params': {'maxiter': 5}}, - 'strategies': [AdaptivityStrategy(useMPI=True)], + configurations[0] = { + 'custom_description': {'sweeper_class': parallel_sweeper}, + 'strategies': [me(useMPI=True, **inexactness) for me in strategies], + 'handle': 'inexact', + 'num_procs_sweeper': 3, } - - # strategies2 = [AdaptivityCollocationTypeStrategy(useMPI=True), AdaptivityCollocationRefinementStrategy(useMPI=True)] - # restol = 1e-6 - # for strategy in strategies2: - # strategy.restol = restol - # configurations[3] = {'strategies':strategies2, 'handle': 'low restol', 'plotting_params': {'ls': '--'}} - - elif mode == 'quench': + elif mode == 'compare_adaptivity': + """ + Compare various modes of adaptivity + """ + # TODO: configurations not final! from pySDC.projects.Resilience.strategies import ( + # AdaptivityCollocationTypeStrategy, + # AdaptivityCollocationRefinementStrategy, AdaptivityStrategy, - DoubleAdaptivityStrategy, - IterateStrategy, - BaseStrategy, + # AdaptivityExtrapolationWithinQStrategy, + ESDIRKStrategy, + ARKStrategy, + AdaptivityPolynomialError, ) - dumbledoresarmy = DoubleAdaptivityStrategy(useMPI=True) - # dumbledoresarmy.residual_e_tol_ratio = 1e2 - dumbledoresarmy.residual_e_tol_abs = 1e-3 + if problem.__name__ in ['run_Schroedinger']: + from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper + else: + from pySDC.implementations.sweeper_classes.generic_implicit_MPI import ( + generic_implicit_MPI as parallel_sweeper, + ) + + inexactness_params = { + # 'double_adaptivity': True, + 'newton_inexactness': True, + 'linear_inexactness': True, + } strategies = [ - AdaptivityStrategy(useMPI=True), - IterateStrategy(useMPI=True), - BaseStrategy(useMPI=True), - dumbledoresarmy, + AdaptivityPolynomialError, + # AdaptivityCollocationTypeStrategy, + # AdaptivityExtrapolationWithinQStrategy, ] - configurations[1] = {'strategies': strategies} + + # restol = None + # for strategy in strategies: + # strategy.restol = restol + + configurations[1] = { + 'custom_description': {'sweeper_class': parallel_sweeper}, + 'strategies': [me(useMPI=True, **inexactness_params) for me in strategies], + 'handle': 'parallel', + 'num_procs_sweeper': 3, + 'plotting_params': {'ls': '-', 'label': '3 procs'}, + } configurations[2] = { - 'strategies': strategies, - 'problem_args': {'imex': True}, - 'handle': 'IMEX', + 'strategies': [me(useMPI=True, **inexactness_params) for me in strategies], 'plotting_params': {'ls': '--'}, } - inexact = {'problem_params': {'newton_iter': 30}} - configurations[3] = { - 'strategies': strategies, - 'custom_description': inexact, - 'handle': 'inexact', - 'plotting_params': {'ls': ':'}, - } - LU = {'sweeper_params': {'QI': 'LU'}} configurations[4] = { - 'strategies': strategies, - 'custom_description': LU, - 'handle': 'LU', - 'plotting_params': {'ls': '-.'}, + 'custom_description': {'step_params': {'maxiter': 5}}, + 'strategies': [AdaptivityStrategy(useMPI=True)], } + + desc_RK = {} + configurations[-1] = { + 'strategies': [ + ARKStrategy(useMPI=True) if problem.__name__ == 'run_Schroedinger' else ESDIRKStrategy(useMPI=True), + ], + 'num_procs': 1, + 'custom_description': desc_RK, + } + elif mode == 'preconditioners': + """ + Compare different preconditioners + """ from pySDC.projects.Resilience.strategies import ( AdaptivityStrategy, IterateStrategy, BaseStrategy, ESDIRKStrategy, ERKStrategy, + AdaptivityPolynomialError, ) - strategies = [AdaptivityStrategy(useMPI=True)] + inexacness = True + strategies = [ + AdaptivityPolynomialError( + useMPI=True, SDC_maxiter=29, newton_inexactness=inexacness, linear_inexactness=inexacness + ), + BaseStrategy(useMPI=True), + ] desc = {} - desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'} - desc['step_params'] = {'maxiter': 5} + desc['sweeper_params'] = { + 'num_nodes': 3, + } + # desc['step_params'] = {'maxiter': 5} precons = ['IE', 'LU'] - ls = ['-', '--', '-.', ':'] - for i in range(len(precons)): - desc['sweeper_params']['QI'] = precons[i] + ls = ['-.', '--', '-', ':'] + for i in range(len(precons) + 1): + if i < len(precons): + desc['sweeper_params']['QI'] = precons[i] + handle = precons[i] + else: + handle = None configurations[i] = { 'strategies': strategies, 'custom_description': copy.deepcopy(desc), - 'handle': precons[i], + 'handle': handle, 'plotting_params': {'ls': ls[i]}, } - configurations[-1] = { - 'strategies': [ - ERKStrategy(useMPI=True), - ESDIRKStrategy(useMPI=True), - ], - 'num_procs': 1, - } - elif mode == 'RK_comp': + """ + Compare parallel adaptive SDC to Runge-Kutta + """ from pySDC.projects.Resilience.strategies import ( AdaptivityStrategy, - BaseStrategy, - IterateStrategy, ERKStrategy, - DIRKStrategy, ESDIRKStrategy, + ARKStrategy, + AdaptivityPolynomialError, ) + if problem.__name__ in ['run_Schroedinger']: + from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper + else: + from pySDC.implementations.sweeper_classes.generic_implicit_MPI import ( + generic_implicit_MPI as parallel_sweeper, + ) + desc = {} - desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'} + desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': "EE"} desc['step_params'] = {'maxiter': 5} - descIterate = {} - descIterate['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'} + desc_poly = {} + desc_poly['sweeper_class'] = parallel_sweeper ls = { 1: '-', @@ -932,30 +1101,56 @@ def get_configs(mode, problem): desc_RK = {} if problem.__name__ in ['run_Schroedinger']: - desc_RK['problem_params'] = {'imex': False} + desc_RK['problem_params'] = {'imex': True} + configurations[3] = { + 'custom_description': desc_poly, + 'strategies': [AdaptivityPolynomialError(useMPI=True)], + 'num_procs': 1, + 'num_procs_sweeper': 3, + 'plotting_params': {'label': r'$\Delta t$-$k$ adaptivity $N$=1x3'}, + } configurations[-1] = { 'strategies': [ ERKStrategy(useMPI=True), - ESDIRKStrategy(useMPI=True), + ARKStrategy(useMPI=True) if problem.__name__ in ['run_Schroedinger'] else ESDIRKStrategy(useMPI=True), ], 'num_procs': 1, 'custom_description': desc_RK, } - for num_procs, label in zip([4, 1], ['GSSDC 4 procs', 'SDC']): - configurations[num_procs] = { - 'strategies': [AdaptivityStrategy(True)], - 'custom_description': desc, - 'num_procs': num_procs, - 'plotting_params': {'ls': ls[num_procs], 'label': label}, - } - elif mode == 'imex': - from pySDC.projects.Resilience.strategies import AdaptivityStrategy, ESDIRKStrategy + configurations[2] = { + 'strategies': [AdaptivityStrategy(useMPI=True)], + 'custom_description': desc, + 'num_procs': 4, + 'plotting_params': {'label': r'$\Delta t$ adaptivity $N$=4x1'}, + } + elif mode == 'RK_comp_high_order': + """ + Compare higher order SDC than we can get with RKM to RKM + """ + from pySDC.projects.Resilience.strategies import ( + AdaptivityStrategy, + ERKStrategy, + ESDIRKStrategy, + ARKStrategy, + AdaptivityPolynomialError, + ) + + if problem.__name__ in ['run_Schroedinger']: + from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper + else: + from pySDC.implementations.sweeper_classes.generic_implicit_MPI import ( + generic_implicit_MPI as parallel_sweeper, + ) desc = {} - desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'} - desc['step_params'] = {'maxiter': 5} + desc['sweeper_params'] = {'num_nodes': 4, 'QI': 'IE', 'QE': "EE"} + desc['step_params'] = {'maxiter': 7} + + desc_poly = {} + desc_poly['sweeper_class'] = parallel_sweeper + desc_poly['sweeper_params'] = {'num_nodes': 4} ls = { 1: '-', @@ -966,60 +1161,43 @@ def get_configs(mode, problem): } desc_RK = {} - desc_RK['problem_params'] = {'imex': False} + if problem.__name__ in ['run_Schroedinger']: + desc_RK['problem_params'] = {'imex': True} + + configurations[3] = { + 'custom_description': desc_poly, + 'strategies': [AdaptivityPolynomialError(useMPI=True)], + 'num_procs': 1, + 'num_procs_sweeper': 4, + } configurations[-1] = { 'strategies': [ - ESDIRKStrategy(useMPI=True), + ERKStrategy(useMPI=True), + ARKStrategy(useMPI=True) if problem.__name__ in ['run_Schroedinger'] else ESDIRKStrategy(useMPI=True), ], 'num_procs': 1, 'custom_description': desc_RK, } - i = 10 - for imex, label, ls in zip([True, False], ['IMEX', 'fully implicit'], ['-', ':']): - desc['problem_params'] = {'imex': imex} - plotting_params = {'ls': ls, 'label': f'{label} SDC'} - configurations[i] = { - 'strategies': [AdaptivityStrategy(True)], - 'custom_description': copy.deepcopy(desc), - 'num_procs': 1, - 'plotting_params': plotting_params, - 'handle': label, - } - i += 1 - elif mode == 'newton_tol': - from pySDC.projects.Resilience.strategies import AdaptivityStrategy, BaseStrategy, IterateStrategy - tol_range = [1e-7, 1e-9, 1e-11] - ls = ['-', '--', '-.', ':'] - for i in range(len(tol_range)): - configurations[i] = { - 'strategies': [AdaptivityStrategy(useMPI=True), BaseStrategy(useMPI=True)], - 'custom_description': { - 'problem_params': {'newton_tol': tol_range[i]}, - 'step_params': {'maxiter': 5}, - }, - 'handle': f"Newton tol={tol_range[i]:.1e}", - 'plotting_params': {'ls': ls[i]}, - } - configurations[i + len(tol_range)] = { - 'strategies': [IterateStrategy(useMPI=True)], - 'custom_description': { - 'problem_params': {'newton_tol': tol_range[i]}, - }, - 'handle': f"Newton tol={tol_range[i]:.1e}", - 'plotting_params': {'ls': ls[i]}, - } + configurations[2] = { + 'strategies': [AdaptivityStrategy(useMPI=True)], + 'custom_description': desc, + 'num_procs': 4, + } elif mode == 'avoid_restarts': + """ + Test how well avoiding restarts works. + """ from pySDC.projects.Resilience.strategies import ( AdaptivityStrategy, AdaptivityAvoidRestartsStrategy, - AdaptivityInterpolationStrategy, + AdaptivityPolynomialStrategy, ) desc = {'sweeper_params': {'QI': 'IE'}, 'step_params': {'maxiter': 3}} param_range = [1e-3, 1e-5] configurations[0] = { - 'strategies': [AdaptivityInterpolationStrategy(useMPI=True)], + 'strategies': [AdaptivityPolynomialStrategy(useMPI=True)], 'plotting_params': {'ls': '--'}, 'custom_description': desc, 'param_range': param_range, @@ -1063,7 +1241,7 @@ def get_fig(x=1, y=1, **kwargs): # pragma: no cover def save_fig( - fig, name, work_key, precision_key, legend=True, format='pdf', base_path='data', **kwargs + fig, name, work_key, precision_key, legend=True, format='pdf', base_path='data', squares=True, ncols=None, **kwargs ): # pragma: no cover """ Save a figure with a legend on the bottom. @@ -1075,17 +1253,26 @@ def save_fig( precision_key (str): The key in the recorded data you want on the y-axis legend (bool): Put a legend or not format (str): Format to store the figure with + base_path (str): Some path where all the files are stored + squares (bool): Adjust aspect ratio to squares if true Returns: None """ - handles, labels = fig.get_axes()[0].get_legend_handles_labels() + handles = [] + labels = [] + for ax in fig.get_axes(): + h, l = ax.get_legend_handles_labels() + handles += [h[i] for i in range(len(h)) if l[i] not in labels] + labels += [me for me in l if me not in labels] + if squares: + ax.set_box_aspect(1) order = np.argsort([me[0] for me in labels]) fig.legend( [handles[i] for i in order], [labels[i] for i in order], loc='outside lower center', - ncols=3 if len(handles) % 3 == 0 else 4, + ncols=ncols if ncols else 3 if len(handles) % 3 == 0 else 4, frameon=False, fancybox=True, ) @@ -1130,6 +1317,7 @@ def all_problems(mode='compare_strategies', plotting=True, base_path='data', **k ax=axs.flatten()[i], decorate=True, configurations=get_configs(mode, problems[i]), + mode=mode, ) if plotting and shared_params['comm_world'].rank == 0: @@ -1140,6 +1328,7 @@ def all_problems(mode='compare_strategies', plotting=True, base_path='data', **k precision_key=shared_params['precision_key'], legend=True, base_path=base_path, + ncols=3 if mode in ['parallel_efficiency'] else None, ) @@ -1198,7 +1387,7 @@ def single_problem(mode, problem, plotting=True, base_path='data', **kwargs): # mode (str): What you want to look at problem (function): A problem to run """ - fig, ax = get_fig(1, 1, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1, 0.5)) + fig, ax = get_fig(1, 1, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1, 0.8)) params = { 'work_key': 'k_SDC', @@ -1212,7 +1401,9 @@ def single_problem(mode, problem, plotting=True, base_path='data', **kwargs): # } logger.log(26, f"Doing single problem {mode}") - execute_configurations(**params, problem=problem, ax=ax, decorate=True, configurations=get_configs(mode, problem)) + execute_configurations( + **params, problem=problem, ax=ax, decorate=True, configurations=get_configs(mode, problem), mode=mode + ) if plotting: save_fig( @@ -1261,117 +1452,89 @@ def vdp_stiffness_plot(base_path='data', format='pdf', **kwargs): # pragma: no ) -def ERK_stiff_weirdness(): - from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostStep, LogLocalErrorPostStep - from pySDC.projects.Resilience.strategies import ERKStrategy, AdaptivityStrategy - from pySDC.implementations.convergence_controller_classes.step_size_limiter import ( - StepSizeLimiter, - StepSizeSlopeLimiter, - ) - import matplotlib.pyplot as plt - from pySDC.projects.Resilience.Schroedinger import live_plotting, log_number_of_waves +def aggregate_parallel_efficiency_plot(): # pragma: no cover + """ + Make a "weak scaling" plot for diagonal SDC + """ + from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError - fig, axs = plt.subplots(4, 1, sharex=True) + fig, axs = plt.subplots(2, 2) - strategy = AdaptivityStrategy(True) - # strategy = ERKStrategy(True) + _fig, _ax = plt.subplots(1, 1) + num_procs = 1 + num_procs_sweeper = 2 + problem = run_quench - problem = run_Schroedinger - # epsilons = [1e-1, 1e-3, 1e-5, 1e-7, 1e-8] - epsilons = [1e-2, 1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8] - # epsilons = [1e-2, 1e-8] - Tend = 0.1 + num_procs_sweeper_list = [2, 3, 4] - data = {} - desc = strategy.get_custom_description(problem, 1) - where = strategy.precision_parameter_loc - # desc['problem_params'] = {'mu': 30} - # desc['problem_params'] = {'imex': False} - desc['level_params'] = {'dt': Tend * 1e-3} - desc['sweeper_params']['QI'] = 'EE' - # desc['convergence_controllers'][StepSizeLimiter]= {'dt_slope_max': 10., 'dt_max': 5e-4} - - from pySDC.projects.Resilience.strategies import cmap - - for eps, color in zip(epsilons, list(cmap.values())): - set_parameter(desc, where, eps) - data[eps] = {key: [] for key in MAPPINGS.keys()} - data[eps]['stats'] = single_run( - strategy=strategy, - data=data[eps], - hooks=[LogLocalErrorPostStep, log_number_of_waves], - problem=problem, - custom_description=desc, - num_procs=1, - comm_world=MPI.COMM_WORLD, - Tend=Tend, - ) - stats = data[eps]['stats'] - # e = get_sorted(stats, type='e_global_post_step', recomputed=False) - # axs[0].plot([me[0] for me in e], [me[1] for me in e]) - # - e = get_sorted(stats, type='e_local_post_step', recomputed=False) - axs[0].plot( - [me[0] for me in e], [me[1] for me in e], ls='-', label=rf'$\epsilon_\mathrm{{TOL}} = {{{eps:.0e}}}$' + for problem, ax in zip([run_vdp, run_Lorenz, run_quench], axs.flatten()): + speedup = [] + for num_procs_sweeper in num_procs_sweeper_list: + s, e = plot_parallel_efficiency_diagonalSDC( + ax=_ax, + work_key='t', + precision_key='e_global_rel', + num_procs=num_procs, + num_procs_sweeper=num_procs_sweeper, + problem=problem, + strategy=AdaptivityPolynomialError(), + mode='diagonal_SDC', + handle=f'{num_procs_sweeper} nodes', + ) + speedup += [s] + decorate_panel(ax, problem, work_key='nprocs', precision_key='') + + ax.plot(num_procs_sweeper_list, speedup, label='speedup') + ax.plot( + num_procs_sweeper_list, + [speedup[i] / num_procs_sweeper_list[i] for i in range(len(speedup))], + label='parallel efficiency', ) - axs[0].set_ylabel('local error') - - dt = get_sorted(stats, type='dt', recomputed=False) - axs[1].plot([me[0] for me in dt], [me[1] for me in dt], ls='-') - axs[1].set_ylabel('dt') - axs[1].set_yscale('log') - - rhs = get_sorted(stats, type='work_rhs', recomputed=None) - - axs[2].plot([me[0] for me in rhs], np.cumsum([me[1] for me in rhs]), ls='-') - axs[2].set_ylabel('rhs evaluations') - nwaves = get_sorted(stats, type='fastest_wave', recomputed=False) - axs[3].plot([me[0] for me in nwaves], [me[1] for me in nwaves], ls='-', color=color) - slowest_wave = get_sorted(stats, type='slowest_wave', recomputed=False) - # print(slowest_wave) - axs[3].plot([me[0] for me in slowest_wave], [me[1] for me in slowest_wave], ls='--', color=color) - axs[3].set_ylabel('waves') - # axs[3].set_yscale('log') - - axs[0].set_yscale('log') - axs[2].set_yscale('log') - axs[-1].set_xlabel('t') - axs[0].legend(frameon=False) - # print([(eps, data[eps]['t']) for eps in epsilons]) - plt.show() + fig.tight_layout() + save_fig(fig, 'parallel_efficiency', 'nprocs', 'speedup') if __name__ == "__main__": comm_world = MPI.COMM_WORLD - # vdp_stiffness_plot(runs=1, record=False) - # ERK_stiff_weirdness() - params = { - 'mode': 'RK_comp', - 'runs': 5, - #'num_procs': 1, # min(comm_world.size, 5), - 'plotting': comm_world.rank == 0, - } - params_single = { - **params, - 'problem': run_quench, - } record = True - single_problem(**params_single, work_key='t', precision_key='e_global', record=record) - # single_problem(**params_single, work_key='t', precision_key='e_global', record=False) + for mode in [ + # 'compare_strategies', + #'RK_comp_high_order', + # 'step_size_limiting', + # 'parallel_efficiency', + 'diagonal_SDC', + ]: + params = { + 'mode': mode, + 'runs': 1, + #'num_procs': 1, # min(comm_world.size, 5), + 'plotting': comm_world.rank == 0, + } + params_single = { + **params, + 'problem': run_Schroedinger, + } + single_problem(**params_single, work_key='t', precision_key='e_global_rel', record=record) all_params = { - 'record': True, - 'runs': 5, - 'work_key': 't', + 'record': False, + 'runs': 1, + 'work_key': 'k_linear', 'precision_key': 'e_global_rel', 'plotting': comm_world.rank == 0, #'num_procs': 4, } - for mode in ['RK_comp']: # ['parallel_efficiency', 'compare_strategies']: - break + for mode in [ + # 'RK_comp', + # 'parallel_efficiency', + # 'compare_adaptivity', + # 'compare_strategies', + # 'preconditioners', + # 'diagonal_SDC', + ]: all_problems(**all_params, mode=mode) comm_world.Barrier() diff --git a/pySDC/tests/test_projects/test_resilience/test_strategies.py b/pySDC/tests/test_projects/test_resilience/test_strategies.py index a8d7cce14a..1812f217cf 100644 --- a/pySDC/tests/test_projects/test_resilience/test_strategies.py +++ b/pySDC/tests/test_projects/test_resilience/test_strategies.py @@ -14,15 +14,17 @@ 'DIRK', 'explicitRK', 'ESDIRK', + 'AdaptivityPolynomialError', + 'kAdaptivity', ] STRATEGY_NAMES_NONMPIONLY = ['adaptiveHR', 'HotRod'] +STRATEGY_NAMES_MPIONLY = ['ARK'] LOGGER_LEVEL = 30 def single_test_vdp(strategy_name, useMPI, num_procs): import numpy as np from pySDC.helpers.stats_helper import get_sorted - from pySDC.projects.Resilience.vdp import run_vdp import pySDC.projects.Resilience.strategies as strategies from pySDC.implementations.hooks.log_work import LogWork from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun @@ -36,32 +38,45 @@ def single_test_vdp(strategy_name, useMPI, num_procs): # load the strategy avail_strategies = { - 'adaptivity': strategies.AdaptivityStrategy(useMPI=useMPI), - 'DIRK': strategies.DIRKStrategy(useMPI=useMPI), - 'adaptiveHR': strategies.AdaptiveHotRodStrategy(useMPI=useMPI), - 'iterate': strategies.IterateStrategy(useMPI=useMPI), - 'HotRod': strategies.HotRodStrategy(useMPI=useMPI), - 'explicitRK': strategies.ERKStrategy(useMPI=useMPI), - 'doubleAdaptivity': strategies.DoubleAdaptivityStrategy(useMPI=useMPI), - 'collocationRefinement': strategies.AdaptivityCollocationRefinementStrategy(useMPI=useMPI), - 'collocationDerefinement': strategies.AdaptivityCollocationDerefinementStrategy(useMPI=useMPI), - 'collocationType': strategies.AdaptivityCollocationTypeStrategy(useMPI=useMPI), - 'adaptivityAvoidRestarts': strategies.AdaptivityAvoidRestartsStrategy(useMPI=useMPI), - 'adaptivityInterpolation': strategies.AdaptivityInterpolationStrategy(useMPI=useMPI), - 'adaptivityQExtrapolation': strategies.AdaptivityExtrapolationWithinQStrategy(useMPI=useMPI), - 'base': strategies.BaseStrategy(useMPI=useMPI), - 'ESDIRK': strategies.ESDIRKStrategy(useMPI=useMPI), + 'adaptivity': strategies.AdaptivityStrategy, + 'DIRK': strategies.DIRKStrategy, + 'adaptiveHR': strategies.AdaptiveHotRodStrategy, + 'iterate': strategies.IterateStrategy, + 'HotRod': strategies.HotRodStrategy, + 'explicitRK': strategies.ERKStrategy, + 'doubleAdaptivity': strategies.DoubleAdaptivityStrategy, + 'collocationRefinement': strategies.AdaptivityCollocationRefinementStrategy, + 'collocationDerefinement': strategies.AdaptivityCollocationDerefinementStrategy, + 'collocationType': strategies.AdaptivityCollocationTypeStrategy, + 'adaptivityAvoidRestarts': strategies.AdaptivityAvoidRestartsStrategy, + 'adaptivityInterpolation': strategies.AdaptivityInterpolationStrategy, + 'adaptivityQExtrapolation': strategies.AdaptivityExtrapolationWithinQStrategy, + 'base': strategies.BaseStrategy, + 'ESDIRK': strategies.ESDIRKStrategy, + 'ARK': strategies.ARKStrategy, + 'AdaptivityPolynomialError': strategies.AdaptivityPolynomialError, + 'kAdaptivity': strategies.kAdaptivityStrategy, } if strategy_name in avail_strategies.keys(): - strategy = avail_strategies[strategy_name] + strategy = avail_strategies[strategy_name](useMPI=useMPI) else: raise NotImplementedError(f'Strategy \"{strategy_name}\" not implemented for this test!') - prob = run_vdp + if strategy_name in ['ARK']: + from pySDC.projects.Resilience.Schroedinger import run_Schroedinger as prob + else: + from pySDC.projects.Resilience.vdp import run_vdp as prob controller_params = {'logger_level': LOGGER_LEVEL} + custom_description = strategy.get_custom_description(problem=prob, num_procs=num_procs) + + if strategy_name == 'AdaptivityPolynomialError': + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError + + custom_description['convergence_controllers'][AdaptivityPolynomialError]['e_tol'] = 1e-4 + stats, _, Tend = prob( - custom_description=strategy.get_custom_description(problem=prob, num_procs=num_procs), + custom_description=custom_description, hook_class=[LogGlobalErrorPostRun, LogWork], use_MPI=useMPI, custom_controller_params=controller_params, @@ -80,14 +95,14 @@ def single_test_vdp(strategy_name, useMPI, num_procs): assert np.isclose( ref, act, rtol=1e-2 - ), f'Error in \"{strategy.name}\" strategy ({strategy_name})! Expected {key}={ref} but got {act}!' + ), f'Error in \"{strategy.name}\" strategy ({strategy_name})! Got {key}={act} but expected {ref}!' if comm: comm.Free() @pytest.mark.mpi4py -@pytest.mark.parametrize('strategy_name', STRATEGY_NAMES) +@pytest.mark.parametrize('strategy_name', STRATEGY_NAMES + STRATEGY_NAMES_MPIONLY) def test_strategy_with_vdp_MPI(strategy_name, num_procs=1): single_test_vdp(strategy_name=strategy_name, useMPI=True, num_procs=num_procs)