diff --git a/pySDC/implementations/problem_classes/LeakySuperconductor.py b/pySDC/implementations/problem_classes/Quench.py similarity index 69% rename from pySDC/implementations/problem_classes/LeakySuperconductor.py rename to pySDC/implementations/problem_classes/Quench.py index b3442d8afa..8f846127c4 100644 --- a/pySDC/implementations/problem_classes/LeakySuperconductor.py +++ b/pySDC/implementations/problem_classes/Quench.py @@ -9,7 +9,7 @@ # noinspection PyUnusedLocal -class LeakySuperconductor(ptype): +class Quench(ptype): """ This is a toy problem to emulate a magnet that has been cooled to temperatures where superconductivity is possible. However, there is a leak! Some point in the domain is constantly heated and when this has heated up its environment @@ -34,6 +34,7 @@ def __init__( Q_max=1.0, leak_range=(0.45, 0.55), leak_type='linear', + leak_transition='step', order=2, stencil_type='center', bc='neumann-zero', @@ -43,6 +44,7 @@ def __init__( lintol=1e-8, liniter=99, direct_solver=True, + reference_sol_type='scipy', ): """ Initialization routine @@ -62,6 +64,7 @@ def __init__( 'Q_max', 'leak_range', 'leak_type', + 'leak_transition', 'order', 'stencil_type', 'bc', @@ -71,6 +74,7 @@ def __init__( 'lintol', 'liniter', 'direct_solver', + 'reference_sol_type', localVars=locals(), readOnly=True, ) @@ -130,15 +134,17 @@ def eval_f_non_linear(self, u, t): elif self.leak_type == 'exponential': me[:] = Q_max * (np.exp(u) - np.exp(u_thresh)) / (np.exp(u_max) - np.exp(u_thresh)) else: - raise NotImplementedError(f'Leak type {self.leak_type} not implemented!') + raise NotImplementedError(f'Leak type \"{self.leak_type}\" not implemented!') me[u < u_thresh] = 0 - me[self.leak] = Q_max - me[u >= u_max] = Q_max + if self.leak_transition == 'step': + me[self.leak] = Q_max + elif self.leak_transition == 'Gaussian': + me[:] = np.max([me, Q_max * np.exp(-((self.xv - 0.5) ** 2) / 3e-2)], axis=0) + else: + raise NotImplementedError(f'Leak transition \"{self.leak_transition}\" not implemented!') - # boundary conditions - me[0] = 0.0 - me[-1] = 0.0 + me[u >= u_max] = Q_max me[:] /= self.Cv @@ -183,12 +189,14 @@ def get_non_linear_Jacobian(self, u): raise NotImplementedError(f'Leak type {self.leak_type} not implemented!') me[u < u_thresh] = 0 + if self.leak_transition == 'step': + me[self.leak] = 0 + elif self.leak_transition == 'Gaussian': + me[self.leak] = 0 + me[self.leak][u[self.leak] > Q_max * np.exp(-((self.xv[self.leak] - 0.5) ** 2) / 3e-2)] = 1 + else: + raise NotImplementedError(f'Leak transition \"{self.leak_transition}\" not implemented!') me[u > u_max] = 0 - me[self.leak] = 0 - - # boundary conditions - me[0] = 0.0 - me[-1] = 0.0 me[:] /= self.Cv @@ -270,38 +278,92 @@ def u_exact(self, t, u_init=None, t_init=None): me = self.dtype_u(self.init, val=0.0) if t > 0: + if self.reference_sol_type == 'scipy': + + def jac(t, u): + """ + Get the Jacobian for the implicit BDF method to use in `scipy.solve_ivp` + + Args: + t (float): The current time + u (dtype_u): Current solution + + Returns: + scipy.sparse.csc: The derivative of the non-linear part of the solution w.r.t. to the solution. + """ + return self.A + self.get_non_linear_Jacobian(u) + + def eval_rhs(t, u): + """ + Function to pass to `scipy.solve_ivp` to evaluate the full RHS + + Args: + t (float): Current time + u (numpy.1darray): Current solution + + Returns: + (numpy.1darray): RHS + """ + return self.eval_f(u.reshape(self.init[0]), t).flatten() + + me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init, method='BDF', jac=jac) + + elif self.reference_sol_type in ['DIRK', 'SDC']: + from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + from pySDC.implementations.hooks.log_solution import LogSolution + from pySDC.helpers.stats_helper import get_sorted + + description = {} + description['problem_class'] = Quench + description['problem_params'] = { + 'newton_tol': 1e-10, + 'newton_iter': 99, + 'nvars': 2**10, + **self.params, + } + + if self.reference_sol_type == 'DIRK': + from pySDC.implementations.sweeper_classes.Runge_Kutta import DIRK34 + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK + + description['sweeper_class'] = DIRK34 + description['sweeper_params'] = {} + description['step_params'] = {'maxiter': 1} + description['level_params'] = {'dt': 1e-4} + description['convergence_controllers'] = {AdaptivityRK: {'e_tol': 1e-9, 'update_order': 4}} + elif self.reference_sol_type == 'SDC': + from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + + description['sweeper_class'] = generic_implicit + description['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'quad_type': 'RADAU-RIGHT'} + description['step_params'] = {'maxiter': 99} + description['level_params'] = {'dt': 0.5, 'restol': 1e-10} + + controller_params = {'hook_class': LogSolution, 'mssdc_jac': False, 'logger_level': 99} + + controller = controller_nonMPI( + description=description, controller_params=controller_params, num_procs=1 + ) - def jac(t, u): - """ - Get the Jacobian for the implicit BDF method to use in `scipy.odeint` - - Args: - t (float): The current time - u (dtype_u): Current solution - - Returns: - scipy.sparse.csc: The derivative of the non-linear part of the solution w.r.t. to the solution. - """ - return self.A + self.get_non_linear_Jacobian(u) + uend, stats = controller.run( + u0=u_init if u_init is not None else self.u_exact(t=0.0), + t0=t_init if t_init is not None else 0, + Tend=t, + ) - def eval_rhs(t, u): - """ - Function to pass to `scipy.solve_ivp` to evaluate the full RHS + u_last = get_sorted(stats, type='u', recomputed=False)[-1] - Args: - t (float): Current time - u (numpy.1darray): Current solution + if abs(u_last[0] - t) > 1e-2: + self.logger.warning( + f'Time difference between reference solution and requested time is {abs(u_last[0]-t):.2e}!' + ) - Returns: - (numpy.1darray): RHS - """ - return self.eval_f(u.reshape(self.init[0]), t).flatten() + me[:] = u_last[1] - me[:] = self.generate_scipy_reference_solution(eval_rhs, t, u_init, t_init, method='BDF', jac=jac) return me -class LeakySuperconductorIMEX(LeakySuperconductor): +class QuenchIMEX(Quench): dtype_f = imex_mesh def eval_f(self, u, t): @@ -360,8 +422,7 @@ def u_exact(self, t, u_init=None, t_init=None): def jac(t, u): """ - Get the Jacobian for the implicit BDF method to use in `scipy.odeint` - + Get the Jacobian for the implicit BDF method to use in `scipy.solve_ivp` Args: t (float): The current time u (dtype_u): Current solution diff --git a/pySDC/projects/Resilience/Lorenz.py b/pySDC/projects/Resilience/Lorenz.py index a10b998655..83c5e1e756 100644 --- a/pySDC/projects/Resilience/Lorenz.py +++ b/pySDC/projects/Resilience/Lorenz.py @@ -9,6 +9,7 @@ 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 def run_Lorenz( @@ -18,7 +19,6 @@ def run_Lorenz( hook_class=LogData, fault_stuff=None, custom_controller_params=None, - custom_problem_params=None, use_MPI=False, **kwargs, ): @@ -32,7 +32,6 @@ def run_Lorenz( hook_class (pySDC.Hook): A hook to store data fault_stuff (dict): A dictionary with information on how to add faults custom_controller_params (dict): Overwrite presets - custom_problem_params (dict): Overwrite presets use_MPI (bool): Whether or not to use MPI Returns: @@ -56,9 +55,6 @@ def run_Lorenz( 'newton_maxiter': 99, } - if custom_problem_params is not None: - problem_params = {**problem_params, **custom_problem_params} - # initialize step parameters step_params = dict() step_params['maxiter'] = 4 @@ -82,11 +78,7 @@ def run_Lorenz( description['step_params'] = step_params if custom_description is not None: - for k in custom_description.keys(): - if k == 'sweeper_class': - description[k] = custom_description[k] - continue - description[k] = {**description.get(k, {}), **custom_description.get(k, {})} + description = merge_descriptions(description, custom_description) # set time parameters t0 = 0.0 @@ -98,18 +90,13 @@ def run_Lorenz( comm = kwargs.get('comm', MPI.COMM_WORLD) controller = controller_MPI(controller_params=controller_params, description=description, comm=comm) - - # get initial values on finest level P = controller.S.levels[0].prob - uinit = P.u_exact(t0) else: controller = controller_nonMPI( num_procs=num_procs, controller_params=controller_params, description=description ) - - # get initial values on finest level P = controller.MS[0].levels[0].prob - uinit = P.u_exact(t0) + uinit = P.u_exact(t0) # insert faults if fault_stuff is not None: diff --git a/pySDC/projects/Resilience/Schroedinger.py b/pySDC/projects/Resilience/Schroedinger.py index 4a90674da4..69bf81f8b0 100644 --- a/pySDC/projects/Resilience/Schroedinger.py +++ b/pySDC/projects/Resilience/Schroedinger.py @@ -1,6 +1,5 @@ -import numpy as np -from pathlib import Path from mpi4py import MPI +import numpy as np from pySDC.helpers.stats_helper import get_sorted @@ -9,6 +8,58 @@ from pySDC.implementations.problem_classes.NonlinearSchroedinger_MPIFFT import nonlinearschroedinger_imex 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.core.Hooks import hooks + +import matplotlib.pyplot as plt +from mpl_toolkits.axes_grid1 import make_axes_locatable + + +class live_plotting_with_error(hooks): # pragma: no cover + def __init__(self): + super().__init__() + self.fig, self.axs = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(12, 7)) + + divider = make_axes_locatable(self.axs[1]) + self.cax_right = divider.append_axes('right', size='5%', pad=0.05) + divider = make_axes_locatable(self.axs[0]) + self.cax_left = divider.append_axes('right', size='5%', pad=0.05) + + def post_step(self, step, level_number): + lvl = step.levels[level_number] + lvl.sweep.compute_end_point() + + self.axs[0].cla() + im1 = self.axs[0].imshow(np.abs(lvl.uend), vmin=0, vmax=2.0) + self.fig.colorbar(im1, cax=self.cax_left) + + self.axs[1].cla() + im = self.axs[1].imshow(np.abs(lvl.prob.u_exact(lvl.time + lvl.dt) - lvl.uend)) + self.fig.colorbar(im, cax=self.cax_right) + + self.fig.suptitle(f't={lvl.time:.2f}') + self.axs[0].set_title('solution') + self.axs[1].set_title('error') + plt.pause(1e-9) + + +class live_plotting(hooks): # pragma: no cover + def __init__(self): + super().__init__() + self.fig, self.ax = plt.subplots() + divider = make_axes_locatable(self.ax) + self.cax = divider.append_axes('right', size='5%', pad=0.05) + + def post_step(self, step, level_number): + lvl = step.levels[level_number] + lvl.sweep.compute_end_point() + + self.ax.cla() + im = self.ax.imshow(np.abs(lvl.uend), vmin=0.2, vmax=1.8) + self.ax.set_title(f't={lvl.time + lvl.dt:.2f}') + self.fig.colorbar(im, cax=self.cax) + plt.pause(1e-9) def run_Schroedinger( @@ -18,7 +69,6 @@ def run_Schroedinger( hook_class=LogData, fault_stuff=None, custom_controller_params=None, - custom_problem_params=None, use_MPI=False, space_comm=None, **kwargs, @@ -33,7 +83,6 @@ def run_Schroedinger( hook_class (pySDC.Hook): A hook to store data fault_stuff (dict): A dictionary with information on how to add faults custom_controller_params (dict): Overwrite presets - custom_problem_params (dict): Overwrite presets use_MPI (bool): Whether or not to use MPI Returns: @@ -41,13 +90,14 @@ def run_Schroedinger( controller: The controller Tend: The time that was supposed to be integrated to """ + from mpi4py import MPI - space_comm = MPI.COMM_WORLD if space_comm is None else space_comm + space_comm = MPI.COMM_SELF if space_comm is None else space_comm rank = space_comm.Get_rank() # initialize level parameters level_params = dict() - level_params['restol'] = 1e-08 + level_params['restol'] = 1e-8 level_params['dt'] = 1e-01 / 2 level_params['nsweeps'] = 1 @@ -62,11 +112,9 @@ def run_Schroedinger( problem_params = dict() problem_params['nvars'] = (128, 128) problem_params['spectral'] = False + problem_params['c'] = 1.0 problem_params['comm'] = space_comm - if custom_problem_params is not None: - problem_params = {**problem_params, **custom_problem_params} - # initialize step parameters step_params = dict() step_params['maxiter'] = 50 @@ -90,21 +138,26 @@ def run_Schroedinger( description['step_params'] = step_params if custom_description is not None: - for k in custom_description.keys(): - if type(custom_description[k]) == dict: - description[k] = {**description.get(k, {}), **custom_description.get(k, {})} - else: - description[k] = custom_description[k] + description = merge_descriptions(description, custom_description) # set time parameters t0 = 0.0 # instantiate controller - assert use_MPI == False, "MPI version in time not implemented" - controller = controller_nonMPI(num_procs=num_procs, controller_params=controller_params, description=description) + controller_args = { + 'controller_params': controller_params, + 'description': description, + } + if use_MPI: + from pySDC.implementations.controller_classes.controller_MPI import controller_MPI + + comm = kwargs.get('comm', MPI.COMM_WORLD) + controller = controller_MPI(**controller_args, comm=comm) + P = controller.S.levels[0].prob + else: + controller = controller_nonMPI(**controller_args, num_procs=num_procs) + P = controller.MS[0].levels[0].prob - # get initial values on finest level - P = controller.MS[0].levels[0].prob uinit = P.u_exact(t0) # insert faults @@ -124,25 +177,9 @@ def run_Schroedinger( return stats, controller, Tend -def plot_solution(stats): # pragma: no cover - import matplotlib.pyplot as plt - - u = get_sorted(stats, type='u') - fig, axs = plt.subplots(1, 2, figsize=(12, 5)) - axs[0].imshow(np.abs(u[0][1])) - axs[0].set_title(f't={u[0][0]:.2f}') - for i in range(len(u)): - axs[1].cla() - axs[1].imshow(np.abs(u[i][1])) - axs[1].set_title(f't={u[i][0]:.2f}') - plt.pause(1e-1) - fig.tight_layout() - plt.show() - - def main(): - stats, _, _ = run_Schroedinger(space_comm=MPI.COMM_WORLD) - plot_solution(stats) + stats, _, _ = run_Schroedinger(space_comm=MPI.COMM_WORLD, hook_class=live_plotting) + plt.show() if __name__ == "__main__": diff --git a/pySDC/projects/Resilience/accuracy_check.py b/pySDC/projects/Resilience/accuracy_check.py index fdfbc2af75..d79d8c0085 100644 --- a/pySDC/projects/Resilience/accuracy_check.py +++ b/pySDC/projects/Resilience/accuracy_check.py @@ -118,12 +118,14 @@ def multiple_runs( num_procs = 1 if serial else 5 + embedded_error_flavor = 'standard' if avoid_restarts else 'linearized' + # perform rest of the tests for i in range(0, len(dt_list)): desc = { 'step_params': {'maxiter': k}, 'convergence_controllers': { - EstimateEmbeddedError: {}, + EstimateEmbeddedError.get_implementation(flavor=embedded_error_flavor, useMPI=False): {}, EstimateExtrapolationErrorNonMPI: {'no_storage': not serial}, }, } @@ -134,7 +136,11 @@ def multiple_runs( elif var == 'e_tol': from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity - desc['convergence_controllers'][Adaptivity] = {'e_tol': dt_list[i], 'avoid_restarts': avoid_restarts} + desc['convergence_controllers'][Adaptivity] = { + 'e_tol': dt_list[i], + 'avoid_restarts': avoid_restarts, + 'embedded_error_flavor': embedded_error_flavor, + } if custom_description is not None: desc = {**desc, **custom_description} @@ -379,7 +385,16 @@ def check_order_with_adaptivity(): ks = [3, 2] for serial in [True, False]: fig, ax = plt.subplots(1, 1, figsize=(3.5, 3)) - plot_all_errors(ax, ks, serial, Tend_fixed=5e-1, var='e_tol', dt_list=[1e-5, 1e-6, 1e-7], avoid_restarts=True) + plot_all_errors( + ax, + ks, + serial, + Tend_fixed=5e-1, + var='e_tol', + dt_list=[1e-5, 5e-6], + avoid_restarts=False, + custom_controller_params={'logger_level': 30}, + ) if serial: fig.savefig('data/error_estimate_order_adaptivity.png', dpi=300, bbox_inches='tight') else: diff --git a/pySDC/projects/Resilience/advection.py b/pySDC/projects/Resilience/advection.py index b267dc4e72..300c81b8d2 100644 --- a/pySDC/projects/Resilience/advection.py +++ b/pySDC/projects/Resilience/advection.py @@ -7,6 +7,7 @@ import numpy as np from pySDC.projects.Resilience.hook import LogData, hook_collection from pySDC.projects.Resilience.fault_injection import prepare_controller_for_faults +from pySDC.projects.Resilience.strategies import merge_descriptions def plot_embedded(stats, ax): @@ -28,7 +29,6 @@ def run_advection( hook_class=LogData, fault_stuff=None, custom_controller_params=None, - custom_problem_params=None, ): # initialize level parameters level_params = dict() @@ -40,10 +40,7 @@ def run_advection( sweeper_params['num_nodes'] = 3 sweeper_params['QI'] = 'IE' - problem_params = {'freq': 2, 'nvars': 2**9, 'c': 1.0, 'order': 5, 'bc': 'periodic'} - - if custom_problem_params is not None: - problem_params = {**problem_params, **custom_problem_params} + problem_params = {'freq': 2, 'nvars': 2**9, 'c': 1.0, 'stencil_type': 'center', 'order': 4, 'bc': 'periodic'} # initialize step parameters step_params = dict() @@ -60,19 +57,15 @@ def run_advection( # fill description dictionary for easy step instantiation description = dict() - description['problem_class'] = advectionNd # pass problem class - description['problem_params'] = problem_params # pass problem parameters - description['sweeper_class'] = generic_implicit # pass sweeper - description['sweeper_params'] = sweeper_params # pass sweeper parameters - description['level_params'] = level_params # pass level parameters + description['problem_class'] = advectionNd + description['problem_params'] = problem_params + description['sweeper_class'] = generic_implicit + description['sweeper_params'] = sweeper_params + description['level_params'] = level_params description['step_params'] = step_params if custom_description is not None: - for k in custom_description.keys(): - if k == 'sweeper_class': - description[k] = custom_description[k] - continue - description[k] = {**description.get(k, {}), **custom_description.get(k, {})} + description = merge_descriptions(description, custom_description) # set time parameters t0 = 0.0 diff --git a/pySDC/projects/Resilience/dahlquist.py b/pySDC/projects/Resilience/dahlquist.py index 8837d5ce87..5b1ebc2189 100644 --- a/pySDC/projects/Resilience/dahlquist.py +++ b/pySDC/projects/Resilience/dahlquist.py @@ -11,6 +11,7 @@ from pySDC.implementations.hooks.log_solution import LogSolutionAfterIteration from pySDC.implementations.hooks.log_step_size import LogStepSize +from pySDC.projects.Resilience.strategies import merge_descriptions class LogLambdas(hooks): @@ -34,7 +35,6 @@ def run_dahlquist( hook_class=hooks, fault_stuff=None, custom_controller_params=None, - custom_problem_params=None, **kwargs, ): """ @@ -47,7 +47,6 @@ def run_dahlquist( hook_class (pySDC.Hook): A hook to store data fault_stuff (dict): A dictionary with information on how to add faults custom_controller_params (dict): Overwrite presets - custom_problem_params (dict): Overwrite presets Returns: dict: The stats object @@ -78,9 +77,6 @@ def run_dahlquist( 'u0': 1.0 + 0.0j, } - if custom_problem_params is not None: - problem_params = {**problem_params, **custom_problem_params} - # initialize step parameters step_params = dict() step_params['maxiter'] = 5 @@ -104,11 +100,7 @@ def run_dahlquist( description['step_params'] = step_params if custom_description is not None: - for k in custom_description.keys(): - if k == 'sweeper_class': - description[k] = custom_description[k] - continue - description[k] = {**description.get(k, {}), **custom_description.get(k, {})} + description = merge_descriptions(description, custom_description) # set time parameters t0 = 0.0 diff --git a/pySDC/projects/Resilience/extrapolation_within_Q.py b/pySDC/projects/Resilience/extrapolation_within_Q.py index de92b5541d..506d31e69a 100644 --- a/pySDC/projects/Resilience/extrapolation_within_Q.py +++ b/pySDC/projects/Resilience/extrapolation_within_Q.py @@ -32,6 +32,9 @@ def multiple_runs(prob, dts, num_nodes, quad_type='RADAU-RIGHT'): description['sweeper_params'] = {'num_nodes': num_nodes, 'quad_type': quad_type} description['convergence_controllers'] = {EstimateExtrapolationErrorWithinQ: {}} + if prob.__name__ == 'run_advection': + description['problem_params'] = {'order': 6, 'stencil_type': 'center'} + res = {} for dt in dts: @@ -104,9 +107,9 @@ def check_order(ax, prob, dts, num_nodes, quad_type): def main(): fig, ax = plt.subplots() - num_nodes = 2 + num_nodes = 3 quad_type = 'RADAU-RIGHT' - check_order(ax, run_vdp, [5e-1, 1e-1, 5e-2, 1e-2], num_nodes, quad_type) + check_order(ax, run_advection, [5e-1, 1e-1, 5e-2, 1e-2], num_nodes, quad_type) plt.show() diff --git a/pySDC/projects/Resilience/fault_injection.py b/pySDC/projects/Resilience/fault_injection.py index f70adacbe3..84e513827d 100644 --- a/pySDC/projects/Resilience/fault_injection.py +++ b/pySDC/projects/Resilience/fault_injection.py @@ -254,6 +254,8 @@ def inject_fault(self, step, f): None ''' L = step.levels[f.level_number] + _abs_before = None + _abs_after = None # insert the fault in some target if f.target == 0: @@ -268,14 +270,18 @@ def inject_fault(self, step, f): fault happens in the last iteration, it will not show up in the residual and the iteration is wrongly stopped. ''' + _abs_before = abs(L.u[f.node][tuple(f.problem_pos)]) L.u[f.node][tuple(f.problem_pos)] = self.flip_bit(L.u[f.node][tuple(f.problem_pos)], f.bit) L.f[f.node] = L.prob.eval_f(L.u[f.node], L.time + L.dt * L.sweep.coll.nodes[max([0, f.node - 1])]) L.sweep.compute_residual() + _abs_after = abs(L.u[f.node][tuple(f.problem_pos)]) else: raise NotImplementedError(f'Target {f.target} for faults not implemented!') # log what happened to stats and screen - self.logger.info(f'Flipping bit {f.bit} {f.when} iteration {f.iteration} in node {f.node}. Target: {f.target}') + self.logger.info( + f'Flipping bit {f.bit} {f.when} iteration {f.iteration} in node {f.node}. Target: {f.target}. Abs: {_abs_before:.2e} -> {_abs_after:.2e}' + ) self.add_to_stats( process=step.status.slot, time=L.time, @@ -326,6 +332,7 @@ def pre_run(self, step, level_number): 'iteration': step.params.maxiter, 'problem_pos': step.levels[level_number].u[0].shape, 'bit': bit, # change manually if you ever have something else + **self.rnd_params, } # initialize the faults have been added before we knew the random parameters @@ -420,7 +427,8 @@ def post_iteration(self, step, level_number): return None - def to_binary(self, f): + @classmethod + def to_binary(cls, f): ''' Converts a single float in a string containing its binary representation in memory following IEEE754 The struct.pack function returns the input with the applied conversion code in 8 bit blocks, which are then @@ -437,13 +445,14 @@ def to_binary(self, f): elif type(f) in [np.float32]: conversion_code = '>f' # big endian, float elif type(f) in [np.complex128]: - return f'{self.to_binary(f.real)}{self.to_binary(f.imag)}' + return f'{cls.to_binary(f.real)}{cls.to_binary(f.imag)}' else: raise NotImplementedError(f'Don\'t know how to convert number of type {type(f)} to binary') return ''.join('{:0>8b}'.format(c) for c in struct.pack(conversion_code, f)) - def to_float(self, s): + @classmethod + def to_float(cls, s): ''' Converts a string of a IEEE754 binary representation in a float. The string is converted to integer with base 2 and converted to bytes, which can be unpacked into a Python float by the struct module. @@ -463,14 +472,15 @@ def to_float(self, s): elif len(s) == 128: # complex floats real = s[0:64] imag = s[64:128] - return self.to_float(real) + self.to_float(imag) * 1j + return cls.to_float(real) + cls.to_float(imag) * 1j else: raise NotImplementedError(f'Don\'t know how to convert string of length {len(s)} to float') return struct.unpack(conversion_code, int(s, 2).to_bytes(byte_count, 'big'))[0] - def flip_bit(self, target, bit): + @classmethod + def flip_bit(cls, target, bit): ''' Flips a bit at position bit in a target using the bitwise xor operator @@ -481,8 +491,8 @@ def flip_bit(self, target, bit): Returns: (float) The floating point number resulting from flipping the respective bit in target ''' - binary = self.to_binary(target) - return self.to_float(f'{binary[:bit]}{int(binary[bit]) ^ 1}{binary[bit+1:]}') + binary = cls.to_binary(target) + return cls.to_float(f'{binary[:bit]}{int(binary[bit]) ^ 1}{binary[bit+1:]}') def prepare_controller_for_faults(controller, fault_stuff, rnd_args, args): @@ -501,10 +511,19 @@ def prepare_controller_for_faults(controller, fault_stuff, rnd_args, args): """ faultHook = get_fault_injector_hook(controller) faultHook.random_generator = fault_stuff['rng'] - faultHook.add_fault( - rnd_args={**rnd_args, **fault_stuff.get('rnd_params', {})}, - args={**args, **fault_stuff.get('args', {})}, - ) + + for key in ['fault_frequency_iter']: + if key in fault_stuff.keys(): + faultHook.__dict__[key] = fault_stuff[key] + + for key, val in fault_stuff.get('rnd_params', {}).items(): + faultHook.rnd_params[key] = val + + if not len(faultHook.rnd_params.keys()) > 0: + faultHook.add_fault( + rnd_args={**rnd_args, **fault_stuff.get('rnd_params', {})}, + args={**args, **fault_stuff.get('args', {})}, + ) def get_fault_injector_hook(controller): diff --git a/pySDC/projects/Resilience/fault_stats.py b/pySDC/projects/Resilience/fault_stats.py index ab1c5d933a..917d616191 100644 --- a/pySDC/projects/Resilience/fault_stats.py +++ b/pySDC/projects/Resilience/fault_stats.py @@ -1,10 +1,7 @@ import numpy as np import pickle import matplotlib.pyplot as plt -from matplotlib.colors import TABLEAU_COLORS from mpi4py import MPI -import sys -import matplotlib as mpl import pySDC.helpers.plot_helper as plot_helper from pySDC.helpers.stats_helper import get_sorted @@ -13,7 +10,6 @@ from pySDC.projects.Resilience.fault_injection import get_fault_injector_hook from pySDC.implementations.convergence_controller_classes.hotrod import HotRod from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity -from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI from pySDC.implementations.hooks.log_errors import LogLocalErrorPostStep from pySDC.implementations.hooks.log_work import LogWork @@ -23,381 +19,11 @@ from pySDC.projects.Resilience.piline import run_piline from pySDC.projects.Resilience.Lorenz import run_Lorenz from pySDC.projects.Resilience.Schroedinger import run_Schroedinger -from pySDC.projects.Resilience.leaky_superconductor import run_leaky_superconductor +from pySDC.projects.Resilience.quench import run_quench -plot_helper.setup_mpl(reset=True) -cmap = TABLEAU_COLORS - - -class Strategy: - ''' - Abstract class for resilience strategies - ''' - - def __init__(self): - ''' - Initialization routine - ''' - - # set default values for plotting - self.linestyle = '-' - self.marker = '.' - self.name = '' - self.bar_plot_x_label = '' - self.color = list(cmap.values())[0] - - # setup custom descriptions - self.custom_description = {} - - # prepare parameters for masks to identify faults that cannot be fixed by this strategy - self.fixable = [] - self.fixable += [ - { - 'key': 'node', - 'op': 'gt', - 'val': 0, - } - ] - self.fixable += [ - { - 'key': 'error', - 'op': 'isfinite', - } - ] - - def get_fixable_params(self, **kwargs): - """ - Return a list containing dictionaries which can be passed to `FaultStats.get_mask` as keyword arguments to - obtain a mask of faults that can be fixed - - Returns: - list: Dictionary of parameters - """ - return self.fixable - - def get_custom_description(self, problem, num_procs): - ''' - Routine to get a custom description that realizes the resilience strategy and tailors it to the problem at hand - - 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: The custom descriptions you can supply to the problem when running it - ''' - - return self.custom_description - - def get_fault_args(self, problem, num_procs): - ''' - Routine to get arguments for the faults that are exempt from randomization - - 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: Arguments for the faults that are exempt from randomization - ''' - - return {} - - 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 - ''' - - return {} - - @property - def style(self): - """ - Get the plotting parameters for the strategy. - Supply them to a plotting function using `**` - - Returns: - (dict): The plotting parameters as a dictionary - """ - return { - 'marker': self.marker, - 'label': self.label, - 'color': self.color, - 'ls': self.linestyle, - } - - @property - def label(self): - """ - Get a label for plotting - """ - return self.name - - -class BaseStrategy(Strategy): - ''' - Do a fixed iteration count - ''' - - def __init__(self): - ''' - Initialization routine - ''' - super(BaseStrategy, self).__init__() - self.color = list(cmap.values())[0] - self.marker = 'o' - self.name = 'base' - self.bar_plot_x_label = 'base' - - -class AdaptivityStrategy(Strategy): - ''' - Adaptivity as a resilience strategy - ''' - - def __init__(self): - ''' - Initialization routine - ''' - super(AdaptivityStrategy, self).__init__() - self.color = list(cmap.values())[1] - self.marker = '*' - self.name = 'adaptivity' - self.bar_plot_x_label = 'adaptivity' - - def get_fixable_params(self, maxiter, **kwargs): - """ - Here faults occurring in the last iteration cannot be fixed. - - Args: - maxiter (int): Max. iterations until convergence is declared - - Returns: - (list): Contains dictionaries of keyword arguments for `FaultStats.get_mask` - """ - self.fixable += [ - { - 'key': 'iteration', - 'op': 'lt', - 'val': maxiter, - } - ] - return self.fixable - - 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 - ''' - custom_description = {} - - dt_max = np.inf - dt_min = 1e-5 - - if problem == run_piline: - e_tol = 1e-7 - dt_min = 1e-2 - elif problem == run_vdp: - e_tol = 2e-5 - dt_min = 1e-3 - elif problem == run_Lorenz: - e_tol = 2e-5 - dt_min = 1e-3 - elif problem == run_Schroedinger: - e_tol = 4e-6 - dt_min = 1e-3 - elif problem == run_leaky_superconductor: - e_tol = 1e-7 - dt_min = 1e-3 - dt_max = 1e2 - else: - raise NotImplementedError( - 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\ - strategy' - ) - - custom_description['convergence_controllers'] = { - Adaptivity: {'e_tol': e_tol, 'dt_min': dt_min, 'dt_max': dt_max} - } - return {**custom_description, **self.custom_description} - - -class AdaptiveHotRodStrategy(Strategy): - ''' - Adaptivity + Hot Rod as a resilience strategy - ''' - - def __init__(self): - ''' - Initialization routine - ''' - super(AdaptiveHotRodStrategy, self).__init__() - self.color = list(cmap.values())[4] - self.marker = '.' - self.name = 'adaptive Hot Rod' - self.bar_plot_x_label = 'adaptive\nHot Rod' - - def get_custom_description(self, problem, num_procs): - ''' - Routine to get a custom description that adds adaptivity and Hot Rod - - 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 description you can supply to the problem when running it - ''' - if problem == run_vdp: - e_tol = 3e-7 - dt_min = 1e-3 - maxiter = 4 - HotRod_tol = 3e-7 - else: - raise NotImplementedError( - 'I don\'t have a tolerance for adaptive Hot Rod for your problem. Please add one \ -to the strategy' - ) - - no_storage = num_procs > 1 - - custom_description = { - 'convergence_controllers': { - Adaptivity: {'e_tol': e_tol, 'dt_min': dt_min}, - HotRod: {'HotRod_tol': HotRod_tol, 'no_storage': no_storage}, - }, - 'step_params': {'maxiter': maxiter}, - } - - return {**custom_description, **self.custom_description} +from pySDC.projects.Resilience.strategies import BaseStrategy, AdaptivityStrategy, IterateStrategy, HotRodStrategy - -class IterateStrategy(Strategy): - ''' - Iterate for as much as you want - ''' - - def __init__(self): - ''' - Initialization routine - ''' - super(IterateStrategy, self).__init__() - self.color = list(cmap.values())[2] - self.marker = 'v' - self.name = 'iterate' - self.bar_plot_x_label = 'iterate' - - @property - def label(self): - return r'$k$ adaptivity' - - def get_custom_description(self, problem, num_procs): - ''' - Routine to get a custom description that allows for adaptive iteration counts - - 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 description you can supply to the problem when running it - ''' - restol = -1 - e_tol = -1 - - if problem == run_piline: - restol = 2.3e-8 - elif problem == run_vdp: - restol = 9e-7 - elif problem == run_Lorenz: - restol = 16e-7 - elif problem == run_Schroedinger: - restol = 6.5e-7 - elif problem == run_leaky_superconductor: - # e_tol = 1e-6 - restol = 1e-11 - else: - raise NotImplementedError( - 'I don\'t have a residual tolerance for your problem. Please add one to the \ -strategy' - ) - - custom_description = { - 'step_params': {'maxiter': 99}, - 'level_params': {'restol': restol, 'e_tol': e_tol}, - } - - return {**custom_description, **self.custom_description} - - -class HotRodStrategy(Strategy): - ''' - Hot Rod as a resilience strategy - ''' - - def __init__(self): - ''' - Initialization routine - ''' - super(HotRodStrategy, self).__init__() - self.color = list(cmap.values())[3] - self.marker = '^' - self.name = 'Hot Rod' - self.bar_plot_x_label = 'Hot Rod' - - def get_custom_description(self, problem, num_procs): - ''' - Routine to get a custom description that adds Hot Rod - - 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 description you can supply to the problem when running it - ''' - if problem == run_vdp: - HotRod_tol = 5e-7 - maxiter = 4 - elif problem == run_Lorenz: - HotRod_tol = 4e-7 - maxiter = 6 - elif problem == run_Schroedinger: - HotRod_tol = 3e-7 - maxiter = 6 - elif problem == run_leaky_superconductor: - HotRod_tol = 3e-5 - maxiter = 6 - else: - raise NotImplementedError( - 'I don\'t have a tolerance for Hot Rod for your problem. Please add one to the\ - strategy' - ) - - no_storage = num_procs > 1 - - custom_description = { - 'convergence_controllers': { - HotRod: {'HotRod_tol': HotRod_tol, 'no_storage': no_storage}, - BasicRestartingNonMPI: {'max_restarts': 2, 'crash_after_max_restarts': False}, - }, - 'step_params': {'maxiter': maxiter}, - } - - return {**custom_description, **self.custom_description} +plot_helper.setup_mpl(reset=True) class FaultStats: @@ -412,10 +38,11 @@ def __init__( faults=None, reload=True, recovery_thresh=1 + 1e-3, - recovery_thresh_abs=1e9, + recovery_thresh_abs=0.0, num_procs=1, mode='combination', stats_path='data/stats', + **kwargs, ): ''' Initialization routine @@ -438,6 +65,10 @@ def __init__( self.num_procs = num_procs self.mode = mode self.stats_path = stats_path + self.kwargs = { + 'fault_frequency_iter': 500, + **kwargs, + } def get_Tend(self): ''' @@ -446,56 +77,9 @@ def get_Tend(self): Returns: float: Tend to put into the run ''' - if self.prob == run_vdp: - return 2.3752559741400825 - elif self.prob == run_piline: - return 20.0 - elif self.prob == run_Lorenz: - return 1.5 - elif self.prob == run_Schroedinger: - return 1.0 - elif self.prob == run_leaky_superconductor: - return 450 - else: - raise NotImplementedError('I don\'t have a final time for your problem!') - - def get_custom_description(self): - ''' - Get a custom description based on the problem - - Returns: - dict: Custom description - ''' - custom_description = {} - if self.prob == run_vdp: - custom_description['step_params'] = {'maxiter': 3} - elif self.prob == run_Lorenz: - custom_description['step_params'] = {'maxiter': 5} - elif self.prob == run_Schroedinger: - custom_description['step_params'] = {'maxiter': 5} - custom_description['level_params'] = {'dt': 1e-2, 'restol': -1} - elif self.prob == run_leaky_superconductor: - custom_description['level_params'] = {'restol': -1, 'dt': 10.0} - custom_description['step_params'] = {'maxiter': 5} - custom_description['problem_params'] = {'newton_iter': 99, 'newton_tol': 1e-10} - return custom_description - - def get_custom_problem_params(self): - ''' - Get a custom problem parameters based on the problem - - Returns: - dict: Custom problem params - ''' - custom_params = {} - if self.prob == run_vdp: - custom_params = { - 'u0': np.array([0.99995, -0.00999985], dtype=np.float64), - 'crash_at_maxiter': False, - } - return custom_params + return self.strategies[0].get_Tend(self.prob, self.num_procs) - def run_stats_generation(self, runs=1000, step=None, comm=None, _reload=False, _runs_partial=0): + def run_stats_generation(self, runs=1000, step=None, comm=None, kwargs_range=None, _reload=False, _runs_partial=0): ''' Run the generation of stats for all strategies in the `self.strategies` variable @@ -503,11 +87,21 @@ def run_stats_generation(self, runs=1000, step=None, comm=None, _reload=False, _ runs (int): Number of runs you want to do step (int): Number of runs you want to do between saving comm (MPI.Communicator): Communicator for distributing runs + kw_args_range (dict): Range for the parameters _reload, _runs_partial: Variables only used for recursion. Do not change! Returns: None ''' + for key, val in kwargs_range.items() if kwargs_range is not None else {}: + if type(val) == int: + self.kwargs[key] = val + else: + for me in val: + kwargs_range_me = {**kwargs_range, key: me} + self.run_stats_generation(runs=runs, step=step, comm=comm, kwargs_range=kwargs_range_me) + 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 _runs_partial = step if _runs_partial == 0 else _runs_partial @@ -519,9 +113,11 @@ def run_stats_generation(self, runs=1000, step=None, comm=None, _reload=False, _ # sort the strategies to do some load balancing sorting_index = None if comm.rank == 0: - already_completed = np.array([self.load(strategy, True).get('runs', 0) for strategy in self.strategies]) + already_completed = np.array( + [self.load(strategy=strategy, faults=True).get('runs', 0) for strategy in self.strategies] + ) sorting_index_ = np.argsort(already_completed) - sorting_index = sorting_index_[already_completed[sorting_index_] < runs] + sorting_index = sorting_index_[already_completed[sorting_index_] < max_runs] # tell all ranks what strategies to use sorting_index = comm.bcast(sorting_index, root=0) @@ -540,7 +136,7 @@ def run_stats_generation(self, runs=1000, step=None, comm=None, _reload=False, _ else: runs_partial = min([5, _runs_partial]) self.generate_stats( - strategy=strategies[j + comm.rank % len(strategies)], + strategy=strategies[j + (comm.rank % len(strategies) % (len(strategies)) - j)], runs=runs_partial, faults=f, reload=reload, @@ -583,11 +179,17 @@ def generate_stats(self, strategy=None, runs=1000, reload=True, faults=True, com 'target': np.zeros(runs), } + # store arguments for storing and loading + identifier_args = { + 'faults': faults, + 'strategy': strategy, + } + # reload previously recorded stats and write them to dat if reload: already_completed_ = None if comm.rank == 0: - already_completed_ = self.load(strategy, faults) + already_completed_ = self.load(**identifier_args) already_completed = comm.bcast(already_completed_, root=0) if already_completed['runs'] > 0 and already_completed['runs'] <= runs and comm.rank == 0: for k in dat.keys(): @@ -625,13 +227,15 @@ def generate_stats(self, strategy=None, runs=1000, reload=True, faults=True, com # record the new data point if faults: - assert len(faults_run) > 0, f'No faults where recorded in run {i} of strategy {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] + 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] + 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 @@ -646,10 +250,10 @@ def generate_stats(self, strategy=None, runs=1000, reload=True, faults=True, com if already_completed['runs'] < runs: if comm.rank == 0: - self.store(strategy, faults, dat_full) + self.store(dat_full, **identifier_args) if self.faults: try: - self.get_recovered(strategy) + self.get_recovered(strategy=strategy) except KeyError: print('Warning: Can\'t compute recovery rate right now') @@ -668,18 +272,9 @@ def get_error(self, u, t, controller, strategy): Returns: float: Error """ - if self.prob == run_leaky_superconductor: - ref = { - AdaptivityStrategy: 0.036832240840408426, - IterateStrategy: 0.0368214748207781, - HotRodStrategy: 0.03682153860683977, - BaseStrategy: 0.03682153860683977, - } - return abs(max(u) - ref[type(strategy)]) - else: - return abs(u - controller.MS[0].levels[0].prob.u_exact(t=t)) + return abs(u - controller.MS[0].levels[0].prob.u_exact(t=t)) - def single_run(self, strategy, run=0, faults=False, force_params=None, hook_class=None, space_comm=None): + def single_run(self, strategy, run=0, faults=False, force_params=None, hook_class=None, space_comm=None, Tend=None): ''' Run the problem once with the specified parameters @@ -699,34 +294,34 @@ def single_run(self, strategy, run=0, faults=False, force_params=None, hook_clas force_params = {} if force_params is None else force_params # build the custom description - custom_description_prob = self.get_custom_description() - custom_description_strategy = strategy.get_custom_description(self.prob, self.num_procs) - custom_description = {} - for key in list(custom_description_strategy.keys()) + list(custom_description_prob.keys()): - custom_description[key] = { - **custom_description_prob.get(key, {}), - **custom_description_strategy.get(key, {}), - } + custom_description = strategy.get_custom_description(self.prob, self.num_procs) for k in force_params.keys(): custom_description[k] = {**custom_description.get(k, {}), **force_params[k]} custom_controller_params = force_params.get('controller_params', {}) - custom_problem_params = self.get_custom_problem_params() if faults: + fault_stuff = { + 'rng': None, + 'args': strategy.get_fault_args(self.prob, self.num_procs), + 'rnd_params': strategy.get_fault_args(self.prob, self.num_procs), + } + # make parameters for faults: if self.mode == 'random': - rng = np.random.RandomState(run) + fault_stuff['rng'] = np.random.RandomState(run) elif self.mode == 'combination': - rng = run + fault_stuff['rng'] = run + elif self.mode == 'regular': + fault_stuff['rng'] = np.random.RandomState(run) + fault_stuff['fault_frequency_iter'] = self.kwargs['fault_frequency_iter'] + fault_stuff['rnd_params'] = { + 'bit': 12, + 'min_node': 1, + } else: raise NotImplementedError(f'Don\'t know how to add faults in mode {self.mode}') - fault_stuff = { - 'rng': rng, - 'args': strategy.get_fault_args(self.prob, self.num_procs), - 'rnd_params': strategy.get_fault_args(self.prob, self.num_procs), - } else: fault_stuff = None @@ -735,13 +330,12 @@ def single_run(self, strategy, run=0, faults=False, force_params=None, hook_clas num_procs=self.num_procs, hook_class=hook_class, fault_stuff=fault_stuff, - Tend=self.get_Tend(), + Tend=self.get_Tend() if Tend is None else Tend, custom_controller_params=custom_controller_params, - custom_problem_params=custom_problem_params, space_comm=space_comm, ) - def compare_strategies(self, run=0, faults=False, ax=None): + def compare_strategies(self, run=0, faults=False, ax=None): # pragma: no cover ''' Take a closer look at how the strategies compare for a specific run @@ -771,7 +365,9 @@ def compare_strategies(self, run=0, faults=False, ax=None): fig.tight_layout() plt.savefig(f'data/{self.get_name()}-comparison.pdf', transparent=True) - def scrutinize_visual(self, strategy, run, faults, ax=None, k_ax=None, ls='-', plot_restarts=False): + def scrutinize_visual( + self, strategy, run, faults, ax=None, k_ax=None, ls='-', plot_restarts=False + ): # pragma: no cover ''' Take a closer look at a specific run with a plot @@ -792,7 +388,7 @@ def scrutinize_visual(self, strategy, run, faults, ax=None, k_ax=None, ls='-', p else: store = False - force_params = dict() + force_params = {} stats, controller, Tend = self.single_run( strategy=strategy, @@ -843,7 +439,7 @@ def scrutinize(self, strategy, run, faults=True): Returns: None ''' - force_params = dict() + force_params = {} force_params['controller_params'] = {'logger_level': 15} stats, controller, Tend = self.single_run(strategy=strategy, run=run, faults=faults, force_params=force_params) @@ -856,7 +452,7 @@ def scrutinize(self, strategy, run, faults=True): print(f'\nOverview for {strategy.name} strategy') # see if we can determine if the faults where recovered - no_faults = self.load(strategy, False) + no_faults = self.load(strategy=strategy, faults=False) e_star = np.mean(no_faults.get('error', [0])) if t < Tend: error = np.inf @@ -889,10 +485,10 @@ def scrutinize(self, strategy, run, faults=True): # print faults faults = get_sorted(stats, type='bitflip') print('\nfaults:') - print(' t | level | iter | node | bit | trgt | pos') - print('------+-------+------+------+-----+------+----') + print(' t | level | iter | node | bit | trgt | pos') + print('--------+-------+------+------+-----+------+----') for f in faults: - print(f' {f[0]:.2f} | {f[1][0]:5d} | {f[1][1]:4d} | {f[1][2]:4d} | {f[1][4]:3d} | {f[1][5]:4d} |', f[1][3]) + print(f' {f[0]:6.2f} | {f[1][0]:5d} | {f[1][1]:4d} | {f[1][2]:4d} | {f[1][4]:3d} | {f[1][5]:4d} |', f[1][3]) return None @@ -919,7 +515,7 @@ def convert_faults(self, faults): bit = [faults[i][1][4] for i in range(len(faults))] return time, level, iteration, node, problem_pos, bit - def get_path(self, strategy, faults): + def get_path(self, **kwargs): ''' Get the path to where the stats are stored @@ -930,9 +526,9 @@ def get_path(self, strategy, faults): Returns: str: The path to what you are looking for ''' - return f'{self.stats_path}/{self.get_name(strategy, faults)}.pickle' + return f'{self.stats_path}/{self.get_name(**kwargs)}.pickle' - def get_name(self, strategy=None, faults=False): + def get_name(self, strategy=None, faults=True, mode=None): ''' Function to get a unique name for a kind of statistics based on the problem and strategy that was used @@ -953,7 +549,7 @@ def get_name(self, strategy=None, faults=False): prob_name = 'Lorenz' elif self.prob == run_Schroedinger: prob_name = 'Schroedinger' - elif self.prob == run_leaky_superconductor: + elif self.prob == run_quench: prob_name = 'Quench' else: raise NotImplementedError(f'Name not implemented for problem {self.prob}') @@ -968,42 +564,41 @@ def get_name(self, strategy=None, faults=False): else: strategy_name = '' - return f'{prob_name}{strategy_name}{fault_name}-{self.num_procs}procs' + mode = self.mode if mode is None else mode + if mode == 'regular': + mode_thing = f'-regular{self.kwargs["fault_frequency_iter"] if faults else ""}' + else: + mode_thing = '' + + return f'{prob_name}{strategy_name}{fault_name}-{self.num_procs}procs{mode_thing}' - def store(self, strategy, faults, dat): + def store(self, dat, **kwargs): ''' Stores the data for a run at a predefined path Args: - strategy (Strategy): Resilience strategy - faults (bool): Whether or not faults where inserted dat (dict): The data of the recorded statistics Returns: None ''' - with open(self.get_path(strategy, faults), 'wb') as f: + with open(self.get_path(**kwargs), 'wb') as f: pickle.dump(dat, f) return None - def load(self, strategy=None, faults=True): + def load(self, **kwargs): ''' Loads the stats belonging to a specific strategy and whether or not faults where inserted. When no data has been generated yet, a dictionary is returned which only contains the number of completed runs, which is 0 of course. - Args: - strategy (Strategy): Resilience strategy - faults (bool): Whether or not faults where inserted - Returns: dict: Data from previous run or if it is not available a placeholder dictionary ''' - if strategy is None: - strategy = self.strategies[MPI.COMM_WORLD.rank % len(self.strategies)] + kwargs['strategy'] = kwargs.get('strategy', self.strategies[MPI.COMM_WORLD.rank % len(self.strategies)]) try: - with open(self.get_path(strategy, faults), 'rb') as f: + with open(self.get_path(**kwargs), 'rb') as f: dat = pickle.load(f) except FileNotFoundError: return {'runs': 0} @@ -1016,27 +611,26 @@ def get_thresh(self, strategy=None): Args: strategy (Strategy): The resilience strategy """ - fault_free = self.load(strategy, False) + fault_free = self.load(strategy=strategy, faults=False) assert fault_free['error'].std() / fault_free['error'].mean() < 1e-5 return self.recovery_thresh_abs + self.recovery_thresh * fault_free["error"].mean() - def get_recovered(self, strategy=None): + def get_recovered(self, **kwargs): ''' Determine the recovery rate for a specific strategy and store it to disk. - Args: - strategy (Strategy): The resilience strategy you want to get the recovery rate for. If left at None, it will - be computed for all available strategies - Returns: None ''' - if strategy is None: - [self.get_recovered(strat) for strat in self.strategies] - - with_faults = self.load(strategy, True) - with_faults['recovered'] = with_faults['error'] < self.get_thresh(strategy) - self.store(strategy, True, with_faults) + if 'strategy' not in kwargs.keys(): + [self.get_recovered(strategy=strat, **kwargs) for strat in self.strategies] + else: + try: + with_faults = self.load(faults=True, **kwargs) + with_faults['recovered'] = with_faults['error'] < self.get_thresh(kwargs['strategy']) + self.store(faults=True, dat=with_faults, **kwargs) + except KeyError: + print("Can\'t compute recovery rate right now") return None @@ -1112,7 +706,9 @@ def extra_mean(self, dat, no_faults, thingA, mask): else: return None - def plot_thingA_per_thingB(self, strategy, thingA, thingB, ax=None, mask=None, recovered=False, op=None): + def plot_thingA_per_thingB( + self, strategy, thingA, thingB, ax=None, mask=None, recovered=False, op=None + ): # pragma: no cover ''' Plot thingA vs. thingB for a single strategy @@ -1129,8 +725,8 @@ def plot_thingA_per_thingB(self, strategy, thingA, thingB, ax=None, mask=None, r None ''' op = self.rec_rate if op is None else op - dat = self.load(strategy, True) - no_faults = self.load(strategy, False) + dat = self.load(strategy=strategy, faults=True) + no_faults = self.load(strategy=strategy, faults=False) if mask is None: mask = np.ones_like(dat[thingB], dtype=bool) @@ -1178,7 +774,7 @@ def plot_things_per_things( store=True, ax=None, fig=None, - ): + ): # pragma: no cover ''' Plot thingA vs thingB for multiple strategies @@ -1221,7 +817,7 @@ def plot_things_per_things( return None - def plot_recovery_thresholds(self, strategies=None, thresh_range=None, ax=None): + def plot_recovery_thresholds(self, strategies=None, thresh_range=None, ax=None): # pragma: no cover ''' Plot the recovery rate for a range of thresholds @@ -1243,8 +839,8 @@ def plot_recovery_thresholds(self, strategies=None, thresh_range=None, ax=None): for strategy_idx in range(len(strategies)): strategy = strategies[strategy_idx] # load the stats - fault_free = self.load(strategy, False) - with_faults = self.load(strategy, True) + fault_free = self.load(strategy=strategy, faults=False) + with_faults = self.load(strategy=strategy, faults=True) for thresh_idx in range(len(thresh_range)): rec_mask = with_faults['error'] < thresh_range[thresh_idx] * fault_free['error'].mean() @@ -1257,7 +853,7 @@ def plot_recovery_thresholds(self, strategies=None, thresh_range=None, ax=None): return None - def analyse_adaptivity(self, mask): + def analyse_adaptivity(self, mask): # pragma: no cover ''' Analyse a set of runs with adaptivity @@ -1286,7 +882,7 @@ def analyse_adaptivity(self, mask): print(f'We only restart when e_em > e_tol = {e_tol:.2e}!') return None - def analyse_adaptivity_single(self, run): + def analyse_adaptivity_single(self, run): # pragma: no cover ''' Compute what the difference in embedded and global error are for a specific run with adaptivity @@ -1322,7 +918,7 @@ def analyse_adaptivity_single(self, run): return [e_em[i][-1] for i in [0, 1]], e_glob - def analyse_HotRod(self, mask): + def analyse_HotRod(self, mask): # pragma: no cover ''' Analyse a set of runs with Hot Rod @@ -1358,7 +954,7 @@ def analyse_HotRod(self, mask): print(f'We only restart when diff > tol = {tol:.2e}!') return None - def analyse_HotRod_single(self, run): + def analyse_HotRod_single(self, run): # pragma: no cover ''' Compute what the difference in embedded, extrapolated and global error are for a specific run with Hot Rod @@ -1400,7 +996,7 @@ def analyse_HotRod_single(self, run): return [e_em[i][-1] for i in [0, 1]], [e_ex[i][-1] for i in [0, 1]], e_glob - def print_faults(self, mask=None): + def print_faults(self, mask=None): # pragma: no cover ''' Print all faults that happened within a certain mask @@ -1441,16 +1037,16 @@ def get_mask(self, strategy=None, key=None, val=None, op='eq', old_mask=None, co Numpy.ndarray with boolean entries that can be used as a mask ''' strategy = self.strategies[0] if strategy is None else strategy - dat = self.load(strategy, True) + dat = self.load(strategy=strategy, faults=True) if compare_faults: if val is not None: raise ValueError('Can\'t use val and compare_faults in get_mask at the same time!') else: - vals = self.load(strategy, False)[key] + vals = self.load(strategy=strategy, faults=False)[key] val = sum(vals) / len(vals) - if None in [key, val] and not op in ['isfinite']: + if None in [key, val] and op not in ['isfinite']: mask = dat['bit'] == dat['bit'] else: if op == 'uneq': @@ -1485,7 +1081,9 @@ def get_fixable_faults_only(self, strategy): Returns: Numpy.ndarray with boolean entries that can be used as a mask """ - fixable = strategy.get_fixable_params(maxiter=self.get_custom_description()['step_params']['maxiter']) + fixable = strategy.get_fixable_params( + maxiter=strategy.get_custom_description(self.prob, self.num_procs)['step_params']['maxiter'] + ) mask = self.get_mask(strategy=strategy) for kwargs in fixable: @@ -1509,7 +1107,7 @@ def get_index(self, mask=None): else: return np.arange(len(mask))[mask] - def get_statistics_info(self, mask=None, strategy=None, print_all=False, ax=None): + def get_statistics_info(self, mask=None, strategy=None, print_all=False, ax=None): # pragma: no cover ''' Get information about how many data points for faults we have given a particular mask @@ -1526,7 +1124,7 @@ def get_statistics_info(self, mask=None, strategy=None, print_all=False, ax=None # load some data from which to infer the number occurrences of some event strategy = self.strategies[0] if strategy is None else strategy - dat = self.load(strategy, True) + dat = self.load(stratagy=strategy, faults=True) # make a dummy mask in case none is supplied if mask is None: @@ -1548,7 +1146,7 @@ def get_statistics_info(self, mask=None, strategy=None, print_all=False, ax=None return None - def combinations_histogram(self, dat=None, keys=None, mask=None, ax=None): + def combinations_histogram(self, dat=None, keys=None, mask=None, ax=None): # pragma: no cover ''' Make a histogram ouf of the occurrences of combinations @@ -1571,7 +1169,7 @@ def combinations_histogram(self, dat=None, keys=None, mask=None, ax=None): return ax - def get_combination_histogram(self, dat=None, keys=None, mask=None): + def get_combination_histogram(self, dat=None, keys=None, mask=None): # pragma: no cover ''' Check how many combinations of values we expect and how many we find to see if we need to do more experiments. It is assumed that each allowed value for each key appears at least once in dat after the mask was applied @@ -1587,7 +1185,7 @@ def get_combination_histogram(self, dat=None, keys=None, mask=None): ''' # load default values - dat = self.load(self.strategies[0], True) if dat is None else dat + dat = self.load(strategy=self.strategies[0], faults=True) if dat is None else dat keys = ['iteration', 'bit', 'node'] if keys is None else keys if mask is None: mask = np.ones_like(dat['error'], dtype=bool) @@ -1673,7 +1271,7 @@ def count_occurrences(self, vals): def bar_plot_thing( self, x=None, thing=None, ax=None, mask=None, store=False, faults=False, name=None, op=None, args=None - ): + ): # pragma: no cover ''' Make a bar plot about something! @@ -1701,8 +1299,8 @@ def bar_plot_thing( strategy = self.strategies[strategy_idx] # load the values - dat = self.load(strategy, faults) - no_faults = self.load(strategy, False) + dat = self.load(strategy=strategy, faults=faults) + no_faults = self.load(strategy=strategy, faults=False) # check if we have a mask if mask is None: @@ -1730,21 +1328,98 @@ def bar_plot_thing( return None + def fault_frequency_plot(self, ax, iter_ax, kwargs_range, strategy=None): # pragma: no cover + func_args = locals() + func_args.pop('self', None) + if strategy is None: + for strat in self.strategies: + args = {**func_args, 'strategy': strat} + self.fault_frequency_plot(**args) + return None + + # load data + all_data = {} + for me in kwargs_range['fault_frequency_iter']: + self.kwargs['fault_frequency_iter'] = me + self.get_recovered() + all_data[me] = self.load(strategy=strategy, faults=True, mode='regular') + + # get_recovery_rate + results = {} + results['frequencies'] = list(all_data.keys()) + results['recovery_rate'] = [ + len(all_data[key]['recovered'][all_data[key]['recovered'] is True]) / len(all_data[key]['recovered']) + for key in all_data.keys() + ] + # results['iterations'] = [np.mean(all_data[key]['total_iteration']) for key in all_data.keys()] + results['iterations'] = [ + np.mean(all_data[key]['total_iteration'][all_data[key]['error'] != np.inf]) for key in all_data.keys() + ] + + ax.plot(results['frequencies'], results['recovery_rate'], **strategy.style) + iter_ax.plot(results['frequencies'], results['iterations'], **{**strategy.style, 'ls': '--'}) + + ax.set_xscale('log') + ax.set_xlabel('iterations between fault') + ax.set_ylabel('recovery rate') + iter_ax.set_ylabel('average total iterations if not crashed (dashed)') + ax.legend(frameon=False) + + return None + + +def check_local_error(): # pragma: no cover + """ + Make a plot of the resolution over time for all problems + """ + problems = [run_vdp, run_Lorenz, run_Schroedinger, run_quench] + problems = [run_quench] + strategies = [BaseStrategy(), AdaptivityStrategy(), IterateStrategy()] + + for i in range(len(problems)): + stats_analyser = FaultStats( + prob=problems[i], + strategies=strategies, + faults=[False], + reload=True, + recovery_thresh=1.1, + num_procs=1, + mode='random', + ) + stats_analyser.compare_strategies() + plt.show() + def main(): stats_analyser = FaultStats( - prob=run_leaky_superconductor, + prob=run_vdp, strategies=[BaseStrategy(), AdaptivityStrategy(), IterateStrategy(), HotRodStrategy()], faults=[False, True], reload=True, recovery_thresh=1.1, - recovery_thresh_abs=5e-5, + # recovery_thresh_abs=1e-5, num_procs=1, mode='random', stats_path='data/stats-jusuf', ) + ######################## + # msk = stats_analyser.get_mask(AdaptivityStrategy(), val=False, key='recovered') + # stats_analyser.print_faults(msk) + fig, ax = plt.subplots() + iter_ax = ax.twinx() + kwargs_range = {'fault_frequency_iter': (10, 100, 1000, 10000)} + stats_analyser.run_stats_generation(runs=10, kwargs_range=kwargs_range) + stats_analyser.fault_frequency_plot(ax=ax, iter_ax=iter_ax, kwargs_range=kwargs_range) + # stats_analyser.scrutinize(AdaptivityStrategy(), 4, True) + plt.show() + return None + ######################## + + stats_analyser.run_stats_generation(runs=5000) + + if MPI.COMM_WORLD.rank > 0: # make sure only one rank accesses the data + return None - stats_analyser.run_stats_generation(runs=1000) stats_analyser.get_recovered() mask = None @@ -1811,4 +1486,5 @@ def main(): if __name__ == "__main__": + # check_local_error() main() diff --git a/pySDC/projects/Resilience/heat.py b/pySDC/projects/Resilience/heat.py index 8d31cb2d6d..02af22c418 100644 --- a/pySDC/projects/Resilience/heat.py +++ b/pySDC/projects/Resilience/heat.py @@ -7,6 +7,7 @@ from pySDC.helpers.stats_helper import get_sorted from pySDC.projects.Resilience.hook import hook_collection, LogData import numpy as np +from pySDC.projects.Resilience.strategies import merge_descriptions def run_heat( @@ -16,7 +17,6 @@ def run_heat( hook_class=LogData, fault_stuff=None, custom_controller_params=None, - custom_problem_params=None, ): """ Run a heat problem with default parameters. @@ -28,7 +28,6 @@ def run_heat( hook_class (pySDC.Hook): A hook to store data fault_stuff (dict): A dictionary with information on how to add faults custom_controller_params (dict): Overwrite presets - custom_problem_params (dict): Overwrite presets Returns: dict: The stats object @@ -58,9 +57,6 @@ def run_heat( 'liniter': None, } - if custom_problem_params is not None: - problem_params = {**problem_params, **custom_problem_params} - # initialize step parameters step_params = dict() step_params['maxiter'] = 5 @@ -84,11 +80,7 @@ def run_heat( description['step_params'] = step_params if custom_description is not None: - for k in custom_description.keys(): - if k == 'sweeper_class': - description[k] = custom_description[k] - continue - description[k] = {**description.get(k, {}), **custom_description.get(k, {})} + description = merge_descriptions(description, custom_description) # set time parameters t0 = 0.0 diff --git a/pySDC/projects/Resilience/jobscript_jusuf.sh b/pySDC/projects/Resilience/jobscript_jusuf.sh new file mode 100755 index 0000000000..639a564beb --- /dev/null +++ b/pySDC/projects/Resilience/jobscript_jusuf.sh @@ -0,0 +1,22 @@ +#!/bin/bash +#SBATCH --nodes=1 +#SBATCH --ntasks=64 +#SBATCH --time=24:00:00 +#SBATCH --output=out/out%j.txt +#SBATCH --error=out/err%j.txt +#SBATCH -A cstma +#SBATCH --mail-type=END,FAIL +#SBATCH --mail-user=t.baumann@fz-juelich.de +#SBATCH -p batch +#SBATCH -J red_wedding + +module --force purge +module load Stages/2023 +module load Intel/2022.1.0 +module load ParaStationMPI/5.8.0-1 + +cd /p/project/ccstma/baumann7/pySDC/pySDC/projects/Resilience + +source /p/project/ccstma/baumann7/miniconda/bin/activate pySDC + +srun -n 64 python ${1} diff --git a/pySDC/projects/Resilience/leaky_superconductor.py b/pySDC/projects/Resilience/leaky_superconductor.py deleted file mode 100644 index a2d0dde360..0000000000 --- a/pySDC/projects/Resilience/leaky_superconductor.py +++ /dev/null @@ -1,273 +0,0 @@ -# script to run a quench problem -from pySDC.implementations.problem_classes.LeakySuperconductor import LeakySuperconductor, LeakySuperconductorIMEX -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 -import numpy as np - -import matplotlib.pyplot as plt -from pySDC.core.Errors import ConvergenceError - - -class live_plot(hooks): # pragma: no cover - """ - This hook plots the solution and the non-linear part of the right hand side after every step. Keep in mind that using adaptivity will result in restarts, which is not marked in these plots. Prepare to see the temperature profile jumping back again after a restart. - """ - - def _plot_state(self, step, level_number): # pragma: no cover - """ - Plot the solution at all collocation nodes and the non-linear part of the right hand side - - Args: - step (pySDC.Step.step): The current step - level_number (int): Number of current level - - Returns: - None - """ - L = step.levels[level_number] - for ax in self.axs: - ax.cla() - [self.axs[0].plot(L.prob.xv, L.u[i], legend=f"node {i}") for i in range(len(L.u))] - self.axs[0].axhline(L.prob.u_thresh, color='black') - self.axs[1].plot(L.prob.xv, L.prob.eval_f_non_linear(L.u[-1], L.time)) - self.axs[0].set_ylim(0, 0.025) - self.fig.suptitle(f"t={L.time:.2e}, k={step.status.iter}") - plt.pause(1e-9) - - def pre_run(self, step, level_number): # pragma: no cover - """ - Setup a figure to plot into - - Args: - step (pySDC.Step.step): The current step - level_number (int): Number of current level - - Returns: - None - """ - self.fig, self.axs = plt.subplots(1, 2, figsize=(10, 4)) - - def post_step(self, step, level_number): # pragma: no cover - """ - Call the plotting function after the step - - Args: - step (pySDC.Step.step): The current step - level_number (int): Number of current level - - Returns: - None - """ - self._plot_state(step, level_number) - - -def run_leaky_superconductor( - custom_description=None, - num_procs=1, - Tend=6e2, - hook_class=LogData, - fault_stuff=None, - custom_controller_params=None, - custom_problem_params=None, - imex=False, - u0=None, - t0=None, - **kwargs, -): - """ - Run a toy problem of a superconducting magnet with a temperature leak with default parameters. - - Args: - custom_description (dict): Overwrite presets - num_procs (int): Number of steps for MSSDC - Tend (float): Time to integrate to - hook_class (pySDC.Hook): A hook to store data - fault_stuff (dict): A dictionary with information on how to add faults - custom_controller_params (dict): Overwrite presets - custom_problem_params (dict): Overwrite presets - imex (bool): Solve the problem IMEX or fully implicit - u0 (dtype_u): Initial value - t0 (float): Starting time - - Returns: - dict: The stats object - controller: The controller - Tend: The time that was supposed to be integrated to - """ - - # initialize level parameters - level_params = {} - level_params['dt'] = 10.0 - - # initialize sweeper parameters - sweeper_params = {} - sweeper_params['quad_type'] = 'RADAU-RIGHT' - sweeper_params['num_nodes'] = 3 - sweeper_params['QI'] = 'IE' - sweeper_params['QE'] = 'PIC' - - problem_params = {} - - if custom_problem_params is not None: - problem_params = {**problem_params, **custom_problem_params} - - # initialize step parameters - step_params = {} - step_params['maxiter'] = 5 - - # initialize controller parameters - controller_params = {} - controller_params['logger_level'] = 30 - controller_params['hook_class'] = hook_collection + (hook_class if type(hook_class) == list else [hook_class]) - controller_params['mssdc_jac'] = False - - if custom_controller_params is not None: - controller_params = {**controller_params, **custom_controller_params} - - # fill description dictionary for easy step instantiation - description = {} - description['problem_class'] = LeakySuperconductorIMEX if imex else LeakySuperconductor - description['problem_params'] = problem_params - description['sweeper_class'] = imex_1st_order if imex else generic_implicit - description['sweeper_params'] = sweeper_params - description['level_params'] = level_params - description['step_params'] = step_params - - if custom_description is not None: - for k in custom_description.keys(): - if k == 'sweeper_class': - description[k] = custom_description[k] - continue - description[k] = {**description.get(k, {}), **custom_description.get(k, {})} - - # set time parameters - t0 = 0.0 if t0 is None else t0 - - # instantiate controller - controller = controller_nonMPI(num_procs=num_procs, controller_params=controller_params, description=description) - - # insert faults - if fault_stuff is not None: - from pySDC.projects.Resilience.fault_injection import prepare_controller_for_faults - - rnd_args = {'iteration': 5, 'min_node': 1} - args = {'time': 21.0, 'target': 0} - prepare_controller_for_faults(controller, fault_stuff, rnd_args, args) - - # get initial values on finest level - P = controller.MS[0].levels[0].prob - uinit = P.u_exact(t0) if u0 is None else u0 - - # call main function to get things done... - try: - uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) - except ConvergenceError: - stats = controller.return_stats() - return stats, controller, Tend - - -def plot_solution(stats, controller): # pragma: no cover - import matplotlib.pyplot as plt - - fig, ax = plt.subplots(1, 1) - u_ax = ax - dt_ax = u_ax.twinx() - - u = get_sorted(stats, type='u', recomputed=False) - u_ax.plot([me[0] for me in u], [max(me[1]) for me in u], label=r'$T$') - - dt = get_sorted(stats, type='dt', recomputed=False) - dt_ax.plot([me[0] for me in dt], [me[1] for me in dt], color='black', ls='--') - u_ax.plot([None], [None], color='black', ls='--', label=r'$\Delta t$') - - P = controller.MS[0].levels[0].prob - u_ax.axhline(P.params.u_thresh, color='grey', ls='-.', label=r'$T_\mathrm{thresh}$') - u_ax.axhline(P.params.u_max, color='grey', ls=':', label=r'$T_\mathrm{max}$') - - u_ax.legend() - u_ax.set_xlabel(r'$t$') - u_ax.set_ylabel(r'$T$') - dt_ax.set_ylabel(r'$\Delta t$') - - -def compare_imex_full(plotting=False, leak_type='linear'): - """ - Compare the results of IMEX and fully implicit runs. For IMEX we need to limit the step size in order to achieve convergence, but for fully implicit, adaptivity can handle itself better. - - Args: - plotting (bool): Plot the solution or not - """ - from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity - from pySDC.implementations.hooks.log_work import LogWork - from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun - - maxiter = 5 - num_nodes = 3 - newton_iter_max = 20 - - res = {} - rhs = {} - error = {} - - custom_description = {} - custom_description['problem_params'] = { - 'newton_tol': 1e-10, - 'newton_iter': newton_iter_max, - 'nvars': 2**9, - 'leak_type': leak_type, - } - custom_description['step_params'] = {'maxiter': maxiter} - custom_description['sweeper_params'] = {'num_nodes': num_nodes} - - custom_controller_params = {'logger_level': 30} - for imex in [False, True]: - custom_description['convergence_controllers'] = {Adaptivity: {'e_tol': 1e-6, 'dt_max': 1e2}} - stats, controller, _ = run_leaky_superconductor( - custom_description=custom_description, - custom_controller_params=custom_controller_params, - imex=imex, - Tend=4.3e2, - hook_class=[LogWork, LogGlobalErrorPostRun], - ) - - res[imex] = get_sorted(stats, type='u')[-1][1] - newton_iter = [me[1] for me in get_sorted(stats, type='work_newton')] - rhs[imex] = np.mean([me[1] for me in get_sorted(stats, type='work_rhs')]) // 1 - error[imex] = get_sorted(stats, type='e_global_post_run')[-1][1] - - if imex: - assert all(me == 0 for me in newton_iter), "IMEX is not supposed to do Newton iterations!" - else: - assert ( - max(newton_iter) / num_nodes / maxiter <= newton_iter_max - ), "Took more Newton iterations than allowed!" - if plotting: # pragma: no cover - plot_solution(stats, controller) - - diff = abs(res[True] - res[False]) - thresh = 3e-3 - assert ( - diff < thresh - ), f"Difference between IMEX and fully-implicit too large! Got {diff:.2e}, allowed is only {thresh:.2e}!" - prob = controller.MS[0].levels[0].prob - assert ( - max(res[True]) > prob.u_max - ), f"Expected runaway to happen, but maximum temperature is {max(res[True]):.2e} < u_max={prob.u_max:.2e}!" - - assert ( - rhs[True] == rhs[False] - ), f"Expected IMEX and fully implicit schemes to take the same number of right hand side evaluations per step, but got {rhs[True]} and {rhs[False]}!" - - assert ( - error[True] > error[False] - ), f"Expected IMEX to be less accurate at the same precision settings than unsplit version, got for IMEX: e={error[True]:.2e} and fully implicit: e={error[False]:.2e}" - assert error[True] < 1.1e-4, f'Expected error of IMEX version to be less than 1.1e-4, but got e={error[True]:.2e}!' - - -if __name__ == '__main__': - compare_imex_full(plotting=True) - plt.show() diff --git a/pySDC/projects/Resilience/notes/Lorenz/compare_strategies.png b/pySDC/projects/Resilience/notes/Lorenz/compare_strategies.png index 1aa165f58f..7921043db4 100644 Binary files a/pySDC/projects/Resilience/notes/Lorenz/compare_strategies.png and b/pySDC/projects/Resilience/notes/Lorenz/compare_strategies.png differ diff --git a/pySDC/projects/Resilience/notes/Lorenz/recovery_rate_compared.png b/pySDC/projects/Resilience/notes/Lorenz/recovery_rate_compared.png index 9b24184fca..893ad19bd4 100644 Binary files a/pySDC/projects/Resilience/notes/Lorenz/recovery_rate_compared.png and b/pySDC/projects/Resilience/notes/Lorenz/recovery_rate_compared.png differ diff --git a/pySDC/projects/Resilience/paper_plots.py b/pySDC/projects/Resilience/paper_plots.py index c1949560b8..7f95c70009 100644 --- a/pySDC/projects/Resilience/paper_plots.py +++ b/pySDC/projects/Resilience/paper_plots.py @@ -11,7 +11,7 @@ run_Lorenz, run_Schroedinger, run_vdp, - run_leaky_superconductor, + run_quench, ) from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal from pySDC.helpers.stats_helper import get_sorted @@ -23,7 +23,7 @@ BASE_PATH = 'data/paper' -def get_stats(problem, path='data/stats'): +def get_stats(problem, path='data/stats-jusuf'): """ Create a FaultStats object for a given problem to use for the plots. Note that the statistics need to be already generated somewhere else, this function will only load them. @@ -41,14 +41,14 @@ def get_stats(problem, path='data/stats'): mode = 'random' recovery_thresh_abs = { - run_leaky_superconductor: 5e-5, + run_quench: 5e-3, } strategies = [BaseStrategy(), AdaptivityStrategy(), IterateStrategy()] if JOURNAL not in ['JSC_beamer']: strategies += [HotRodStrategy()] - return FaultStats( + stats_analyser = FaultStats( prob=problem, strategies=strategies, faults=[False, True], @@ -59,6 +59,8 @@ def get_stats(problem, path='data/stats'): mode=mode, stats_path=path, ) + stats_analyser.get_recovered() + return stats_analyser def my_setup_mpl(**kwargs): @@ -66,17 +68,19 @@ def my_setup_mpl(**kwargs): mpl.rcParams.update({'lines.markersize': 6}) -def savefig(fig, name, format='pdf'): # pragma: no cover +def savefig(fig, name, format='pdf', tight_layout=True): # pragma: no cover """ Save a figure to some predefined location. Args: fig (Matplotlib.Figure): The figure of the plot name (str): The name of the plot + tight_layout (bool): Apply tight layout or leave as is Returns: None """ - fig.tight_layout() + if tight_layout: + fig.tight_layout() path = f'{BASE_PATH}/{name}.{format}' fig.savefig(path, bbox_inches='tight', transparent=True, dpi=200) print(f'saved "{path}"') @@ -95,6 +99,18 @@ def analyse_resilience(problem, path='data/stats', **kwargs): # pragma: no cove """ stats_analyser = get_stats(problem, path) + stats_analyser.get_recovered() + + strategy = IterateStrategy() + not_fixed = stats_analyser.get_mask(strategy=strategy, key='recovered', val=False) + 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) @@ -131,7 +147,7 @@ def plot_recovery_rate(stats_analyser, **kwargs): # pragma: no cover stats_analyser.plot_things_per_things( 'recovered', 'bit', False, op=stats_analyser.rec_rate, args={'ylabel': 'recovery rate'}, ax=axs[0] ) - plot_recovery_rate_recoverable_only(stats_analyser, fig, axs[1], ylabel='', xlabel='') + plot_recovery_rate_recoverable_only(stats_analyser, fig, axs[1], ylabel='') axs[1].get_legend().remove() axs[0].set_title('All faults') axs[1].set_title('Only recoverable faults') @@ -177,15 +193,15 @@ def compare_recovery_rate_problems(): # pragma: no cover stats = [ get_stats(run_vdp), get_stats(run_Lorenz), - get_stats(run_Schroedinger, 'data/stats-jusuf'), - get_stats(run_leaky_superconductor, 'data/stats-jusuf'), + get_stats(run_Schroedinger), + get_stats(run_quench), ] titles = ['Van der Pol', 'Lorenz attractor', r'Schr\"odinger', 'Quench'] my_setup_mpl() - fig, axs = plt.subplots(2, 2, figsize=figsize_by_journal(JOURNAL, 1, 0.7), sharey=True) + fig, axs = plt.subplots(2, 2, figsize=figsize_by_journal(JOURNAL, 1, 0.8), sharey=True) [ - plot_recovery_rate_recoverable_only(stats[i], fig, axs.flatten()[i], ylabel='', xlabel='', title=titles[i]) + plot_recovery_rate_recoverable_only(stats[i], fig, axs.flatten()[i], ylabel='', title=titles[i]) for i in range(len(stats)) ] @@ -193,13 +209,53 @@ def compare_recovery_rate_problems(): # pragma: no cover ax.get_legend().remove() axs[1, 1].legend(frameon=False) - axs[1, 0].set_xlabel('bit') axs[1, 0].set_ylabel('recovery rate') + axs[0, 0].set_ylabel('recovery rate') savefig(fig, 'compare_equations') -def plot_efficiency_polar(problem, path='data/stats'): # pragma: no cover +def plot_efficiency_polar_vdp(problem, path='data/stats'): # pragma: no cover + stats_analyser = get_stats(problem, path) + fig, ax = plt.subplots( + subplot_kw={'projection': 'polar'}, figsize=figsize_by_journal(JOURNAL, 0.7, 0.5), layout='constrained' + ) + theta, norms = plot_efficiency_polar_single(stats_analyser, ax) + + labels = ['fail rate', 'extra iterations\nfor recovery', 'iterations for solution'] + ax.set_xticks(theta[:-1], [f'{labels[i]}\nmax={norms[i]:.2f}' for i in range(len(labels))]) + ax.set_rlabel_position(90) + + fig.legend(frameon=False, loc='outside right', ncols=1) + savefig(fig, 'efficiency', tight_layout=False) + + +def plot_efficiency_polar_other(): # pragma: no cover + problems = [run_Lorenz, run_Schroedinger, run_quench] + paths = ['./data/stats/', './data/stats-jusuf', './data/stats-jusuf'] + titles = ['Lorenz attractor', r'Schr\"odinger', 'Quench'] + + fig, axs = plt.subplots( + 1, 3, subplot_kw={'projection': 'polar'}, figsize=figsize_by_journal(JOURNAL, 0.7, 0.5), layout='constrained' + ) + + for i in range(len(problems)): + stats_analyser = get_stats(problems[i], paths[i]) + ax = axs[i] + theta, norms = plot_efficiency_polar_single(stats_analyser, ax) + + labels = ['fail rate', 'extra iterations\nfor recovery', 'iterations for solution'] + ax.set_rlabel_position(90) + # ax.set_xticks(theta[:-1], [f'max={norms[i]:.2f}' for i in range(len(labels))]) + ax.set_xticks(theta[:-1], ['' for i in range(len(labels))]) + ax.set_title(titles[i]) + + handles, labels = fig.get_axes()[0].get_legend_handles_labels() + fig.legend(handles=handles, labels=labels, frameon=False, loc='outside lower center', ncols=4) + savefig(fig, 'efficiency_other', tight_layout=False) + + +def plot_efficiency_polar_single(stats_analyser, ax): # pragma: no cover """ Plot the recovery rate and the computational cost in a polar plot. @@ -218,12 +274,10 @@ def plot_efficiency_polar(problem, path='data/stats'): # pragma: no cover Returns: None """ - - stats_analyser = get_stats(problem, path) + # TODO: fix docs mask = stats_analyser.get_mask() # get empty mask, potentially put in some other mask later my_setup_mpl() - fig, ax = plt.subplots(subplot_kw={'projection': 'polar'}, figsize=(7 * cm, 7 * cm)) res = {} for strategy in stats_analyser.strategies: @@ -255,13 +309,7 @@ def plot_efficiency_polar(problem, path='data/stats'): # pragma: no cover theta = np.array([30, 150, 270, 30]) * 2 * np.pi / 360 for s in stats_analyser.strategies: ax.plot(theta, res_norm[s.name] + [res_norm[s.name][0]], label=s.label, color=s.color, marker=s.marker) - - labels = ['fail rate', 'extra iterations\nfor recovery', 'iterations for solution'] - ax.set_xticks(theta[:-1], [f'{labels[i]}\nmax={norms[i]:.2f}' for i in range(len(labels))]) - ax.set_rlabel_position(90) - - ax.legend(frameon=True, loc='lower right') - savefig(fig, 'efficiency') + return theta, norms def plot_adaptivity_stuff(): # pragma: no cover @@ -272,7 +320,9 @@ def plot_adaptivity_stuff(): # pragma: no cover Returns: None """ - from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import EstimateEmbeddedErrorNonMPI + from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import EstimateEmbeddedError + from pySDC.implementations.hooks.log_errors import LogLocalErrorPostStep + from pySDC.projects.Resilience.hook import LogData stats_analyser = get_stats(run_vdp, 'data/stats') @@ -293,26 +343,31 @@ def plot_error(stats, ax, iter_ax, strategy, **kwargs): Returns: None """ - e = get_sorted(stats, type='error_embedded_estimate', recomputed=False) - ax.plot([me[0] for me in e], [me[1] for me in e], markevery=15, **strategy.style, **kwargs) + markevery = 40 + e = get_sorted(stats, type='e_local_post_step', recomputed=False) + ax.plot([me[0] for me in e], [me[1] for me in e], markevery=markevery, **strategy.style, **kwargs) k = get_sorted(stats, type='k') - iter_ax.plot([me[0] for me in k], np.cumsum([me[1] for me in k]), **strategy.style, markevery=15, **kwargs) + iter_ax.plot( + [me[0] for me in k], np.cumsum([me[1] for me in k]), **strategy.style, markevery=markevery, **kwargs + ) ax.set_yscale('log') ax.set_ylabel('local error') - iter_ax.set_ylabel(r'iterations') + iter_ax.set_ylabel(r'SDC iterations') - force_params = {'convergence_controllers': {EstimateEmbeddedErrorNonMPI: {}}} + 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]: - stats, _, _ = stats_analyser.single_run(strategy=strategy(), force_params=force_params) + stats, _, _ = stats_analyser.single_run( + strategy=strategy(), force_params=force_params, hook_class=[LogLocalErrorPostStep, LogData] + ) plot_error(stats, axs[1], axs[2], strategy()) - if strategy == AdaptivityStrategy: - u = get_sorted(stats, type='u') + 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[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[1].set_ylim(bottom=1e-9) axs[2].set_xlabel(r'$t$') axs[0].set_ylabel('solution') axs[2].legend(frameon=JOURNAL == 'JSC_beamer') @@ -347,15 +402,16 @@ def plot_fault_vdp(bit=0): # pragma: no cover ) my_setup_mpl() - fig, ax = plt.subplots(1, 1, figsize=(TEXTWIDTH * 3.0 / 4.0, 5 * cm)) + fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.8, 0.5)) colors = ['blue', 'red', 'magenta'] ls = ['--', '-'] - markers = ['*', '.', 'y'] + markers = ['*', '^'] do_faults = [False, True] superscripts = ['*', ''] subscripts = ['', 't', ''] - run = 779 + 12 * bit + run = 779 + 12 * bit # for faults in u_t + # run = 11 + 12 * bit # for faults in u for i in range(len(do_faults)): stats, controller, Tend = stats_analyser.single_run( @@ -366,15 +422,15 @@ def plot_fault_vdp(bit=0): # pragma: no cover ) u = get_sorted(stats, type='u') faults = get_sorted(stats, type='bitflip') - for j in range(len(u[0][1])): + for j in [0, 1]: ax.plot( [me[0] for me in u], [me[1][j] for me in u], ls=ls[i], color=colors[j], label=rf'$u^{{{superscripts[i]}}}_{{{subscripts[j]}}}$', - marker=markers[0], - markevery=15, + marker=markers[j], + markevery=60, ) for idx in range(len(faults)): ax.axvline(faults[idx][0], color='black', label='Fault', ls=':') @@ -383,11 +439,104 @@ def plot_fault_vdp(bit=0): # pragma: no cover ) ax.set_title(f'Fault in bit {faults[idx][1][4]}') - ax.legend(frameon=False) + ax.legend(frameon=True, loc='lower left') ax.set_xlabel(r'$t$') savefig(fig, f'fault_bit_{bit}') +def plot_quench_solution(): # pragma: no cover + """ + Plot the solution of Quench problem over time + + Returns: + None + """ + my_setup_mpl() + if JOURNAL == 'JSC_beamer': + fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 0.9)) + else: + fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 1.0, 0.45)) + + strategy = BaseStrategy() + + custom_description = strategy.get_custom_description(run_quench) + + 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') + + 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) + ax.axhline(prob.u_max, label='$T_\mathrm{max}$', ls=':', color='grey', zorder=-1) + + ax.set_xlabel(r'$t$') + ax.legend(frameon=False) + savefig(fig, 'quench_sol') + + +def plot_Lorenz_solution(): # pragma: no cover + """ + Plot the solution of Lorenz attractor problem over time + + Returns: + None + """ + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + from pySDC.projects.Resilience.strategies import AdaptivityStrategy + from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun + + strategy = AdaptivityStrategy() + + my_setup_mpl() + + fig = plt.figure() + ax = fig.add_subplot(projection='3d') + + # if JOURNAL == 'JSC_beamer': + # fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 0.9)) + # else: + # fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 1.0, 0.33)) + + custom_description = strategy.get_custom_description(run_Lorenz, 1) + custom_description['convergence_controllers'] = {Adaptivity: {'e_tol': 1e-10}} + + stats, _, _ = run_Lorenz( + custom_description=custom_description, + Tend=strategy.get_Tend(run_Lorenz, 1) * 20, + hook_class=LogGlobalErrorPostRun, + ) + + u = get_sorted(stats, type='u') + e = get_sorted(stats, type='e_global_post_run')[-1] + print(u[-1], e) + ax.plot([me[1][0] for me in u], [me[1][1] for me in u], [me[1][2] for me in u]) + + ################## + from pySDC.projects.Resilience.strategies import DIRKStrategy, ERKStrategy, IterateStrategy + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK + + strategy = ERKStrategy() + custom_description = strategy.get_custom_description(run_Lorenz, 1) + custom_description['convergence_controllers'] = {Adaptivity: {'e_tol': 1e-10}} + stats, _, _ = run_Lorenz( + custom_description=custom_description, + Tend=strategy.get_Tend(run_Lorenz, 1) * 20, + hook_class=LogGlobalErrorPostRun, + ) + + u = get_sorted(stats, type='u') + e = get_sorted(stats, type='e_global_post_run')[-1] + print(u[-1], e) + ax.plot([me[1][0] for me in u], [me[1][1] for me in u], [me[1][2] for me in u], ls='--') + ################ + ax.set_xlabel('x') + ax.set_ylabel('y') + ax.set_zlabel('z') + savefig(fig, 'lorenz_sol') + + def plot_vdp_solution(): # pragma: no cover """ Plot the solution of van der Pol problem over time to illustrate the varying time scales. @@ -398,12 +547,14 @@ def plot_vdp_solution(): # pragma: no cover from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity my_setup_mpl() - fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 0.9)) + if JOURNAL == 'JSC_beamer': + fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 0.9)) + else: + fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 1.0, 0.33)) custom_description = {'convergence_controllers': {Adaptivity: {'e_tol': 1e-7}}} - problem_params = {} - stats, _, _ = run_vdp(custom_description=custom_description, custom_problem_params=problem_params, Tend=28.6) + stats, _, _ = run_vdp(custom_description=custom_description, Tend=28.6) u = get_sorted(stats, type='u') ax.plot([me[0] for me in u], [me[1][0] for me in u], color='black') @@ -412,6 +563,80 @@ def plot_vdp_solution(): # pragma: no cover savefig(fig, 'vdp_sol') +def work_precision(): # pragma: no cover + from pySDC.projects.Resilience.work_precision import ( + all_problems, + single_problem, + ODEs, + get_fig, + execute_configurations, + save_fig, + get_configs, + MPI, + vdp_stiffness_plot, + ) + + all_params = { + 'record': False, + 'work_key': 't', + 'precision_key': 'e_global_rel', + 'plotting': True, + 'base_path': 'data/paper', + } + + for mode in ['compare_strategies', 'parallel_efficiency']: + 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}$') + + 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"], + ) + + vdp_stiffness_plot(base_path='data/paper') + + +def make_plots_for_TIME_X_website(): # pragma: no cover + global JOURNAL, BASE_PATH + JOURNAL = 'JSC_beamer' + BASE_PATH = 'data/paper/time-x_website' + + fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 2.0 / 3.0)) + plot_recovery_rate_recoverable_only(get_stats(run_vdp), fig, ax) + savefig(fig, 'recovery_rate', format='png') + + from pySDC.projects.Resilience.work_precision import vdp_stiffness_plot + + vdp_stiffness_plot(base_path=BASE_PATH, format='png') + + def make_plots_for_SIAM_CSE23(): # pragma: no cover """ Make plots for the SIAM talk @@ -437,12 +662,14 @@ def make_plots_for_paper(): # pragma: no cover JOURNAL = 'Springer_Numerical_Algorithms' BASE_PATH = 'data/paper' + plot_vdp_solution() + plot_quench_solution() plot_recovery_rate(get_stats(run_vdp)) plot_fault_vdp(0) plot_fault_vdp(13) plot_adaptivity_stuff() - plot_efficiency_polar(run_vdp) compare_recovery_rate_problems() + work_precision() def make_plots_for_notes(): # pragma: no cover @@ -454,9 +681,11 @@ def make_plots_for_notes(): # pragma: no cover BASE_PATH = 'notes/Lorenz' analyse_resilience(run_Lorenz, format='png') + analyse_resilience(run_quench, format='png') if __name__ == "__main__": - make_plots_for_notes() - make_plots_for_SIAM_CSE23() + # make_plots_for_notes() + # make_plots_for_SIAM_CSE23() + # make_plots_for_TIME_X_website() make_plots_for_paper() diff --git a/pySDC/projects/Resilience/piline.py b/pySDC/projects/Resilience/piline.py index e77eecb1f8..9c231cb01f 100644 --- a/pySDC/projects/Resilience/piline.py +++ b/pySDC/projects/Resilience/piline.py @@ -8,6 +8,7 @@ from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity from pySDC.implementations.convergence_controller_classes.hotrod import HotRod from pySDC.projects.Resilience.hook import LogData, hook_collection +from pySDC.projects.Resilience.strategies import merge_descriptions def run_piline( @@ -17,7 +18,6 @@ def run_piline( hook_class=LogData, fault_stuff=None, custom_controller_params=None, - custom_problem_params=None, ): """ Run a Piline problem with default parameters. @@ -29,7 +29,6 @@ def run_piline( hook_class (pySDC.Hook): A hook to store data fault_stuff (dict): A dictionary with information on how to add faults custom_controller_params (dict): Overwrite presets - custom_problem_params (dict): Overwrite presets Returns: dict: The stats object @@ -58,9 +57,6 @@ def run_piline( 'Rl': 5.0, } - if custom_problem_params is not None: - problem_params = {**problem_params, **custom_problem_params} - # initialize step parameters step_params = dict() step_params['maxiter'] = 4 @@ -84,11 +80,7 @@ def run_piline( description['step_params'] = step_params if custom_description is not None: - for k in custom_description.keys(): - if k == 'sweeper_class': - description[k] = custom_description[k] - continue - description[k] = {**description.get(k, {}), **custom_description.get(k, {})} + description = merge_descriptions(description, custom_description) # set time parameters t0 = 0.0 diff --git a/pySDC/projects/Resilience/quench.py b/pySDC/projects/Resilience/quench.py new file mode 100644 index 0000000000..d9447428f8 --- /dev/null +++ b/pySDC/projects/Resilience/quench.py @@ -0,0 +1,557 @@ +# 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 +import numpy as np + +import matplotlib.pyplot as plt +from pySDC.core.Errors import ConvergenceError + + +class live_plot(hooks): # pragma: no cover + """ + This hook plots the solution and the non-linear part of the right hand side after every step. Keep in mind that using adaptivity will result in restarts, which is not marked in these plots. Prepare to see the temperature profile jumping back again after a restart. + """ + + def _plot_state(self, step, level_number): # pragma: no cover + """ + Plot the solution at all collocation nodes and the non-linear part of the right hand side + + Args: + step (pySDC.Step.step): The current step + level_number (int): Number of current level + + Returns: + None + """ + L = step.levels[level_number] + for ax in self.axs: + ax.cla() + # [self.axs[0].plot(L.prob.xv, L.u[i], label=f"node {i}") for i in range(len(L.u))] + self.axs[0].plot(L.prob.xv, L.u[-1]) + self.axs[0].axhline(L.prob.u_thresh, color='black') + self.axs[1].plot(L.prob.xv, L.prob.eval_f_non_linear(L.u[-1], L.time)) + self.axs[0].set_ylim(0, 0.025) + self.fig.suptitle(f"t={L.time:.2e}, k={step.status.iter}") + plt.pause(1e-1) + + def pre_run(self, step, level_number): # pragma: no cover + """ + Setup a figure to plot into + + Args: + step (pySDC.Step.step): The current step + level_number (int): Number of current level + + Returns: + None + """ + self.fig, self.axs = plt.subplots(1, 2, figsize=(10, 4)) + + def post_step(self, step, level_number): # pragma: no cover + """ + Call the plotting function after the step + + Args: + step (pySDC.Step.step): The current step + level_number (int): Number of current level + + Returns: + None + """ + self._plot_state(step, level_number) + + +def run_quench( + custom_description=None, + num_procs=1, + Tend=6e2, + hook_class=LogData, + fault_stuff=None, + custom_controller_params=None, + imex=False, + u0=None, + t0=None, + use_MPI=False, + **kwargs, +): + """ + Run a toy problem of a superconducting magnet with a temperature leak with default parameters. + + Args: + custom_description (dict): Overwrite presets + num_procs (int): Number of steps for MSSDC + Tend (float): Time to integrate to + hook_class (pySDC.Hook): A hook to store data + fault_stuff (dict): A dictionary with information on how to add faults + custom_controller_params (dict): Overwrite presets + imex (bool): Solve the problem IMEX or fully implicit + u0 (dtype_u): Initial value + t0 (float): Starting time + use_MPI (bool): Whether or not to use MPI + + Returns: + dict: The stats object + controller: The controller + Tend: The time that was supposed to be integrated to + """ + + # initialize level parameters + level_params = {} + level_params['dt'] = 10.0 + + # initialize sweeper parameters + sweeper_params = {} + sweeper_params['quad_type'] = 'RADAU-RIGHT' + sweeper_params['num_nodes'] = 3 + sweeper_params['QI'] = 'IE' + sweeper_params['QE'] = 'PIC' + + problem_params = { + 'newton_tol': 1e-9, + } + + # initialize step parameters + step_params = {} + step_params['maxiter'] = 5 + + # initialize controller parameters + controller_params = {} + controller_params['logger_level'] = 30 + controller_params['hook_class'] = hook_collection + (hook_class if type(hook_class) == list else [hook_class]) + controller_params['mssdc_jac'] = False + + if custom_controller_params is not None: + controller_params = {**controller_params, **custom_controller_params} + + # fill description dictionary for easy step instantiation + 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_params'] = sweeper_params + description['level_params'] = level_params + description['step_params'] = step_params + + if custom_description is not None: + description = merge_descriptions(description, custom_description) + + # set time parameters + t0 = 0.0 if t0 is None else t0 + + # instantiate controller + controller_args = { + 'controller_params': controller_params, + 'description': description, + } + if use_MPI: + from mpi4py import MPI + from pySDC.implementations.controller_classes.controller_MPI import controller_MPI + + comm = kwargs.get('comm', MPI.COMM_WORLD) + controller = controller_MPI(**controller_args, comm=comm) + P = controller.S.levels[0].prob + else: + controller = controller_nonMPI(**controller_args, num_procs=num_procs) + P = controller.MS[0].levels[0].prob + + uinit = P.u_exact(t0) if u0 is None else u0 + + # insert faults + if fault_stuff is not None: + from pySDC.projects.Resilience.fault_injection import prepare_controller_for_faults + + rnd_args = {'iteration': 1, 'min_node': 1} + args = {'time': 31.0, 'target': 0} + prepare_controller_for_faults(controller, fault_stuff, rnd_args, args) + + # call main function to get things done... + try: + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + except ConvergenceError: + print('Warning: Premature termination!') + stats = controller.return_stats() + return stats, controller, Tend + + +def faults(seed=0): # pragma: no cover + import matplotlib.pyplot as plt + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + + fig, ax = plt.subplots(1, 1) + + rng = np.random.RandomState(seed) + fault_stuff = {'rng': rng, 'args': {}, 'rnd_args': {}} + + controller_params = {'logger_level': 30} + description = {'level_params': {'dt': 1e1}, 'step_params': {'maxiter': 5}} + stats, controller, _ = run_quench(custom_controller_params=controller_params, custom_description=description) + plot_solution_faults(stats, controller, ax, plot_lines=True, label='ref') + + stats, controller, _ = run_quench( + fault_stuff=fault_stuff, + custom_controller_params=controller_params, + ) + plot_solution_faults(stats, controller, ax, label='fixed') + + description['convergence_controllers'] = {Adaptivity: {'e_tol': 1e-7, 'dt_max': 1e2, 'dt_min': 1e-3}} + stats, controller, _ = run_quench( + fault_stuff=fault_stuff, custom_controller_params=controller_params, custom_description=description + ) + + plot_solution_faults(stats, controller, ax, label='adaptivity', ls='--') + plt.show() + + +def plot_solution_faults(stats, controller, ax, plot_lines=False, **kwargs): # pragma: no cover + u_ax = ax + + u = get_sorted(stats, type='u', recomputed=False) + u_ax.plot([me[0] for me in u], [np.mean(me[1]) for me in u], **kwargs) + + if plot_lines: + P = controller.MS[0].levels[0].prob + u_ax.axhline(P.u_thresh, color='grey', ls='-.', label=r'$T_\mathrm{thresh}$') + u_ax.axhline(P.u_max, color='grey', ls=':', label=r'$T_\mathrm{max}$') + + [ax.axvline(me[0], color='grey', label=f'fault at t={me[0]:.2f}') for me in get_sorted(stats, type='bitflip')] + + u_ax.legend() + u_ax.set_xlabel(r'$t$') + u_ax.set_ylabel(r'$T$') + + +def get_crossing_time(stats, controller, num_points=5, inter_points=50, temperature_error_thresh=1e-5): + """ + Compute the time when the temperature threshold is crossed based on interpolation. + + Args: + stats (dict): The stats from a pySDC run + controller (pySDC.Controller.controller): The controller + num_points (int): The number of points in the solution you want to use for interpolation + inter_points (int): The resolution of the interpolation + temperature_error_thresh (float): The temperature error compared to the actual threshold you want to allow + + Returns: + float: The time when the temperature threshold is crossed + """ + from pySDC.core.Lagrange import LagrangeApproximation + from pySDC.core.Collocation import CollBase + + P = controller.MS[0].levels[0].prob + u_thresh = P.u_thresh + + u = get_sorted(stats, type='u', recomputed=False) + temp = np.array([np.mean(me[1]) for me in u]) + t = np.array([me[0] for me in u]) + + crossing_index = np.arange(len(temp))[temp > u_thresh][0] + + # interpolation stuff + num_points = min([num_points, crossing_index * 2, len(temp) - crossing_index]) + idx = np.arange(num_points) - num_points // 2 + crossing_index + t_grid = t[idx] + u_grid = temp[idx] + t_inter = np.linspace(t_grid[0], t_grid[-1], inter_points) + interpolator = LagrangeApproximation(points=t_grid) + u_inter = interpolator.getInterpolationMatrix(t_inter) @ u_grid + + crossing_inter = np.arange(len(u_inter))[u_inter > u_thresh][0] + + temperature_error = abs(u_inter[crossing_inter] - u_thresh) + + assert temperature_error < temp[crossing_index], "Temperature error is rising due to interpolation!" + + if temperature_error > temperature_error_thresh and inter_points < 300: + return get_crossing_time(stats, controller, num_points + 4, inter_points + 15, temperature_error_thresh) + + return t_inter[crossing_inter] + + +def plot_solution(stats, controller): # pragma: no cover + import matplotlib.pyplot as plt + + fig, ax = plt.subplots(1, 1) + u_ax = ax + dt_ax = u_ax.twinx() + + u = get_sorted(stats, type='u', recomputed=False) + u_ax.plot([me[0] for me in u], [np.mean(me[1]) for me in u], label=r'$T$') + + dt = get_sorted(stats, type='dt', recomputed=False) + dt_ax.plot([me[0] for me in dt], [me[1] for me in dt], color='black', ls='--') + u_ax.plot([None], [None], color='black', ls='--', label=r'$\Delta t$') + + if controller.useMPI: + P = controller.S.levels[0].prob + else: + P = controller.MS[0].levels[0].prob + u_ax.axhline(P.u_thresh, color='grey', ls='-.', label=r'$T_\mathrm{thresh}$') + u_ax.axhline(P.u_max, color='grey', ls=':', label=r'$T_\mathrm{max}$') + + [ax.axvline(me[0], color='grey', label=f'fault at t={me[0]:.2f}') for me in get_sorted(stats, type='bitflip')] + + u_ax.legend() + u_ax.set_xlabel(r'$t$') + u_ax.set_ylabel(r'$T$') + dt_ax.set_ylabel(r'$\Delta t$') + + +def compare_imex_full(plotting=False, leak_type='linear'): + """ + Compare the results of IMEX and fully implicit runs. For IMEX we need to limit the step size in order to achieve convergence, but for fully implicit, adaptivity can handle itself better. + + Args: + plotting (bool): Plot the solution or not + """ + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + from pySDC.implementations.hooks.log_work import LogWork + from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun + + maxiter = 5 + num_nodes = 3 + newton_iter_max = 99 + + res = {} + rhs = {} + error = {} + + custom_description = {} + custom_description['problem_params'] = { + 'newton_tol': 1e-10, + 'newton_iter': newton_iter_max, + 'nvars': 2**9, + 'leak_type': leak_type, + } + custom_description['step_params'] = {'maxiter': maxiter} + custom_description['sweeper_params'] = {'num_nodes': num_nodes} + custom_description['convergence_controllers'] = { + Adaptivity: {'e_tol': 1e-6, 'dt_max': 50}, + } + + custom_controller_params = {'logger_level': 30} + for imex in [False, True]: + stats, controller, _ = run_quench( + custom_description=custom_description, + custom_controller_params=custom_controller_params, + imex=imex, + Tend=4.3e2, + use_MPI=False, + hook_class=[LogWork, LogGlobalErrorPostRun], + ) + + res[imex] = get_sorted(stats, type='u')[-1][1] + newton_iter = [me[1] for me in get_sorted(stats, type='work_newton')] + rhs[imex] = np.mean([me[1] for me in get_sorted(stats, type='work_rhs')]) // 1 + error[imex] = get_sorted(stats, type='e_global_post_run')[-1][1] + + if imex: + assert all(me == 0 for me in newton_iter), "IMEX is not supposed to do Newton iterations!" + else: + assert ( + max(newton_iter) / num_nodes / maxiter <= newton_iter_max + ), "Took more Newton iterations than allowed!" + if plotting: # pragma: no cover + plot_solution(stats, controller) + + diff = abs(res[True] - res[False]) + thresh = 4e-3 + assert ( + diff < thresh + ), f"Difference between IMEX and fully-implicit too large! Got {diff:.2e}, allowed is only {thresh:.2e}!" + prob = controller.MS[0].levels[0].prob + assert ( + max(res[True]) > prob.u_max + ), f"Expected runaway to happen, but maximum temperature is {max(res[True]):.2e} < u_max={prob.u_max:.2e}!" + + assert ( + rhs[True] == rhs[False] + ), f"Expected IMEX and fully implicit schemes to take the same number of right hand side evaluations per step, but got {rhs[True]} and {rhs[False]}!" + + assert error[True] < 1e-4, f'Expected error of IMEX version to be less than 1e-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}!' + + +def compare_reference_solutions_single(): + from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostStep, LogLocalErrorPostStep + from pySDC.implementations.hooks.log_solution import LogSolution + + types = ['DIRK', 'SDC', 'scipy'] + types = ['scipy'] + fig, ax = plt.subplots() + error_ax = ax.twinx() + Tend = 500 + + colors = ['black', 'teal', 'magenta'] + + from pySDC.projects.Resilience.strategies import AdaptivityStrategy, merge_descriptions, DoubleAdaptivityStrategy + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + + strategy = DoubleAdaptivityStrategy() + + controller_params = {'logger_level': 15} + + for j in range(len(types)): + description = {} + description['level_params'] = {'dt': 5.0, 'restol': 1e-10} + description['sweeper_params'] = {'QI': 'IE', 'num_nodes': 3} + description['problem_params'] = { + 'leak_type': 'linear', + 'leak_transition': 'step', + 'nvars': 2**10, + 'reference_sol_type': types[j], + 'newton_tol': 1e-12, + } + + description['level_params'] = {'dt': 5.0, 'restol': -1} + description = merge_descriptions(description, strategy.get_custom_description(run_quench, 1)) + description['step_params'] = {'maxiter': 5} + description['convergence_controllers'][Adaptivity]['e_tol'] = 1e-4 + + stats, controller, _ = run_quench( + custom_description=description, + hook_class=[LogGlobalErrorPostStep, LogLocalErrorPostStep, LogSolution], + Tend=Tend, + imex=False, + custom_controller_params=controller_params, + ) + e_glob = get_sorted(stats, type='e_global_post_step', recomputed=False) + e_loc = get_sorted(stats, type='e_local_post_step', recomputed=False) + u = get_sorted(stats, type='u', recomputed=False) + + ax.plot([me[0] for me in u], [max(me[1]) for me in u], color=colors[j], label=f'{types[j]} reference') + + error_ax.plot([me[0] for me in e_glob], [me[1] for me in e_glob], color=colors[j], ls='--') + error_ax.plot([me[0] for me in e_loc], [me[1] for me in e_loc], color=colors[j], ls=':') + + prob = controller.MS[0].levels[0].prob + ax.axhline(prob.u_thresh, ls='-.', color='grey') + ax.axhline(prob.u_max, ls='-.', color='grey') + ax.plot([None], [None], ls='--', label=r'$e_\mathrm{global}$', color='grey') + ax.plot([None], [None], ls=':', label=r'$e_\mathrm{local}$', color='grey') + error_ax.set_yscale('log') + ax.legend(frameon=False) + ax.set_xlabel(r'$t$') + ax.set_ylabel('solution') + error_ax.set_ylabel('error') + ax.set_title('Fully implicit quench problem') + fig.tight_layout() + fig.savefig('data/quench_refs_single.pdf', bbox_inches='tight') + + +def compare_reference_solutions(): + from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun, LogLocalErrorPostStep + + types = ['DIRK', 'SDC', 'scipy'] + 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) + for i in range(len(dt_list)): + description = {} + description['level_params'] = {'dt': dt_list[i], 'restol': 1e-10} + description['sweeper_params'] = {'QI': 'IE', 'num_nodes': 3} + description['problem_params'] = { + 'leak_type': 'linear', + 'leak_transition': 'step', + 'nvars': 2**10, + 'reference_sol_type': types[j], + } + + stats, controller, _ = run_quench( + custom_description=description, + hook_class=[LogGlobalErrorPostRun, LogLocalErrorPostStep], + Tend=Tend, + imex=False, + ) + # errors[i] = get_sorted(stats, type='e_global_post_run')[-1][1] + errors[i] = max([me[1] for me in get_sorted(stats, type='e_local_post_step', recomputed=False)]) + print(errors) + ax.loglog(dt_list, errors, label=f'{types[j]} reference') + + ax.legend(frameon=False) + ax.set_xlabel(r'$\Delta t$') + ax.set_ylabel('global error') + ax.set_title('Fully implicit quench problem') + fig.tight_layout() + fig.savefig('data/quench_refs.pdf', bbox_inches='tight') + + +def check_order(reference_sol_type='scipy'): + from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun + from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import EstimateEmbeddedError + + 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]] + + fig, ax = plt.subplots() + + from pySDC.implementations.sweeper_classes.Runge_Kutta import DIRK34 + + colors = ['black', 'teal', 'magenta', 'orange', 'red'] + for j in range(len(maxiter_list)): + errors = [None] * len(dt_list) + + for i in range(len(dt_list)): + description = {} + description['level_params'] = {'dt': dt_list[i]} + description['step_params'] = {'maxiter': maxiter_list[j]} + description['sweeper_params'] = {'QI': 'IE', 'num_nodes': 3} + description['problem_params'] = { + 'leak_type': 'linear', + 'leak_transition': 'step', + 'nvars': 2**10, + 'reference_sol_type': reference_sol_type, + } + description['convergence_controllers'] = {EstimateEmbeddedError: {}} + + # if maxiter_list[j] == 5: + # description['sweeper_class'] = DIRK34 + # description['sweeper_params'] = {'maxiter': 1} + + stats, controller, _ = run_quench( + custom_description=description, hook_class=[LogGlobalErrorPostRun], Tend=Tend, imex=False + ) + # 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( + dt_list, [errors[0] * (me / dt_list[0]) ** maxiter_list[j] for me in dt_list], color=colors[j], ls='--' + ) + + dt_list = np.array(dt_list) + errors = np.array(errors) + orders = np.log(errors[1:] / errors[:-1]) / np.log(dt_list[1:] / dt_list[:-1]) + print(orders, np.mean(orders)) + + # ax.loglog(dt_list, local_errors) + ax.legend(frameon=False) + ax.set_xlabel(r'$\Delta t$') + ax.set_ylabel('global error') + # ax.set_ylabel('max. local error') + ax.set_title('Fully implicit quench problem') + fig.tight_layout() + fig.savefig(f'data/order_quench_{reference_sol_type}.pdf', bbox_inches='tight') + + +if __name__ == '__main__': + compare_reference_solutions_single() + # for reference_sol_type in ['DIRK', 'SDC', 'scipy']: + # check_order(reference_sol_type=reference_sol_type) + ## faults(19) + ## # get_crossing_time() + # compare_imex_full(plotting=True) + plt.show() diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py new file mode 100644 index 0000000000..be899fb3a0 --- /dev/null +++ b/pySDC/projects/Resilience/strategies.py @@ -0,0 +1,819 @@ +import numpy as np +from matplotlib.colors import TABLEAU_COLORS + +cmap = TABLEAU_COLORS + + +def merge_descriptions(descA, descB): + """ + Merge two dictionaries that may contain dictionaries, which happens when merging descriptions, for instance. + + Keys that occur in both dictionaries will be overwritten by the ones from `descB` and `descA` will be modified, not + copied! + + Args: + descA (dict): Dictionary that you want to merge into + descB (dict): Dictionary you want to merge from + + Returns: + dict: decsA with updated parameters + """ + for key in descB.keys(): + if type(descB[key]) == dict: + descA[key] = merge_descriptions(descA.get(key, {}), descB[key]) + else: + descA[key] = descB[key] + return descA + + +class Strategy: + ''' + Abstract class for resilience strategies + ''' + + def __init__(self, useMPI=False): + ''' + Initialization routine + ''' + self.useMPI = useMPI + + # set default values for plotting + self.linestyle = '-' + self.marker = '.' + self.name = '' + self.bar_plot_x_label = '' + self.color = list(cmap.values())[0] + + # setup custom descriptions + self.custom_description = {} + + # prepare parameters for masks to identify faults that cannot be fixed by this strategy + self.fixable = [] + self.fixable += [ + { + 'key': 'node', + 'op': 'gt', + 'val': 0, + } + ] + self.fixable += [ + { + 'key': 'error', + 'op': 'isfinite', + } + ] + + # stuff for work-precision diagrams + self.precision_parameter = None + self.precision_parameter_loc = [] + + def get_fixable_params(self, **kwargs): + """ + Return a list containing dictionaries which can be passed to `FaultStats.get_mask` as keyword arguments to + obtain a mask of faults that can be fixed + + Returns: + list: Dictionary of parameters + """ + return self.fixable + + def get_fault_args(self, problem, num_procs): + ''' + Routine to get arguments for the faults that are exempt from randomization + + 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: Arguments for the faults that are exempt from randomization + ''' + + return {} + + 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 + ''' + + return {} + + @property + def style(self): + """ + Get the plotting parameters for the strategy. + Supply them to a plotting function using `**` + + Returns: + (dict): The plotting parameters as a dictionary + """ + return { + 'marker': self.marker, + 'label': self.label, + 'color': self.color, + 'ls': self.linestyle, + } + + @property + def label(self): + """ + Get a label for plotting + """ + return self.name + + def get_Tend(self, problem, num_procs=1): + ''' + Get the final time of runs for fault stats based on the problem + + Args: + problem (function): A problem to run + num_procs (int): Number of processes + + Returns: + float: Tend to put into the run + ''' + 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": + return 1.5 + elif problem.__name__ == "run_Schroedinger": + return 1.0 + elif problem.__name__ == "run_quench": + return 500.0 + else: + raise NotImplementedError('I don\'t have a final time for your problem!') + + def get_custom_description(self, problem, num_procs=1): + ''' + Get a custom description based on the problem + + Args: + problem (function): A problem to run + num_procs (int): Number of processes + + Returns: + dict: Custom description + ''' + custom_description = {} + 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, + } + 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} + elif problem.__name__ == "run_Schroedinger": + custom_description['step_params'] = {'maxiter': 5} + custom_description['level_params'] = {'dt': 1e-2, 'restol': -1} + 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_iter': 99, 'newton_tol': 1e-11} + return merge_descriptions(custom_description, self.custom_description) + + +class BaseStrategy(Strategy): + ''' + Do a fixed iteration count + ''' + + def __init__(self, useMPI=False): + ''' + Initialization routine + ''' + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[0] + self.marker = 'o' + self.name = 'base' + self.bar_plot_x_label = 'base' + self.precision_parameter = 'dt' + self.precision_parameter_loc = ['level_params', 'dt'] + + @property + def label(self): + return r'fixed' + + +class AdaptivityStrategy(Strategy): + ''' + Adaptivity as a resilience strategy + ''' + + def __init__(self, useMPI=False): + ''' + Initialization routine + ''' + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[1] + self.marker = '*' + self.name = 'adaptivity' + self.bar_plot_x_label = 'adaptivity' + self.precision_parameter = 'e_tol' + self.precision_parameter_loc = ['convergence_controllers', Adaptivity, 'e_tol'] + + def get_fixable_params(self, maxiter, **kwargs): + """ + Here faults occurring in the last iteration cannot be fixed. + + Args: + maxiter (int): Max. iterations until convergence is declared + + Returns: + (list): Contains dictionaries of keyword arguments for `FaultStats.get_mask` + """ + self.fixable += [ + { + 'key': 'iteration', + 'op': 'lt', + 'val': maxiter, + } + ] + return self.fixable + + 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 Adaptivity + + 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. + + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestarting + + flavor = 'MPI' if self.useMPI else 'nonMPI' + custom_description['convergence_controllers'][BasicRestarting.get_implementation(flavor)] = { + 'max_restarts': 15 + } + else: + raise NotImplementedError( + 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\ + strategy' + ) + + custom_description['convergence_controllers'][Adaptivity] = { + 'e_tol': e_tol, + 'dt_min': dt_min, + 'dt_max': dt_max, + 'dt_slope_max': dt_slope_max, + } + return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + + +class AdaptiveHotRodStrategy(Strategy): + ''' + Adaptivity + Hot Rod as a resilience strategy + ''' + + def __init__(self, useMPI=False): + ''' + Initialization routine + ''' + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[4] + self.marker = '.' + self.name = 'adaptive Hot Rod' + self.bar_plot_x_label = 'adaptive\nHot Rod' + self.precision_parameter = 'e_tol' + self.precision_parameter_loc = ['convergence_controllers', Adaptivity, 'e_tol'] + + def get_custom_description(self, problem, num_procs): + ''' + Routine to get a custom description that adds adaptivity and Hot Rod + + 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 description you can supply to the problem when running it + ''' + from pySDC.implementations.convergence_controller_classes.hotrod import HotRod + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + + if problem.__name__ == "run_vdp": + e_tol = 3e-7 + dt_min = 1e-3 + maxiter = 4 + HotRod_tol = 3e-7 + else: + raise NotImplementedError( + 'I don\'t have a tolerance for adaptive Hot Rod for your problem. Please add one \ +to the strategy' + ) + + no_storage = num_procs > 1 + + custom_description = { + 'convergence_controllers': { + Adaptivity: {'e_tol': e_tol, 'dt_min': dt_min}, + HotRod: {'HotRod_tol': HotRod_tol, 'no_storage': no_storage}, + }, + 'step_params': {'maxiter': maxiter}, + } + + return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + + +class IterateStrategy(Strategy): + ''' + Iterate for as much as you want + ''' + + def __init__(self, useMPI=False): + ''' + Initialization routine + ''' + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[2] + self.marker = 'v' + self.name = 'iterate' + self.bar_plot_x_label = 'iterate' + self.precision_parameter = 'restol' + self.precision_parameter_loc = ['level_params', 'restol'] + + @property + def label(self): + return r'$k$ adaptivity' + + def get_custom_description(self, problem, num_procs): + ''' + Routine to get a custom description that allows for adaptive iteration counts + + 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 description you can supply to the problem when running it + ''' + restol = -1 + e_tol = -1 + + if problem.__name__ == "run_piline": + restol = 2.3e-8 + elif problem.__name__ == "run_vdp": + restol = 9e-7 + elif problem.__name__ == "run_Lorenz": + restol = 16e-7 + elif problem.__name__ == "run_Schroedinger": + restol = 6.5e-7 + elif problem.__name__ == "run_quench": + restol = 1e-7 + else: + raise NotImplementedError( + 'I don\'t have a residual tolerance for your problem. Please add one to the \ +strategy' + ) + + custom_description = { + 'step_params': {'maxiter': 99}, + 'level_params': {'restol': restol, 'e_tol': e_tol}, + } + + if problem.__name__ == "run_quench": + custom_description['level_params']['dt'] = 1 + + return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + + +class HotRodStrategy(Strategy): + ''' + Hot Rod as a resilience strategy + ''' + + def __init__(self, useMPI=False): + ''' + Initialization routine + ''' + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[3] + self.marker = '^' + self.name = 'Hot Rod' + self.bar_plot_x_label = 'Hot Rod' + self.precision_parameter = 'dt' + self.precision_parameter_loc = ['level_params', 'dt'] + + def get_custom_description(self, problem, num_procs): + ''' + Routine to get a custom description that adds Hot Rod + + 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 description you can supply to the problem when running it + ''' + from pySDC.implementations.convergence_controller_classes.hotrod import HotRod + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI + + if problem.__name__ == "run_vdp": + HotRod_tol = 5e-7 + maxiter = 4 + elif problem.__name__ == "run_Lorenz": + HotRod_tol = 4e-7 + maxiter = 6 + elif problem.__name__ == "run_Schroedinger": + HotRod_tol = 3e-7 + maxiter = 6 + elif problem.__name__ == "run_quench": + HotRod_tol = 3e-5 + maxiter = 6 + else: + raise NotImplementedError( + 'I don\'t have a tolerance for Hot Rod for your problem. Please add one to the\ + strategy' + ) + + no_storage = num_procs > 1 + + custom_description = { + 'convergence_controllers': { + HotRod: {'HotRod_tol': HotRod_tol, 'no_storage': no_storage}, + BasicRestartingNonMPI: {'max_restarts': 2, 'crash_after_max_restarts': False}, + }, + 'step_params': {'maxiter': maxiter}, + } + + return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + + +class AdaptivityCollocationStrategy(Strategy): + ''' + Adaptivity based on collocation as a resilience strategy + ''' + + def __init__(self, useMPI=False): + ''' + Initialization routine + ''' + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityCollocation + + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[1] + self.marker = '*' + self.name = 'adaptivity_coll' + self.bar_plot_x_label = 'adaptivity collocation' + 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): + ''' + 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 AdaptivityCollocation + + custom_description = {} + custom_description['step_params'] = {'maxiter': self.maxiter} + + dt_max = np.inf + dt_min = 1e-5 + + 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 = 1e2 + 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, + } + } + return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + + +class AdaptivityCollocationTypeStrategy(AdaptivityCollocationStrategy): + def __init__(self, useMPI=False): + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[4] + self.marker = '.' + self.adaptive_coll_params = { + 'quad_type': ['RADAU-RIGHT', 'GAUSS'], + 'do_coll_update': [False, True], + } + # self.adaptive_coll_params = {'quad_type': ['RADAU-LEFT', 'RADAU-RIGHT'], 'num_nodes': [3, 3]} + + @property + def label(self): + return 'adaptivity type' + + +class AdaptivityCollocationRefinementStrategy(AdaptivityCollocationStrategy): + def __init__(self, useMPI=False): + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[5] + self.marker = '^' + self.adaptive_coll_params = { + 'num_nodes': [2, 3], + 'quad_type': ['GAUSS', 'RADAU-RIGHT'], + 'do_coll_update': [True, False], + } + + @property + def label(self): + return 'adaptivity refinement' + + +class AdaptivityCollocationDerefinementStrategy(AdaptivityCollocationStrategy): + def __init__(self, useMPI=False): + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[6] + self.marker = '^' + self.adaptive_coll_params = {'num_nodes': [4, 3]} + + @property + def label(self): + return 'adaptivity de-refinement' + + +class DIRKStrategy(AdaptivityStrategy): + ''' + DIRK4(3) + ''' + + def __init__(self, useMPI=False): + ''' + Initialization routine + ''' + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK + + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[7] + self.marker = '^' + self.name = 'DIRK' + self.bar_plot_x_label = 'DIRK4(3)' + self.precision_parameter = 'e_tol' + self.precision_parameter_loc = ['convergence_controllers', AdaptivityRK, 'e_tol'] + + @property + def label(self): + return 'DIRK4(3)' + + 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 DIRK34 + + adaptivity_description = super().get_custom_description(problem, num_procs) + + e_tol = adaptivity_description['convergence_controllers'][Adaptivity]['e_tol'] + adaptivity_description['convergence_controllers'].pop(Adaptivity, None) + adaptivity_description.pop('sweeper_params', None) + + rk_params = { + 'step_params': {'maxiter': 1}, + 'sweeper_class': DIRK34, + 'convergence_controllers': {AdaptivityRK: {'e_tol': e_tol}}, + } + + custom_description = merge_descriptions(adaptivity_description, rk_params) + + return custom_description + + +class ERKStrategy(DIRKStrategy): + """ + Explicit embedded RK using Cash-Karp's method + """ + + def __init__(self, useMPI=False): + ''' + Initialization routine + ''' + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[9] + self.marker = 'x' + self.name = 'ERK' + self.bar_plot_x_label = 'ERK5(4)' + + @property + def label(self): + return 'CP5(4)' + + """ + Explicit Cash-Karp's method + """ + + def get_custom_description(self, problem, num_procs=1): + from pySDC.implementations.sweeper_classes.Runge_Kutta import Cash_Karp + + desc = super().get_custom_description(problem, num_procs) + desc['sweeper_class'] = Cash_Karp + return desc + + +class DoubleAdaptivityStrategy(AdaptivityStrategy): + ''' + Adaptivity based both on embedded estimate and on residual + ''' + + def __init__(self, useMPI=False): + ''' + Initialization routine + ''' + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + + super().__init__(useMPI=useMPI) + self.color = list(cmap.values())[7] + self.marker = '^' + self.name = 'double_adaptivity' + self.bar_plot_x_label = 'double adaptivity' + self.precision_parameter = 'e_tol' + self.precision_parameter_loc = ['convergence_controllers', Adaptivity, 'e_tol'] + self.residual_e_tol_ratio = 1.0 + self.residual_e_tol_abs = None + + @property + def label(self): + return 'double adaptivity' + + 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 AdaptivityResidual, Adaptivity + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestarting + + custom_description = super().get_custom_description(problem, num_procs) + + if self.residual_e_tol_abs: + e_tol = self.residual_e_tol_abs + else: + e_tol = custom_description['convergence_controllers'][Adaptivity]['e_tol'] * self.residual_e_tol_ratio + custom_description['convergence_controllers'][AdaptivityResidual] = { + 'e_tol': e_tol, + 'allowed_modifications': ['decrease'], + } + + flavor = 'MPI' if self.useMPI else 'nonMPI' + custom_description['convergence_controllers'][BasicRestarting.get_implementation(flavor)] = {'max_restarts': 15} + + return custom_description + + +class AdaptivityAvoidRestartsStrategy(AdaptivityStrategy): + """ + Adaptivity with the avoid restarts option + """ + + @property + def label(self): + return 'adaptivity (avoid restarts)' + + 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 Adaptivity + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestarting + + custom_description = super().get_custom_description(problem, num_procs) + + custom_description['convergence_controllers'][Adaptivity]['avoid_restarts'] = True + + flavor = 'MPI' if self.useMPI else 'nonMPI' + custom_description['convergence_controllers'][BasicRestarting.get_implementation(flavor)] = {'max_restarts': 15} + + return custom_description + + +class AdaptivityInterpolationStrategy(AdaptivityStrategy): + """ + Adaptivity with interpolation between restarts + """ + + @property + def label(self): + return 'adaptivity+interpolation' + + 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 Adaptivity + from pySDC.implementations.convergence_controller_classes.interpolate_between_restarts import ( + InterpolateBetweenRestarts, + ) + from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestarting + + custom_description = super().get_custom_description(problem, num_procs) + + custom_description['convergence_controllers'][Adaptivity]['avoid_restarts'] = False + custom_description['convergence_controllers'][InterpolateBetweenRestarts] = {} + + flavor = 'MPI' if self.useMPI else 'nonMPI' + custom_description['convergence_controllers'][BasicRestarting.get_implementation(flavor)] = {'max_restarts': 15} + + return custom_description diff --git a/pySDC/projects/Resilience/vdp.py b/pySDC/projects/Resilience/vdp.py index ce7dd5a132..2c7177e03d 100644 --- a/pySDC/projects/Resilience/vdp.py +++ b/pySDC/projects/Resilience/vdp.py @@ -9,6 +9,7 @@ from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity from pySDC.core.Errors import ProblemError, ConvergenceError from pySDC.projects.Resilience.hook import LogData, hook_collection +from pySDC.projects.Resilience.strategies import merge_descriptions def plot_step_sizes(stats, ax, e_em_key='error_embedded_estimate'): @@ -88,7 +89,6 @@ def run_vdp( hook_class=LogData, fault_stuff=None, custom_controller_params=None, - custom_problem_params=None, use_MPI=False, **kwargs, ): @@ -102,7 +102,6 @@ def run_vdp( hook_class (pySDC.Hook): A hook to store data fault_stuff (dict): A dictionary with information on how to add faults custom_controller_params (dict): Overwrite presets - custom_problem_params (dict): Overwrite presets use_MPI (bool): Whether or not to use MPI Returns: @@ -128,9 +127,6 @@ def run_vdp( 'u0': np.array([2.0, 0.0]), } - if custom_problem_params is not None: - problem_params = {**problem_params, **custom_problem_params} - # initialize step parameters step_params = {} step_params['maxiter'] = 4 @@ -154,11 +150,7 @@ def run_vdp( description['step_params'] = step_params if custom_description is not None: - for k in custom_description.keys(): - if k == 'sweeper_class': - description[k] = custom_description[k] - continue - description[k] = {**description.get(k, {}), **custom_description.get(k, {})} + description = merge_descriptions(description, custom_description) # set time parameters t0 = 0.0 @@ -188,7 +180,8 @@ def run_vdp( from pySDC.projects.Resilience.fault_injection import prepare_controller_for_faults rnd_args = {'iteration': 3} - args = {'time': 1.0, 'target': 0} + # args = {'time': 0.9, 'target': 0} + args = {'time': 5.25, 'target': 0} prepare_controller_for_faults(controller, fault_stuff, rnd_args, args) # call main function to get things done... @@ -442,6 +435,12 @@ def check_step_size_limiter(size=4, comm=None): assert ( dt_min >= expect['dt_min'] ), f"Exceeded minimum allowed step size! Got {dt_min:.4e}, allowed {params['dt_min']:.4e}." + assert ( + dt_slope_max <= expect['dt_slope_max'] + ), f"Exceeded maximum allowed step size slope! Got {dt_slope_max:.4e}, allowed {params['dt_slope_max']:.4e}." + assert ( + dt_slope_min >= expect['dt_slope_min'] + ), f"Exceeded minimum allowed step size slope! Got {dt_slope_min:.4e}, allowed {params['dt_slope_min']:.4e}." assert ( dt_slope_max <= expect['dt_slope_max'] diff --git a/pySDC/projects/Resilience/work_precision.py b/pySDC/projects/Resilience/work_precision.py new file mode 100644 index 0000000000..05a9151049 --- /dev/null +++ b/pySDC/projects/Resilience/work_precision.py @@ -0,0 +1,1000 @@ +from mpi4py import MPI +import numpy as np +import matplotlib.pyplot as plt +import pickle + +from pySDC.projects.Resilience.strategies import merge_descriptions +from pySDC.projects.Resilience.Lorenz import run_Lorenz +from pySDC.projects.Resilience.vdp import run_vdp +from pySDC.projects.Resilience.Schroedinger import run_Schroedinger +from pySDC.projects.Resilience.quench import run_quench + +from pySDC.helpers.stats_helper import get_sorted +from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal + +setup_mpl(reset=True) +LOGGER_LEVEL = 30 +VERBOSE = True + +MAPPINGS = { + 'e_global': ('e_global_post_run', max, False), + 'e_global_rel': ('e_global_rel_post_run', max, False), + 't': ('timing_run', max, False), + # 'e_local_max': ('e_local_post_step', max, False), + 'k_SDC': ('k', sum, None), + 'k_SDC_no_restart': ('k', sum, False), + 'k_Newton': ('work_newton', 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), + 'e_embedded_max': ('error_embedded_estimate', max, False), +} + + +def single_run(problem, strategy, data, custom_description, num_procs=1, comm_world=None, problem_args=None): + """ + Make a single run of a particular problem with a certain strategy. + + Args: + problem (function): A problem to run + strategy (Strategy): SDC strategy + data (dict): Put the results in here + custom_description (dict): Overwrite presets + num_procs (int): Number of processes for the time communicator + comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script + + Returns: + None + """ + from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRunMPI + from pySDC.implementations.hooks.log_work import LogWork + from pySDC.projects.Resilience.hook import LogData + + comm = comm_world.Split(comm_world.rank < num_procs) + if comm_world.rank >= num_procs: + return None + + strategy_description = strategy.get_custom_description(problem, num_procs) + description = merge_descriptions(strategy_description, custom_description) + + controller_params = {'logger_level': LOGGER_LEVEL} + problem_args = {} if problem_args is None else problem_args + + stats, controller, _ = problem( + custom_description=description, + Tend=strategy.get_Tend(problem, num_procs), + hook_class=[LogData, LogWork, LogGlobalErrorPostRunMPI], + custom_controller_params=controller_params, + use_MPI=True, + comm=comm, + **problem_args, + ) + + # record all the metrics + for key, mapping in MAPPINGS.items(): + me = get_sorted(stats, type=mapping[0], recomputed=mapping[2], comm=comm) + if len(me) == 0: + data[key] += [np.nan] + else: + data[key] += [mapping[1]([you[1] for you in me])] + return None + + +def get_parameter(dictionary, where): + """ + Get a parameter at a certain position in a dictionary of dictionaries. + + Args: + dictionary (dict): The dictionary + where (list): The list of keys leading to the value you want + + Returns: + The value of the dictionary + """ + if len(where) == 1: + return dictionary[where[0]] + else: + return get_parameter(dictionary[where[0]], where[1:]) + + +def set_parameter(dictionary, where, parameter): + """ + Set a parameter at a certain position in a dictionary of dictionaries + + Args: + dictionary (dict): The dictionary + where (list): The list of keys leading to the value you want to set + parameter: Whatever you want to set the parameter to + + Returns: + None + """ + if len(where) == 1: + dictionary[where[0]] = parameter + else: + set_parameter(dictionary[where[0]], where[1:], parameter) + + +def get_path(problem, strategy, num_procs, handle='', base_path='data/work_precision'): + """ + Get the path to a certain data. + + Args: + problem (function): A problem to run + strategy (Strategy): SDC strategy + 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 + + 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' + + +def record_work_precision( + problem, + strategy, + num_procs=1, + custom_description=None, + handle='', + runs=1, + comm_world=None, + problem_args=None, + param_range=None, +): + """ + Run problem with strategy and record the cost parameters. + + Args: + problem (function): A problem to run + strategy (Strategy): SDC strategy + num_procs (int): Number of processes for the time communicator + custom_description (dict): Overwrite presets + 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 + + Returns: + None + """ + from pySDC.implementations.convergence_controller_classes.estimate_embedded_error import EstimateEmbeddedError + + data = {} + + estimate_embedded_error = {'convergence_controllers': {EstimateEmbeddedError: {}}} + + # prepare precision parameters + param = strategy.precision_parameter + description = merge_descriptions( + merge_descriptions( + strategy.get_custom_description(problem, num_procs), + {} if custom_description is None else custom_description, + ), + estimate_embedded_error, + ) + 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] + if problem.__name__ == 'run_vdp': + exponents = [-4, -3, -2, -1, 0, 1, 2] + elif param == 'dt': + power = 2.0 + exponents = [-1, 0, 1, 2, 3] + elif param == 'restol': + power = 10.0 + exponents = [-2, -1, 0, 1, 2, 3] + if problem.__name__ == 'run_vdp': + exponents = [-4, -3, -2, -1, 0, 1] + else: + raise NotImplementedError(f"I don't know how to get default value for parameter \"{param}\"") + + where = strategy.precision_parameter_loc + default = get_parameter(description, where) + param_range = [default * power**i for i in exponents] if param_range is None else param_range + + 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]] + + # 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'], 1e-9) + set_parameter(description, ['level_params', 'restol'], param_range[i] / 10.0) + + 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]] + for _j in range(runs): + single_run( + problem, + strategy, + data[param_range[i]], + custom_description=description, + comm_world=comm_world, + problem_args=problem_args, + num_procs=num_procs, + ) + + comm_world.Barrier() + + if VERBOSE and comm_world.rank == 0: + print( + f'{problem.__name__} {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]}' + ) + + if comm_world.rank == 0: + import socket + import time + + data['meta'] = { + 'hostname': socket.gethostname(), + 'time': time.time, + 'runs': runs, + } + with open(get_path(problem, strategy, num_procs, handle), 'wb') as f: + pickle.dump(data, f) + + +def plot_work_precision( + problem, + strategy, + num_procs, + ax, + work_key='k_SDC', + precision_key='e_global', + handle='', + plotting_params=None, + comm_world=None, +): # pragma: no cover + """ + Plot data from running a problem with a strategy. + + Args: + problem (function): A problem to run + strategy (Strategy): SDC strategy + num_procs (int): Number of processes for the time communicator + ax (matplotlib.pyplot.axes): Somewhere to plot + work_key (str): The key in the recorded data you want on the x-axis + precision_key (str): The key in the recorded data you want on the y-axis + 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 + + Returns: + None + """ + if comm_world.rank > 0: + return None + + with open(get_path(problem, strategy, num_procs, handle=handle), 'rb') as f: + data = pickle.load(f) + + keys = [key for key in data.keys() if key not in ['meta']] + work = [np.nanmean(data[key][work_key]) for key in keys] + precision = [np.nanmean(data[key][precision_key]) for key in keys] + + for key in [work_key, precision_key]: + 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): + print( + f"WARNING: Variance in \"{key}\" for {get_path(problem, strategy, num_procs, handle)} too large! Got {rel_variance}" + ) + + style = merge_descriptions( + {**strategy.style, 'label': f'{strategy.style["label"]}{f" {handle}" if handle else ""}'}, + plotting_params if plotting_params else {}, + ) + + ax.loglog(work, precision, **style) + + if 't' in [work_key, precision_key]: + meta = 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) + + +def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only=False): # pragma: no cover + """ + Decorate a plot + + Args: + ax (matplotlib.pyplot.axes): Somewhere to plot + problem (function): A problem to run + work_key (str): The key in the recorded data you want on the x-axis + precision_key (str): The key in the recorded data you want on the y-axis + num_procs (int): Number of processes for the time communicator + title_only (bool): Put only the title on top, or do the whole shebang + + Returns: + None + """ + labels = { + 'k_SDC': 'SDC iterations', + 'k_SDC_no_restart': 'SDC iterations (restarts excluded)', + 'k_Newton': 'Newton iterations', + 'k_Newton_no_restart': 'Newton iterations (restarts excluded)', + 'k_rhs': 'right hand side evaluations', + 't': 'wall clock time / s', + 'e_global': 'global error', + 'e_global_rel': 'relative global error', + 'e_local_max': 'max. local error', + 'restart': 'restarts', + 'dt_max': r'$\Delta t_\mathrm{max}$', + 'dt_mean': r'$\bar{\Delta t}$', + 'param': 'parameter', + } + + if not title_only: + ax.set_xlabel(labels.get(work_key, 'work')) + ax.set_ylabel(labels.get(precision_key, 'precision')) + # ax.legend(frameon=False) + + titles = { + 'run_vdp': 'Van der Pol', + 'run_Lorenz': 'Lorenz attractor', + 'run_Schroedinger': r'Schr\"odinger', + 'run_quench': 'Quench', + } + ax.set_title(titles.get(problem.__name__, '')) + + +def execute_configurations( + problem, + configurations, + work_key, + precision_key, + num_procs, + ax, + decorate, + record, + runs, + comm_world, + plotting, +): + """ + Run for multiple configurations. + + Args: + problem (function): A problem to run + configurations (dict): The configurations you want to run with + work_key (str): The key in the recorded data you want on the x-axis + precision_key (str): The key in the recorded data you want on the y-axis + num_procs (int): Number of processes for the time communicator + ax (matplotlib.pyplot.axes): Somewhere to plot + decorate (bool): Whether to decorate fully or only put the title + record (bool): Whether to only plot or also record the data first + 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 + + Returns: + None + """ + for _, config in configurations.items(): + for strategy in config['strategies']: + shared_args = { + 'problem': problem, + 'strategy': strategy, + 'handle': config.get('handle', ''), + 'num_procs': config.get('num_procs', num_procs), + } + if record: + record_work_precision( + **shared_args, + custom_description=config.get('custom_description', {}), + runs=runs, + comm_world=comm_world, + problem_args=config.get('problem_args', {}), + param_range=config.get('param_range', None), + ) + if plotting and comm_world.rank == 0: + plot_work_precision( + **shared_args, + work_key=work_key, + precision_key=precision_key, + ax=ax, + plotting_params=config.get('plotting_params', {}), + comm_world=comm_world, + ) + + decorate_panel( + ax=ax, + problem=problem, + work_key=work_key, + precision_key=precision_key, + num_procs=num_procs, + title_only=not decorate, + ) + + +def get_configs(mode, problem): + """ + Get configurations for work-precision plots. These are dictionaries containing strategies and handles and so on. + + Args: + mode (str): The of the configurations you want to retrieve + problem (function): A problem to run + + Returns: + dict: Configurations + """ + configurations = {} + if mode == 'regular': + from pySDC.projects.Resilience.strategies import AdaptivityStrategy, BaseStrategy, IterateStrategy + + handle = 'regular' + configurations[0] = { + 'handle': handle, + 'strategies': [AdaptivityStrategy(useMPI=True), BaseStrategy(useMPI=True), IterateStrategy(useMPI=True)], + } + elif mode == 'step_size_limiting': + from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter + from pySDC.projects.Resilience.strategies import AdaptivityStrategy + + configurations[0] = { + 'custom_description': {'convergence_controllers': {StepSizeLimiter: {'dt_max': 25}}}, + 'handle': 'step limiter', + 'strategies': [AdaptivityStrategy(useMPI=True)], + 'plotting_params': {'color': 'teal', 'marker': 'v'}, + } + configurations[1] = { + 'custom_description': {'convergence_controllers': {StepSizeLimiter: {'dt_slope_max': 2}}}, + 'handle': 'slope limiter', + 'strategies': [AdaptivityStrategy(useMPI=True)], + 'plotting_params': {'color': 'magenta', 'marker': 'x'}, + } + configurations[2] = { + 'custom_description': {}, + 'handle': 'no limits', + 'plotting_params': {'label': 'adaptivity'}, + 'strategies': [AdaptivityStrategy(useMPI=True)], + } + 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': '--'} + + 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, + } + + 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, + } + 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)], + } + 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)'}, + } + elif mode == 'parallel_efficiency': + from pySDC.projects.Resilience.strategies import AdaptivityStrategy, BaseStrategy, IterateStrategy, ERKStrategy + + desc = {} + desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'} + desc['step_params'] = {'maxiter': 5} + + descIterate = {} + descIterate['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'} + + ls = { + 1: '-', + 2: '--', + 3: '-.', + 4: ':', + 5: 'loosely dashdotted', + } + + # configurations[-1] = { + # 'strategies': [ERKStrategy(useMPI=False)], 'num_procs':1, + # } + + for num_procs in [4, 2, 1]: + plotting_params = {'ls': ls[num_procs], 'label': f'adaptivity {num_procs} procs'} + configurations[num_procs] = { + 'strategies': [AdaptivityStrategy(True)], + 'custom_description': desc, + 'num_procs': num_procs, + 'plotting_params': plotting_params, + } + plotting_params = {'ls': ls[num_procs], 'label': fr'$k$ adaptivity {num_procs} procs'} + configurations[num_procs + 100] = { + 'strategies': [IterateStrategy(True)], + 'custom_description': descIterate, + 'num_procs': num_procs, + 'plotting_params': plotting_params, + } + + elif mode[:13] == 'vdp_stiffness': + from pySDC.projects.Resilience.strategies import AdaptivityStrategy, ERKStrategy, DIRKStrategy + + mu = float(mode[14:]) + + problem_desc = {'problem_params': {'mu': mu}} + + desc = {} + desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'} + desc['step_params'] = {'maxiter': 5} + desc['problem_params'] = problem_desc['problem_params'] + + ls = { + 1: '-', + 2: '--', + 3: '-.', + 4: ':', + 5: 'loosely dashdotted', + } + + for num_procs in [4, 1]: + plotting_params = {'ls': ls[num_procs], 'label': f'SDC {num_procs} procs'} + configurations[num_procs] = { + 'strategies': [AdaptivityStrategy(True)], + 'custom_description': desc, + 'num_procs': num_procs, + 'plotting_params': plotting_params, + 'handle': mode, + } + + configurations[2] = { + 'strategies': [ERKStrategy(useMPI=True)], + 'num_procs': 1, + '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, + } + + elif mode == 'compare_adaptivity': + # TODO: configurations not final! + from pySDC.projects.Resilience.strategies import ( + AdaptivityCollocationTypeStrategy, + AdaptivityCollocationRefinementStrategy, + AdaptivityStrategy, + ) + + strategies = [ + AdaptivityCollocationTypeStrategy(useMPI=True), + AdaptivityCollocationRefinementStrategy(useMPI=True), + ] + + restol = None + for strategy in strategies: + strategy.restol = restol + + configurations[1] = {'strategies': strategies} + configurations[2] = { + 'custom_description': {'step_params': {'maxiter': 5}}, + 'strategies': [AdaptivityStrategy(useMPI=True)], + } + + # 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': + from pySDC.projects.Resilience.strategies import ( + AdaptivityStrategy, + DoubleAdaptivityStrategy, + IterateStrategy, + BaseStrategy, + ) + + dumbledoresarmy = DoubleAdaptivityStrategy(useMPI=True) + # dumbledoresarmy.residual_e_tol_ratio = 1e2 + dumbledoresarmy.residual_e_tol_abs = 1e-3 + + strategies = [ + AdaptivityStrategy(useMPI=True), + IterateStrategy(useMPI=True), + BaseStrategy(useMPI=True), + dumbledoresarmy, + ] + configurations[1] = {'strategies': strategies} + configurations[2] = { + 'strategies': strategies, + 'problem_args': {'imex': True}, + 'handle': 'IMEX', + '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': '-.'}, + } + elif mode == 'preconditioners': + from pySDC.projects.Resilience.strategies import AdaptivityStrategy, IterateStrategy, BaseStrategy + + strategies = [AdaptivityStrategy(useMPI=True), IterateStrategy(useMPI=True), BaseStrategy(useMPI=True)] + + precons = ['IE', 'LU', 'MIN'] + ls = ['-', '--', '-.', ':'] + for i in range(len(precons)): + configurations[i] = { + 'strategies': strategies, + 'custom_description': {'sweeper_params': {'QI': precons[i]}}, + 'handle': precons[i], + 'plotting_params': {'ls': ls[i]}, + } + + 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]}, + } + elif mode == 'avoid_restarts': + from pySDC.projects.Resilience.strategies import ( + AdaptivityStrategy, + AdaptivityAvoidRestartsStrategy, + AdaptivityInterpolationStrategy, + ) + + desc = {'sweeper_params': {'QI': 'IE'}, 'step_params': {'maxiter': 3}} + param_range = [1e-3, 1e-5] + configurations[0] = { + 'strategies': [AdaptivityInterpolationStrategy(useMPI=True)], + 'plotting_params': {'ls': '--'}, + 'custom_description': desc, + 'param_range': param_range, + } + configurations[1] = { + 'strategies': [AdaptivityAvoidRestartsStrategy(useMPI=True)], + 'plotting_params': {'ls': '-.'}, + 'custom_description': desc, + 'param_range': param_range, + } + configurations[2] = { + 'strategies': [AdaptivityStrategy(useMPI=True)], + 'custom_description': desc, + 'param_range': param_range, + } + else: + raise NotImplementedError(f'Don\'t know the mode "{mode}"!') + + return configurations + + +def get_fig(x=1, y=1, **kwargs): # pragma: no cover + """ + Get a figure to plot in. + + Args: + x (int): How many panels in horizontal direction you want + y (int): How many panels in vertical direction you want + + Returns: + matplotlib.pyplot.Figure + """ + width = 1.0 + ratio = 1.0 if y == 2 else 0.5 + keyword_arguments = { + 'figsize': figsize_by_journal('Springer_Numerical_Algorithms', width, ratio), + 'layout': 'constrained', + **kwargs, + } + return plt.subplots(y, x, **keyword_arguments) + + +def save_fig( + fig, name, work_key, precision_key, legend=True, format='pdf', base_path='data', **kwargs +): # pragma: no cover + """ + Save a figure with a legend on the bottom. + + Args: + fig (matplotlib.pyplot.Figure): Figure you want to save + name (str): Name of the plot to put in the path + work_key (str): The key in the recorded data you want on the x-axis + 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 + + Returns: + None + """ + handles, labels = fig.get_axes()[0].get_legend_handles_labels() + 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, + frameon=False, + fancybox=True, + ) + + path = f'{base_path}/wp-{name}-{work_key}-{precision_key}.{format}' + fig.savefig(path, bbox_inches='tight', **kwargs) + print(f'Stored figure \"{path}\"') + + +def all_problems(mode='compare_strategies', plotting=True, base_path='data', **kwargs): # pragma: no cover + """ + Make a plot comparing various strategies for all problems. + + Args: + work_key (str): The key in the recorded data you want on the x-axis + precision_key (str): The key in the recorded data you want on the y-axis + + Returns: + None + """ + + fig, axs = get_fig(2, 2) + + shared_params = { + 'work_key': 'k_SDC', + 'precision_key': 'e_global', + 'num_procs': 1, + 'runs': 1, + 'comm_world': MPI.COMM_WORLD, + 'record': False, + 'plotting': plotting, + **kwargs, + } + + problems = [run_vdp, run_Lorenz, run_Schroedinger, run_quench] + + for i in range(len(problems)): + execute_configurations( + **shared_params, + problem=problems[i], + ax=axs.flatten()[i], + decorate=True, + configurations=get_configs(mode, problems[i]), + ) + + if plotting and shared_params['comm_world'].rank == 0: + save_fig( + fig=fig, + name=mode, + work_key=shared_params['work_key'], + precision_key=shared_params['precision_key'], + legend=True, + base_path=base_path, + ) + + +def ODEs(mode='compare_strategies', plotting=True, base_path='data', **kwargs): # pragma: no cover + """ + Make a plot comparing various strategies for the two ODEs. + + Args: + work_key (str): The key in the recorded data you want on the x-axis + precision_key (str): The key in the recorded data you want on the y-axis + + Returns: + None + """ + + fig, axs = get_fig(x=2, y=1) + + shared_params = { + 'work_key': 'k_SDC', + 'precision_key': 'e_global', + 'num_procs': 1, + 'runs': 1, + 'comm_world': MPI.COMM_WORLD, + 'record': False, + 'plotting': plotting, + **kwargs, + } + + problems = [run_vdp, run_Lorenz] + + for i in range(len(problems)): + execute_configurations( + **shared_params, + problem=problems[i], + ax=axs.flatten()[i], + decorate=i == 0, + configurations=get_configs(mode, problems[i]), + ) + + if plotting and shared_params['comm_world'].rank == 0: + save_fig( + fig=fig, + name=f'ODEs-{mode}', + work_key=shared_params['work_key'], + precision_key=shared_params['precision_key'], + legend=True, + base_path=base_path, + ) + + +def single_problem(mode, problem, plotting=True, base_path='data', **kwargs): # pragma: no cover + """ + Make a plot for a single problem + + Args: + 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)) + + params = { + 'work_key': 'k_SDC', + 'precision_key': 'e_global', + 'num_procs': 1, + 'runs': 1, + 'comm_world': MPI.COMM_WORLD, + 'record': False, + 'plotting': plotting, + **kwargs, + } + + execute_configurations(**params, problem=problem, ax=ax, decorate=True, configurations=get_configs(mode, problem)) + + if plotting: + save_fig( + fig=fig, + name=f'{problem.__name__}-{mode}', + work_key=params['work_key'], + precision_key=params['precision_key'], + legend=False, + base_path=base_path, + ) + + +def vdp_stiffness_plot(base_path='data', format='pdf', **kwargs): # pragma: no cover + fig, axs = get_fig(2, 2, sharex=True) + + mus = [0, 5, 10, 15] + + for i in range(len(mus)): + params = { + 'runs': 1, + 'problem': run_vdp, + 'record': False, + 'work_key': 't', + 'precision_key': 'e_global_rel', + 'comm_world': MPI.COMM_WORLD, + **kwargs, + } + params['num_procs'] = min(params['comm_world'].size, 5) + params['plotting'] = params['comm_world'].rank == 0 + + configurations = get_configs(mode=f'vdp_stiffness-{mus[i]}', problem=run_vdp) + execute_configurations(**params, ax=axs.flatten()[i], decorate=True, configurations=configurations) + axs.flatten()[i].set_title(rf'$\mu={{{mus[i]}}}$') + + fig.suptitle('Van der Pol') + if params['comm_world'].rank == 0: + save_fig( + fig=fig, + name='vdp-stiffness', + work_key=params['work_key'], + precision_key=params['precision_key'], + legend=False, + base_path=base_path, + format=format, + ) + + +if __name__ == "__main__": + comm_world = MPI.COMM_WORLD + + params = { + 'mode': 'parallel_efficiency', + 'runs': 1, + 'num_procs': min(comm_world.size, 5), + 'plotting': comm_world.rank == 0, + } + params_single = { + **params, + 'problem': run_Schroedinger, + } + record = False + single_problem(**params_single, work_key='t', precision_key='e_global_rel', record=record) + # single_problem(**params_single, work_key='k_Newton_no_restart', precision_key='e_global_rel', record=False) + # single_problem(**params_single, work_key='param', precision_key='e_global_rel', record=False) + # ODEs(**params, work_key='t', precision_key='e_global_rel', record=record) + + all_params = { + 'record': False, + 'runs': 1, + 'work_key': 't', + 'precision_key': 'e_global_rel', + 'plotting': comm_world.rank == 0, + } + + for _mode in ['parallel_efficiency']: # , 'preconditioners', 'compare_adaptivity']: + # all_problems(**all_params, mode=mode) + comm_world.Barrier() + + if comm_world.rank == 0: + # parallel_efficiency(**params_single, work_key='k_SDC', precision_key='e_global_rel') + plt.show() diff --git a/pySDC/tests/test_projects/test_resilience/test_fault_injection.py b/pySDC/tests/test_projects/test_resilience/test_fault_injection.py index 94e252c4cf..48586eec28 100644 --- a/pySDC/tests/test_projects/test_resilience/test_fault_injection.py +++ b/pySDC/tests/test_projects/test_resilience/test_fault_injection.py @@ -60,7 +60,7 @@ def test_complex_conversion(): injector = FaultInjector() num_tests = int(1e3) - for i in range(num_tests): + for _i in range(num_tests): rand_complex = get_random_float() + get_random_float() * 1j # convert to bytes and back @@ -144,20 +144,13 @@ def test_fault_injection(): @pytest.mark.mpi4py +@pytest.mark.slow @pytest.mark.parametrize("numprocs", [5]) def test_fault_stats(numprocs): """ Test generation of fault statistics and their recovery rates """ import numpy as np - from pySDC.projects.Resilience.fault_stats import ( - FaultStats, - BaseStrategy, - AdaptivityStrategy, - IterateStrategy, - HotRodStrategy, - run_vdp, - ) # Set python path once my_env = os.environ.copy() @@ -174,27 +167,35 @@ def test_fault_stats(numprocs): numprocs, ) - vdp_stats = generate_stats(True) + stats = generate_stats(True) # test number of possible combinations for faults + expected_max_combinations = 3840 assert ( - vdp_stats.get_max_combinations() == 1536 - ), f"Expected 1536 possible combinations for faults in van der Pol problem, but got {vdp_stats.get_max_combinations()}!" + stats.get_max_combinations() == expected_max_combinations + ), f"Expected {expected_max_combinations} possible combinations for faults in van der Pol problem, but got {stats.get_max_combinations()}!" recovered_reference = { 'base': 1, 'adaptivity': 2, 'iterate': 1, 'Hot Rod': 2, + 'adaptivity_coll': 0, + 'double_adaptivity': 0, } - vdp_stats.get_recovered() + stats.get_recovered() - for strategy in vdp_stats.strategies: - dat = vdp_stats.load(strategy, True) - fixable_mask = vdp_stats.get_fixable_faults_only(strategy) - recovered_mask = vdp_stats.get_mask(strategy=strategy, key='recovered', op='eq', val=True) + for strategy in stats.strategies: + dat = stats.load(strategy=strategy, faults=True) + fixable_mask = stats.get_fixable_faults_only(strategy) + recovered_mask = stats.get_mask(strategy=strategy, key='recovered', op='eq', val=True) + index = stats.get_index(mask=fixable_mask) assert all(fixable_mask[:-1] == [False, True, False]), "Error in generating mask of fixable faults" + assert all(index == [1, 3]), "Error when converting to index" + + combinations = np.array(stats.get_combination_counts(dat, keys=['bit'], mask=fixable_mask)) + assert all(combinations == [1.0, 1.0]), "Error when counting combinations" recovered = len(dat['recovered'][recovered_mask]) crashed = len(dat['error'][dat['error'] == np.inf]) # on some systems the last run crashes... @@ -213,30 +214,44 @@ def generate_stats(load=False): Returns: Object containing the stats """ - from pySDC.projects.Resilience.fault_stats import ( - FaultStats, + from pySDC.projects.Resilience.strategies import ( BaseStrategy, AdaptivityStrategy, IterateStrategy, HotRodStrategy, - run_vdp, + AdaptivityCollocationRefinementStrategy, + DoubleAdaptivityStrategy, + AdaptivityAvoidRestartsStrategy, + AdaptivityInterpolationStrategy, + ) + from pySDC.projects.Resilience.fault_stats import ( + FaultStats, ) - import matplotlib.pyplot as plt + from pySDC.projects.Resilience.Lorenz import run_Lorenz np.seterr(all='warn') # get consistent behaviour across platforms - vdp_stats = FaultStats( - prob=run_vdp, + stats = FaultStats( + prob=run_Lorenz, faults=[False, True], reload=load, recovery_thresh=1.1, num_procs=1, mode='random', - strategies=[BaseStrategy(), AdaptivityStrategy(), IterateStrategy(), HotRodStrategy()], + strategies=[ + BaseStrategy(), + AdaptivityStrategy(), + IterateStrategy(), + HotRodStrategy(), + AdaptivityCollocationRefinementStrategy(), + DoubleAdaptivityStrategy(), + AdaptivityAvoidRestartsStrategy(), + AdaptivityInterpolationStrategy(), + ], stats_path='data', ) - vdp_stats.run_stats_generation(runs=4, step=2) - return vdp_stats + stats.run_stats_generation(runs=4, step=2) + return stats if __name__ == "__main__": diff --git a/pySDC/tests/test_projects/test_resilience/test_leaky_superconductor.py b/pySDC/tests/test_projects/test_resilience/test_leaky_superconductor.py deleted file mode 100644 index 4b41d366c4..0000000000 --- a/pySDC/tests/test_projects/test_resilience/test_leaky_superconductor.py +++ /dev/null @@ -1,12 +0,0 @@ -import pytest - - -@pytest.mark.base -@pytest.mark.parametrize('leak_type', ['linear', 'exponential']) -def test_imex_vs_fully_implicit_leaky_superconductor(leak_type): - """ - Test if the IMEX and fully implicit schemes get the same solution and that the runaway process has started. - """ - from pySDC.projects.Resilience.leaky_superconductor import compare_imex_full - - compare_imex_full(plotting=False, leak_type=leak_type) diff --git a/pySDC/tests/test_projects/test_resilience/test_order.py b/pySDC/tests/test_projects/test_resilience/test_order.py index 60138ee16a..52d5f4b5d4 100644 --- a/pySDC/tests/test_projects/test_resilience/test_order.py +++ b/pySDC/tests/test_projects/test_resilience/test_order.py @@ -2,7 +2,21 @@ @pytest.mark.base -def test_main(): - from pySDC.projects.Resilience.accuracy_check import main +@pytest.mark.parametrize("ks", [[2], [3], [4]]) +@pytest.mark.parametrize("serial", [True, False]) +def test_order_fixed_step_size(ks, serial): + from pySDC.projects.Resilience.accuracy_check import plot_all_errors, plt - main() + fig, ax = plt.subplots() + plot_all_errors(ax, ks, serial, Tend_fixed=1.0) + + +@pytest.mark.base +@pytest.mark.parametrize("ks", [[2], [3]]) +@pytest.mark.parametrize("serial", [True, False]) +def test_order_adaptive_step_size(ks, serial): + print(locals()) + from pySDC.projects.Resilience.accuracy_check import plot_all_errors, plt + + fig, ax = plt.subplots() + plot_all_errors(ax, ks, serial, Tend_fixed=5e-1, var='e_tol', dt_list=[1e-5, 5e-6], avoid_restarts=False) diff --git a/pySDC/tests/test_projects/test_resilience/test_quench.py b/pySDC/tests/test_projects/test_resilience/test_quench.py new file mode 100644 index 0000000000..d020e2c240 --- /dev/null +++ b/pySDC/tests/test_projects/test_resilience/test_quench.py @@ -0,0 +1,26 @@ +import pytest + + +@pytest.mark.base +@pytest.mark.parametrize('leak_type', ['linear', 'exponential']) +def test_imex_vs_fully_implicit_quench(leak_type): + """ + Test if the IMEX and fully implicit schemes get the same solution and that the runaway process has started. + """ + from pySDC.projects.Resilience.quench import compare_imex_full + + compare_imex_full(plotting=False, leak_type=leak_type) + + +@pytest.mark.base +def test_crossing_time_computation(): + from pySDC.projects.Resilience.quench import run_quench, get_crossing_time + + controller_params = {'logger_level': 30} + description = {'level_params': {'dt': 2.5e1}, 'step_params': {'maxiter': 5}} + stats, controller, _ = run_quench( + custom_controller_params=controller_params, + custom_description=description, + Tend=400, + ) + _ = get_crossing_time(stats, controller, num_points=5, inter_points=155)