From ca48726167cce9a9e057cb7cfde0c51228204d0a Mon Sep 17 00:00:00 2001 From: Thomas Baumann <39156931+brownbaerchen@users.noreply.github.com> Date: Fri, 8 Dec 2023 09:31:21 +0100 Subject: [PATCH] Allen-Cahn work-precision and resilience plots (#385) --- pySDC/core/Problem.py | 10 +- .../adaptivity.py | 3 + .../convergence_controller_classes/crash.py | 2 +- .../problem_classes/AllenCahn_2D_FFT.py | 12 +- pySDC/projects/Resilience/AC.py | 346 ++++++++++++++++++ pySDC/projects/Resilience/README.rst | 20 + pySDC/projects/Resilience/fault_stats.py | 7 +- pySDC/projects/Resilience/paper_plots.py | 39 +- pySDC/projects/Resilience/strategies.py | 74 +++- pySDC/projects/Resilience/work_precision.py | 105 +++--- .../test_polynomial_error.py | 2 +- .../test_resilience/test_strategies.py | 14 +- 12 files changed, 540 insertions(+), 94 deletions(-) create mode 100644 pySDC/projects/Resilience/AC.py diff --git a/pySDC/core/Problem.py b/pySDC/core/Problem.py index 72ca45556e..c9a6f87322 100644 --- a/pySDC/core/Problem.py +++ b/pySDC/core/Problem.py @@ -125,11 +125,13 @@ def generate_scipy_reference_solution(self, eval_rhs, t, u_init=None, t_init=Non import numpy as np from scipy.integrate import solve_ivp - tol = 100 * np.finfo(float).eps + kwargs = { + 'atol': 100 * np.finfo(float).eps, + 'rtol': 100 * np.finfo(float).eps, + **kwargs, + } u_init = self.u_exact(t=0) if u_init is None else u_init * 1.0 t_init = 0 if t_init is None else t_init u_shape = u_init.shape - return ( - solve_ivp(eval_rhs, (t_init, t), u_init.flatten(), rtol=tol, atol=tol, **kwargs).y[:, -1].reshape(u_shape) - ) + return solve_ivp(eval_rhs, (t_init, t), u_init.flatten(), **kwargs).y[:, -1].reshape(u_shape) diff --git a/pySDC/implementations/convergence_controller_classes/adaptivity.py b/pySDC/implementations/convergence_controller_classes/adaptivity.py index e9f74f8ec1..f4efe12211 100644 --- a/pySDC/implementations/convergence_controller_classes/adaptivity.py +++ b/pySDC/implementations/convergence_controller_classes/adaptivity.py @@ -59,6 +59,9 @@ def dependencies(self, controller, description, **kwargs): step_limiter_params = {key: self.params.__dict__[key] for key in available_keys} controller.add_convergence_controller(StepSizeLimiter, params=step_limiter_params, description=description) + if self.params.useMPI: + self.prepare_MPI_logical_operations() + return None def get_new_step_size(self, controller, S, **kwargs): diff --git a/pySDC/implementations/convergence_controller_classes/crash.py b/pySDC/implementations/convergence_controller_classes/crash.py index 521935136d..f9729e89a3 100644 --- a/pySDC/implementations/convergence_controller_classes/crash.py +++ b/pySDC/implementations/convergence_controller_classes/crash.py @@ -85,7 +85,7 @@ def prepare_next_block(self, controller, S, *args, **kwargs): for u in lvl.u: if u is None: break - isfinite = all(np.isfinite(u.flatten())) + isfinite = np.all(np.isfinite(u)) below_limit = abs(u) < self.params.thresh diff --git a/pySDC/implementations/problem_classes/AllenCahn_2D_FFT.py b/pySDC/implementations/problem_classes/AllenCahn_2D_FFT.py index b000a21472..9cf1a3da52 100644 --- a/pySDC/implementations/problem_classes/AllenCahn_2D_FFT.py +++ b/pySDC/implementations/problem_classes/AllenCahn_2D_FFT.py @@ -59,11 +59,19 @@ class allencahn2d_imex(ptype): dtype_u = mesh dtype_f = imex_mesh - def __init__(self, nvars=None, nu=2, eps=0.04, radius=0.25, L=1.0, init_type='circle'): + def __init__( + self, + nvars=None, + nu=2, + eps=0.04, + radius=0.25, + L=1.0, + init_type='circle', + ): """Initialization routine""" if nvars is None: - nvars = [(256, 256), (64, 64)] + nvars = (128, 128) # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on if len(nvars) != 2: diff --git a/pySDC/projects/Resilience/AC.py b/pySDC/projects/Resilience/AC.py new file mode 100644 index 0000000000..ff77e6266f --- /dev/null +++ b/pySDC/projects/Resilience/AC.py @@ -0,0 +1,346 @@ +# script to run an Allen-Cahn problem +from pySDC.implementations.problem_classes.AllenCahn_2D_FD import allencahn_fullyimplicit, allencahn_semiimplicit +from pySDC.implementations.problem_classes.AllenCahn_2D_FFT import allencahn2d_imex +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI +from pySDC.core.Hooks import hooks +from pySDC.projects.Resilience.hook import hook_collection, LogData +from pySDC.projects.Resilience.strategies import merge_descriptions +from pySDC.projects.Resilience.sweepers import imex_1st_order_efficient, generic_implicit_efficient +import matplotlib.pyplot as plt +import numpy as np + +from pySDC.core.Errors import ConvergenceError + + +def run_AC( + custom_description=None, + num_procs=1, + Tend=1e-2, + hook_class=LogData, + fault_stuff=None, + custom_controller_params=None, + imex=False, + u0=None, + t0=None, + use_MPI=False, + live_plot=False, + FFT=True, + **kwargs, +): + """ + 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 + bool: If the code crashed + """ + if custom_description is not None: + problem_params = custom_description.get('problem_params', {}) + if 'imex' in problem_params.keys(): + imex = problem_params['imex'] + problem_params.pop('imex', None) + if 'FFT' in problem_params.keys(): + FFT = problem_params['FFT'] + problem_params.pop('FFT', None) + + # import problem and sweeper class + if FFT: + from pySDC.implementations.problem_classes.AllenCahn_2D_FFT import allencahn2d_imex as problem_class + from pySDC.projects.Resilience.sweepers import imex_1st_order_efficient as sweeper_class + elif imex: + from pySDC.implementations.problem_classes.AllenCahn_2D_FD import allencahn_semiimplicit as problem_class + from pySDC.projects.Resilience.sweepers import imex_1st_order_efficient as sweeper_class + else: + from pySDC.implementations.problem_classes.AllenCahn_2D_FD import allencahn_fullyimplicit as problem_class + from pySDC.projects.Resilience.sweepers import generic_implicit_efficient as sweeper_class + + level_params = {} + level_params['dt'] = 1e-4 + level_params['restol'] = 1e-8 + + sweeper_params = {} + sweeper_params['quad_type'] = 'RADAU-RIGHT' + sweeper_params['num_nodes'] = 3 + sweeper_params['QI'] = 'LU' + sweeper_params['QE'] = 'PIC' + + # problem params + fd_params = { + 'newton_tol': 1e-9, + 'order': 2, + } + problem_params = { + 'nvars': (128, 128), + 'init_type': 'circle', + } + if not FFT: + problem_params = {**problem_params, **fd_params} + + step_params = {} + step_params['maxiter'] = 5 + + controller_params = {} + controller_params['logger_level'] = 30 + controller_params['hook_class'] = ( + hook_collection + (hook_class if type(hook_class) == list else [hook_class]) + ([LivePlot] if live_plot else []) + ) + controller_params['mssdc_jac'] = False + + if custom_controller_params is not None: + controller_params = {**controller_params, **custom_controller_params} + + description = {} + description['problem_class'] = problem_class + description['problem_params'] = problem_params + description['sweeper_class'] = sweeper_class + 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) + + t0 = 0.0 if t0 is None else t0 + + 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 + + if fault_stuff is not None: + from pySDC.projects.Resilience.fault_injection import prepare_controller_for_faults + + prepare_controller_for_faults(controller, fault_stuff) + + crash = False + try: + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + except ConvergenceError as e: + print(f'Warning: Premature termination!: {e}') + stats = controller.return_stats() + crash = True + return stats, controller, crash + + +def plot_solution(stats): # pragma: no cover + import matplotlib.pyplot as plt + from pySDC.helpers.stats_helper import get_sorted + + fig, ax = plt.subplots(1, 1) + + u = get_sorted(stats, type='u', recomputed=False) + for me in u: # pun intended + ax.imshow(me[1], vmin=-1, vmax=1) + ax.set_title(f't={me[0]:.2e}') + plt.pause(1e-1) + + plt.show() + + +class LivePlot(hooks): # pragma: no cover + def __init__(self): + super().__init__() + self.fig, self.axs = plt.subplots(1, 3, figsize=(12, 4)) + self.radius = [] + self.exact_radius = [] + self.t = [] + self.dt = [] + + def post_step(self, step, level_number): + super().post_step(step, level_number) + L = step.levels[level_number] + self.t += [step.time + step.dt] + + # plot solution + self.axs[0].cla() + if len(L.uend.shape) > 1: + self.axs[0].imshow(L.uend, vmin=0.0, vmax=1.0) + + # plot radius + self.axs[1].cla() + radius, _ = LogRadius.compute_radius(step.levels[level_number]) + exact_radius = LogRadius.exact_radius(step.levels[level_number]) + + self.radius += [radius] + self.exact_radius += [exact_radius] + self.axs[1].plot(self.t, self.exact_radius, label='exact') + self.axs[1].plot(self.t, self.radius, label='numerical') + self.axs[1].set_ylim([0, 0.26]) + self.axs[1].set_xlim([0, 0.03]) + self.axs[1].legend(frameon=False) + self.axs[1].set_title(r'Radius') + else: + self.axs[0].plot(L.prob.xvalues, L.prob.u_exact(t=L.time + L.dt), label='exact') + self.axs[0].plot(L.prob.xvalues, L.uend, label='numerical') + self.axs[0].set_title(f't = {step.time + step.dt:.2e}') + + # plot step size + self.axs[2].cla() + self.dt += [step.dt] + self.axs[2].plot(self.t, self.dt) + self.axs[2].set_yscale('log') + self.axs[2].axhline(step.levels[level_number].prob.eps ** 2, label=r'$\epsilon^2$', color='black', ls='--') + self.axs[2].legend(frameon=False) + self.axs[2].set_xlim([0, 0.03]) + self.axs[2].set_title(r'$\Delta t$') + + if step.status.restart: + for me in [self.radius, self.exact_radius, self.t, self.dt]: + try: + me.pop(-1) + except (TypeError, IndexError): + pass + + plt.pause(1e-9) + + +class LogRadius(hooks): + @staticmethod + def compute_radius(L): + c = np.count_nonzero(L.u[0] > 0.0) + radius = np.sqrt(c / np.pi) * L.prob.dx + + rows, cols = np.where(L.u[0] > 0.0) + + rows1 = np.where(L.u[0][int((L.prob.init[0][0]) / 2), : int((L.prob.init[0][0]) / 2)] > -0.99) + rows2 = np.where(L.u[0][int((L.prob.init[0][0]) / 2), : int((L.prob.init[0][0]) / 2)] < 0.99) + interface_width = (rows2[0][-1] - rows1[0][0]) * L.prob.dx / L.prob.eps + + return radius, interface_width + + @staticmethod + def exact_radius(L): + init_radius = L.prob.radius + return np.sqrt(max(init_radius**2 - 2.0 * (L.time + L.dt), 0)) + + def pre_run(self, step, level_number): + """ + Overwrite standard pre run hook + + Args: + step (pySDC.Step.step): the current step + level_number (int): the current level number + """ + super().pre_run(step, level_number) + L = step.levels[0] + + radius, interface_width = self.compute_radius(L) + exact_radius = self.exact_radius(L) + + if L.time == 0.0: + self.add_to_stats( + process=step.status.slot, + time=L.time, + level=-1, + iter=step.status.iter, + sweep=L.status.sweep, + type='computed_radius', + value=radius, + ) + self.add_to_stats( + process=step.status.slot, + time=L.time, + level=-1, + iter=step.status.iter, + sweep=L.status.sweep, + type='exact_radius', + value=exact_radius, + ) + self.add_to_stats( + process=step.status.slot, + time=L.time, + level=-1, + iter=step.status.iter, + sweep=L.status.sweep, + type='interface_width', + value=interface_width, + ) + + def post_run(self, step, level_number): + """ + Args: + step (pySDC.Step.step): the current step + level_number (int): the current level number + """ + super().post_run(step, level_number) + + L = step.levels[0] + + exact_radius = self.exact_radius(L) + radius, interface_width = self.compute_radius(L) + + self.add_to_stats( + process=step.status.slot, + time=L.time + L.dt, + level=-1, + iter=step.status.iter, + sweep=L.status.sweep, + type='computed_radius', + value=radius, + ) + self.add_to_stats( + process=step.status.slot, + time=L.time + L.dt, + level=-1, + iter=step.status.iter, + sweep=L.status.sweep, + type='exact_radius', + value=exact_radius, + ) + self.add_to_stats( + process=step.status.slot, + time=L.time + L.dt, + level=-1, + iter=step.status.iter, + sweep=L.status.sweep, + type='interface_width', + value=interface_width, + ) + self.add_to_stats( + process=step.status.slot, + time=L.time + L.dt, + level=level_number, + iter=step.status.iter, + sweep=L.status.sweep, + type='e_global_post_run', + value=abs(radius - exact_radius), + ) + self.add_to_stats( + process=step.status.slot, + time=L.time + L.dt, + level=level_number, + iter=step.status.iter, + sweep=L.status.sweep, + type='e_global_rel_post_run', + value=abs(radius - exact_radius) / abs(exact_radius), + ) + + +if __name__ == '__main__': + from pySDC.implementations.hooks.log_errors import LogLocalErrorPostStep + + stats, _, _ = run_AC(imex=True, hook_class=LogLocalErrorPostStep) + plot_solution(stats) diff --git a/pySDC/projects/Resilience/README.rst b/pySDC/projects/Resilience/README.rst index db5fdedcb0..9e292c0db6 100644 --- a/pySDC/projects/Resilience/README.rst +++ b/pySDC/projects/Resilience/README.rst @@ -39,3 +39,23 @@ You can see the results below, except for the solution, which looks the same as :width: 24% .. image:: ../../../data/error_estimate_order_parallel.png :width: 23% + + +Reproduction of the plots in the adaptive SDC paper +--------------------------------------------------- +To reproduce the plots you need to install pySDC with all packages in the mpi4py environment. +Then, navigate to this directory, `pySDC/projects/Resilience/` and run the following commands: + + +```bash + +mpirun -np 4 python work_precision.py +mpirun -np 4 python fault_stats.py prob run_vdp +mpirun -np 4 python fault_stats.py prob run_quench +mpirun -np 4 python fault_stats.py prob run_AC +mpirun -np 4 python fault_stats.py prob run_Schroedinger +python paper_plots.py + +``` + +Possibly, you need to create some directories in this one to store and load things, if path errors occur. diff --git a/pySDC/projects/Resilience/fault_stats.py b/pySDC/projects/Resilience/fault_stats.py index 4558002ea1..f28aa11746 100644 --- a/pySDC/projects/Resilience/fault_stats.py +++ b/pySDC/projects/Resilience/fault_stats.py @@ -20,6 +20,7 @@ from pySDC.projects.Resilience.Lorenz import run_Lorenz from pySDC.projects.Resilience.Schroedinger import run_Schroedinger from pySDC.projects.Resilience.quench import run_quench +from pySDC.projects.Resilience.AC import run_AC from pySDC.projects.Resilience.strategies import BaseStrategy, AdaptivityStrategy, IterateStrategy, HotRodStrategy import logging @@ -629,6 +630,8 @@ def get_name(self, strategy=None, faults=True, mode=None): prob_name = 'Schroedinger' elif self.prob == run_quench: prob_name = 'Quench' + elif self.prob.__name__ == 'run_AC': + prob_name = 'Allen-Cahn' else: raise NotImplementedError(f'Name not implemented for problem {self.prob}') @@ -1651,11 +1654,11 @@ def compare_adaptivity_modes(): def main(): kwargs = { - 'prob': run_quench, + 'prob': run_AC, 'num_procs': 1, 'mode': 'default', 'runs': 2000, - 'reload': False, + 'reload': True, **parse_args(), } diff --git a/pySDC/projects/Resilience/paper_plots.py b/pySDC/projects/Resilience/paper_plots.py index 90e73020f1..9804ebc2b9 100644 --- a/pySDC/projects/Resilience/paper_plots.py +++ b/pySDC/projects/Resilience/paper_plots.py @@ -8,6 +8,7 @@ run_Schroedinger, run_vdp, run_quench, + run_AC, RECOVERY_THRESH_ABS, ) from pySDC.projects.Resilience.strategies import ( @@ -196,11 +197,11 @@ def compare_recovery_rate_problems(**kwargs): # pragma: no cover """ stats = [ get_stats(run_vdp, **kwargs), - get_stats(run_Lorenz, **kwargs), - get_stats(run_Schroedinger, **kwargs), get_stats(run_quench, **kwargs), + get_stats(run_Schroedinger, **kwargs), + get_stats(run_AC, **kwargs), ] - titles = ['Van der Pol', 'Lorenz attractor', r'Schr\"odinger', 'Quench'] + titles = ['Van der Pol', 'Quench', r'Schr\"odinger', 'Allen-Cahn'] my_setup_mpl() fig, axs = plt.subplots(2, 2, figsize=figsize_by_journal(JOURNAL, 1, 0.8), sharey=True) @@ -409,6 +410,37 @@ def plot_quench_solution(): # pragma: no cover savefig(fig, 'quench_sol') +def plot_AC_solution(): # pragma: no cover + from pySDC.projects.TOMS.AllenCahn_monitor import monitor + + my_setup_mpl() + if JOURNAL == 'JSC_beamer': + raise NotImplementedError + fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.5, 0.9)) + else: + fig, axs = plt.subplots(1, 2, figsize=figsize_by_journal(JOURNAL, 1.0, 0.45)) + + stats, _, _ = run_AC(Tend=0.032, hook_class=monitor) + + u = get_sorted(stats, type='u') + + computed_radius = get_sorted(stats, type='computed_radius') + exact_radius = get_sorted(stats, type='exact_radius') + axs[1].plot([me[0] for me in computed_radius], [me[1] for me in computed_radius], ls='-', label='numerical') + axs[1].plot([me[0] for me in exact_radius], [me[1] for me in exact_radius], ls='--', color='black', label='exact') + axs[1].axvline(0.025, ls=':', label=r'$t=0.025$', color='grey') + axs[1].set_title('Radius over time') + axs[1].set_xlabel('$t$') + axs[1].legend(frameon=False) + + im = axs[0].imshow(u[0][1], extent=(-0.5, 0.5, -0.5, 0.5)) + fig.colorbar(im) + axs[0].set_title(r'$u_0$') + axs[0].set_xlabel('$x$') + axs[0].set_ylabel('$y$') + savefig(fig, 'AC_sol') + + def plot_vdp_solution(): # pragma: no cover """ Plot the solution of van der Pol problem over time to illustrate the varying time scales. @@ -545,6 +577,7 @@ def make_plots_for_paper(): # pragma: no cover work_precision() plot_vdp_solution() + plot_AC_solution() plot_quench_solution() plot_recovery_rate(get_stats(run_vdp)) diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py index 8b235251a8..ecc98c5203 100644 --- a/pySDC/projects/Resilience/strategies.py +++ b/pySDC/projects/Resilience/strategies.py @@ -128,6 +128,8 @@ def get_fault_args(self, problem, num_procs): args['time'] = 41.0 elif problem.__name__ == "run_Lorenz": args['time'] = 0.3 + elif problem.__name__ == "run_AC": + args['time'] = 1e-2 return args @@ -148,7 +150,7 @@ def get_random_params(self, problem, num_procs): rnd_params['iteration'] = base_params['step_params']['maxiter'] rnd_params['rank'] = num_procs - if problem.__name__ in ['run_Schroedinger', 'run_quench']: + if problem.__name__ in ['run_Schroedinger', 'run_quench', 'run_AC']: rnd_params['min_node'] = 1 if problem.__name__ == "run_quench": @@ -203,7 +205,7 @@ def get_Tend(cls, problem, num_procs=1): elif problem.__name__ == "run_quench": return 500.0 elif problem.__name__ == "run_AC": - return 1e-2 + return 0.025 else: raise NotImplementedError('I don\'t have a final time for your problem!') @@ -256,9 +258,16 @@ def get_base_parameters(self, problem, num_procs=1): 'order': 6, } elif problem.__name__ == "run_AC": - custom_description['level_params'] = {'restol': -1, 'dt': 1e-4} + eps = 4e-2 custom_description['step_params'] = {'maxiter': 5} - custom_description['problem_params'] = {'newton_maxiter': 99, 'newton_tol': 1e-11} + custom_description['problem_params'] = { + 'nvars': (128,) * 2, + 'init_type': 'circle', + 'eps': eps, + 'radius': 0.25, + 'nu': 2, + } + custom_description['level_params'] = {'restol': -1, 'dt': 0.5 * eps**2} custom_description['convergence_controllers'] = { # StepSizeLimiter: {'dt_min': self.get_Tend(problem=problem, num_procs=num_procs) / self.max_steps} @@ -276,6 +285,7 @@ def get_base_parameters(self, problem, num_procs=1): 'run_Lorenz': 60, 'run_Schroedinger': 150, 'run_quench': 150, + 'run_AC': 150, } custom_description['convergence_controllers'][StopAtMaxRuntime] = { @@ -360,7 +370,7 @@ def get_custom_description(self, problem, num_procs=1): 'maxiter': 15, } - if self.newton_inexactness and problem.__name__ not in ['run_Schroedinger']: + if self.newton_inexactness and problem.__name__ not in ['run_Schroedinger', 'run_AC']: if problem.__name__ == 'run_quench': inexactness_params['ratio'] = 1e-1 inexactness_params['min_tol'] = 1e-11 @@ -370,7 +380,7 @@ def get_custom_description(self, problem, num_procs=1): if problem.__name__ in ['run_vdp']: desc['problem_params']['stop_at_nan'] = False - if self.linear_inexactness and problem.__name__ in ['run_AC', 'run_quench']: + if self.linear_inexactness and problem.__name__ in ['run_quench']: desc['problem_params']['inexact_linear_ratio'] = 1e-1 if problem.__name__ in ['run_quench']: desc['problem_params']['direct_solver'] = False @@ -407,6 +417,12 @@ def __init__(self, skip_residual_computation='all', **kwargs): def label(self): return r'fixed' + def get_custom_description(self, problem, num_procs): + desc = super().get_custom_description(problem, num_procs) + if problem.__name__ == "run_AC": + desc['level_params']['dt'] = 0.8 * desc['problem_params']['eps'] ** 2 / 8.0 + return desc + def get_custom_description_for_faults(self, problem, *args, **kwargs): desc = self.get_custom_description(problem, *args, **kwargs) if problem.__name__ == "run_quench": @@ -492,6 +508,7 @@ def get_custom_description(self, problem, num_procs): from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter + base_params = super().get_custom_description(problem, num_procs) custom_description = {} custom_description['convergence_controllers'] = {} @@ -505,12 +522,12 @@ def get_custom_description(self, problem, num_procs): elif problem.__name__ == "run_Lorenz": e_tol = 2e-5 elif problem.__name__ == "run_Schroedinger": - e_tol = 4e-6 + e_tol = 4e-7 elif problem.__name__ == "run_quench": e_tol = 1e-8 custom_description['problem_params'] = { - 'newton_tol': 1e-9, - 'lintol': 1e-10, + 'newton_tol': 1e-10, + 'lintol': 1e-11, } from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestarting @@ -519,7 +536,9 @@ def get_custom_description(self, problem, num_procs): 'max_restarts': 99, } elif problem.__name__ == "run_AC": - e_tol = 1e-4 + e_tol = 1e-6 + dt_max = 0.8 * base_params['problem_params']['eps'] ** 2 + else: raise NotImplementedError( 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\ @@ -533,7 +552,7 @@ def get_custom_description(self, problem, num_procs): custom_description['convergence_controllers'][StepSizeLimiter] = { 'dt_max': dt_max, } - return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + return merge_descriptions(base_params, custom_description) def get_reference_value(self, problem, key, op, num_procs=1): """ @@ -723,6 +742,8 @@ def get_custom_description(self, problem, num_procs): restol = 6.5e-7 elif problem.__name__ == "run_quench": restol = 1e-7 + elif problem.__name__ == "run_AC": + restol = 1e-11 else: raise NotImplementedError( 'I don\'t have a residual tolerance for your problem. Please add one to the \ @@ -791,12 +812,17 @@ def get_custom_description(self, problem, num_procs, *args, **kwargs): desc['problem_params']['newton_tol'] = 1e-9 desc['problem_params']['lintol'] = 1e-9 desc['level_params']['dt'] = 2.5 + elif problem.__name__ == "run_AC": + desc['level_params']['restol'] = 1e-11 + desc['level_params']['dt'] = 0.8 * desc['problem_params']['eps'] ** 2 / 8.0 return desc def get_custom_description_for_faults(self, problem, *args, **kwargs): desc = self.get_custom_description(problem, *args, **kwargs) if problem.__name__ == 'run_quench': desc['level_params']['dt'] = 5.0 + elif problem.__name__ == 'run_AC': + desc['level_params']['dt'] = 0.6 * desc['problem_params']['eps'] ** 2 return desc def get_reference_value(self, problem, key, op, num_procs=1): @@ -853,6 +879,7 @@ def get_custom_description(self, problem, num_procs): from pySDC.implementations.convergence_controller_classes.hotrod import HotRod from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestartingNonMPI + base_params = super().get_custom_description(problem, num_procs) if problem.__name__ == "run_vdp": if num_procs == 4: HotRod_tol = 1.800804e-04 @@ -886,6 +913,9 @@ def get_custom_description(self, problem, num_procs): else: HotRod_tol = 5.198620e-04 maxiter = 6 + elif problem.__name__ == 'run_AC': + HotRod_tol = 9.564437e-06 + maxiter = 6 else: raise NotImplementedError( 'I don\'t have a tolerance for Hot Rod for your problem. Please add one to the\ @@ -904,9 +934,11 @@ def get_custom_description(self, problem, num_procs): }, }, 'step_params': {'maxiter': maxiter}, + 'level_params': {}, } - - return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + if problem.__name__ == "run_AC": + custom_description['level_params']['dt'] = 0.8 * base_params['problem_params']['eps'] ** 2 / 8.0 + return merge_descriptions(base_params, custom_description) def get_custom_description_for_faults(self, problem, *args, **kwargs): desc = self.get_custom_description(problem, *args, **kwargs) @@ -1300,7 +1332,7 @@ def get_reference_value(self, problem, key, op, num_procs=1): if key == 'work_newton' and op == sum: return 0 elif key == 'e_global_post_run' and op == max: - return 2.444605657738645e-07 + return 3.1786601531890356e-08 raise NotImplementedError('The reference value you are looking for is not implemented for this strategy!') @@ -1432,6 +1464,7 @@ def get_description_for_tolerance(self, problem, param, **kwargs): desc = {} if problem.__name__ == 'run_Schroedinger': desc['problem_params'] = {'lintol': param} + return desc @property @@ -1459,6 +1492,9 @@ def get_custom_description(self, problem, num_procs=1): desc = super().get_custom_description(problem, num_procs) desc['sweeper_class'] = Cash_Karp + + if problem.__name__ == "run_AC": + desc['level_params']['dt'] = 2e-5 return desc def get_reference_value(self, problem, key, op, num_procs=1): @@ -1809,6 +1845,7 @@ def get_custom_description(self, problem, num_procs): from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter + base_params = super().get_custom_description(problem, num_procs) custom_description = {} dt_max = np.inf @@ -1825,7 +1862,7 @@ def get_custom_description(self, problem, num_procs): elif problem.__name__ == "run_Lorenz": e_tol = 7e-4 elif problem.__name__ == "run_Schroedinger": - e_tol = 3e-4 + e_tol = 3e-5 elif problem.__name__ == "run_quench": e_tol = 1e-7 level_params['dt'] = 50.0 @@ -1833,6 +1870,7 @@ def get_custom_description(self, problem, num_procs): restol_rel = 1e-1 elif problem.__name__ == "run_AC": e_tol = 1e-4 + dt_max = 0.8 * base_params['problem_params']['eps'] ** 2 else: raise NotImplementedError( 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\ @@ -1854,7 +1892,7 @@ def get_custom_description(self, problem, num_procs): }, } custom_description['level_params'] = level_params - return merge_descriptions(super().get_custom_description(problem, num_procs), custom_description) + return merge_descriptions(base_params, custom_description) def get_custom_description_for_faults(self, problem, *args, **kwargs): desc = self.get_custom_description(problem, *args, **kwargs) @@ -1863,6 +1901,10 @@ def get_custom_description_for_faults(self, problem, *args, **kwargs): desc['convergence_controllers'][AdaptivityPolynomialError]['e_tol'] = 1e-7 * 11 desc['level_params']['dt'] = 4.0 + elif problem.__name__ == "run_AC": + from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError + + desc['convergence_controllers'][AdaptivityPolynomialError]['e_tol'] = 1e-3 return desc def get_random_params(self, problem, num_procs): diff --git a/pySDC/projects/Resilience/work_precision.py b/pySDC/projects/Resilience/work_precision.py index 7bde3d8463..09dc3acc9c 100644 --- a/pySDC/projects/Resilience/work_precision.py +++ b/pySDC/projects/Resilience/work_precision.py @@ -11,6 +11,7 @@ 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.projects.Resilience.AC import run_AC from pySDC.helpers.stats_helper import get_sorted, filter_stats from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal @@ -32,10 +33,12 @@ 'k_linear': ('work_linear', sum, None), 'k_Newton_no_restart': ('work_newton', sum, False), 'k_rhs': ('work_rhs', sum, None), + 'num_steps': ('dt', len, None), 'restart': ('restart', sum, None), 'dt_mean': ('dt', np.mean, False), 'dt_max': ('dt', max, False), 'dt_min': ('dt', min, False), + 'dt_sigma': ('dt', np.std, False), 'e_embedded_max': ('error_embedded_estimate', max, False), 'u0_increment_max': ('u0_increment', max, None), 'u0_increment_mean': ('u0_increment', np.mean, None), @@ -55,12 +58,8 @@ def get_forbidden_combinations(problem, strategy, **kwargs): problem (function): A problem to run strategy (Strategy): SDC strategy """ - if problem.__name__ == 'run_quench': - if strategy.name in ['ERK']: - return True - - if problem.__name__ == 'run_Schroedinger': - if strategy.name in ['ERK']: + if strategy.name == 'ERK': + if problem.__name__ in ['run_quench', 'run_Schroedinger', 'run_AC']: return True return False @@ -266,7 +265,7 @@ def record_work_precision( exponents = [-4, -3, -2, -1, 0, 1, 2] elif param == 'dt': power = 2.0 - exponents = [-1, 0, 1, 2, 3] + exponents = [-1, 0, 1, 2, 3][::-1] elif param == 'restol': power = 10.0 exponents = [-2, -1, 0, 1, 2, 3] @@ -285,13 +284,9 @@ def record_work_precision( elif param == 'dt': param_range = [1.25, 2.5, 5.0, 10.0, 20.0][::-1] - elif problem.__name__ == 'run_AC': - if param == 'e_tol': - param_range = [1e-4, 1e-5, 1e-6, 1e-7, 1e-8][::-1] - # run multiple times with different parameters for i in range(len(param_range)): - set_parameter(description, where, param_range[i] * strategy.precision_range_fac) + set_parameter(description, where, param_range[i]) data[param_range[i]] = {key: [] for key in MAPPINGS.keys()} data[param_range[i]]['param'] = [param_range[i]] @@ -537,6 +532,8 @@ def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only 'e_local_max': 'max. local error', 'restart': 'restarts', 'dt_max': r'$\Delta t_\mathrm{max}$', + 'dt_min': r'$\Delta t_\mathrm{min}$', + 'dt_sigma': r'$\sigma(\Delta t)$', 'dt_mean': r'$\bar{\Delta t}$', 'param': 'parameter', 'u0_increment_max': r'$\| \Delta u_0 \|_{\infty} $', @@ -559,6 +556,7 @@ def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only 'run_Lorenz': 'Lorenz attractor', 'run_Schroedinger': r'Schr\"odinger', 'run_quench': 'Quench', + 'run_AC': 'Allen-Cahn', } ax.set_title(titles.get(problem.__name__, '')) @@ -744,19 +742,21 @@ def get_configs(mode, problem): BaseStrategy, ) - configurations[1] = { - 'strategies': [AdaptivityPolynomialError(useMPI=True)], - } configurations[2] = { 'strategies': [kAdaptivityStrategy(useMPI=True)], } - + configurations[1] = { + 'strategies': [AdaptivityPolynomialError(useMPI=True)], + } configurations[0] = { 'custom_description': { 'step_params': {'maxiter': 5}, 'sweeper_params': {'num_nodes': 3, 'quad_type': 'RADAU-RIGHT'}, }, - 'strategies': [AdaptivityStrategy(useMPI=True), BaseStrategy(useMPI=True)], + 'strategies': [ + AdaptivityStrategy(useMPI=True), + BaseStrategy(useMPI=True), + ], } elif mode == 'interpolate_between_restarts': @@ -767,7 +767,7 @@ def get_configs(mode, problem): i = 0 for interpolate_between_restarts, handle, ls in zip( - [True, False], ['Interpolation bewteen restarts', 'regular'], ['--', '-'] + [True, False], ['Interpolation between restarts', 'regular'], ['--', '-'] ): configurations[i] = { 'strategies': [ @@ -814,7 +814,7 @@ def get_configs(mode, problem): """ from pySDC.projects.Resilience.strategies import AdaptivityStrategy, AdaptivityPolynomialError - if problem.__name__ in ['run_Schroedinger']: + if problem.__name__ in ['run_Schroedinger', 'run_AC']: from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper else: from pySDC.implementations.sweeper_classes.generic_implicit_MPI import ( @@ -1077,7 +1077,7 @@ def get_configs(mode, problem): AdaptivityPolynomialError, ) - if problem.__name__ in ['run_Schroedinger']: + if problem.__name__ in ['run_Schroedinger', 'run_AC']: from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper else: from pySDC.implementations.sweeper_classes.generic_implicit_MPI import ( @@ -1099,10 +1099,6 @@ def get_configs(mode, problem): 5: ':', } - desc_RK = {} - if problem.__name__ in ['run_Schroedinger']: - desc_RK['problem_params'] = {'imex': True} - configurations[3] = { 'custom_description': desc_poly, 'strategies': [AdaptivityPolynomialError(useMPI=True)], @@ -1113,18 +1109,19 @@ def get_configs(mode, problem): configurations[-1] = { 'strategies': [ ERKStrategy(useMPI=True), - ARKStrategy(useMPI=True) if problem.__name__ in ['run_Schroedinger'] else ESDIRKStrategy(useMPI=True), + ARKStrategy(useMPI=True) + if problem.__name__ in ['run_Schroedinger', 'run_AC'] + else ESDIRKStrategy(useMPI=True), ], 'num_procs': 1, - 'custom_description': desc_RK, } - configurations[2] = { 'strategies': [AdaptivityStrategy(useMPI=True)], 'custom_description': desc, 'num_procs': 4, 'plotting_params': {'label': r'$\Delta t$ adaptivity $N$=4x1'}, } + elif mode == 'RK_comp_high_order': """ Compare higher order SDC than we can get with RKM to RKM @@ -1307,7 +1304,7 @@ def all_problems(mode='compare_strategies', plotting=True, base_path='data', **k **kwargs, } - problems = [run_vdp, run_Lorenz, run_Schroedinger, run_quench] + problems = [run_vdp, run_quench, run_Schroedinger, run_AC] logger.log(26, f"Doing for all problems {mode}") for i in range(len(problems)): @@ -1498,46 +1495,38 @@ def aggregate_parallel_efficiency_plot(): # pragma: no cover if __name__ == "__main__": comm_world = MPI.COMM_WORLD - record = True - for mode in [ - # 'compare_strategies', - #'RK_comp_high_order', - # 'step_size_limiting', - # 'parallel_efficiency', - 'diagonal_SDC', - ]: - params = { - 'mode': mode, - 'runs': 1, - #'num_procs': 1, # min(comm_world.size, 5), - 'plotting': comm_world.rank == 0, - } - params_single = { - **params, - 'problem': run_Schroedinger, - } - single_problem(**params_single, work_key='t', precision_key='e_global_rel', record=record) + # record = False + # for mode in [ + # 'compare_strategies', + # 'RK_comp', + # 'parallel_efficiency', + # ]: + # params = { + # 'mode': mode, + # 'runs': 5, + # 'plotting': comm_world.rank == 0, + # } + # params_single = { + # **params, + # 'problem': run_AC, + # } + # single_problem(**params_single, work_key='t', precision_key='e_global_rel', record=record) all_params = { - 'record': False, - 'runs': 1, - 'work_key': 'k_linear', + 'record': True, + 'runs': 5, + 'work_key': 't', 'precision_key': 'e_global_rel', 'plotting': comm_world.rank == 0, - #'num_procs': 4, } for mode in [ - # 'RK_comp', - # 'parallel_efficiency', - # 'compare_adaptivity', - # 'compare_strategies', - # 'preconditioners', - # 'diagonal_SDC', + 'RK_comp', + 'parallel_efficiency', + 'compare_strategies', ]: 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_convergence_controllers/test_polynomial_error.py b/pySDC/tests/test_convergence_controllers/test_polynomial_error.py index a19a1c3808..a9244275e1 100644 --- a/pySDC/tests/test_convergence_controllers/test_polynomial_error.py +++ b/pySDC/tests/test_convergence_controllers/test_polynomial_error.py @@ -3,7 +3,7 @@ def get_controller(dt, num_nodes, quad_type, useMPI): """ - Runs a single advection problem with certain parameters + Get a controller prepared for polynomial test equation Args: dt (float): Step size diff --git a/pySDC/tests/test_projects/test_resilience/test_strategies.py b/pySDC/tests/test_projects/test_resilience/test_strategies.py index 1812f217cf..5e797e3e34 100644 --- a/pySDC/tests/test_projects/test_resilience/test_strategies.py +++ b/pySDC/tests/test_projects/test_resilience/test_strategies.py @@ -22,7 +22,7 @@ LOGGER_LEVEL = 30 -def single_test_vdp(strategy_name, useMPI, num_procs): +def single_test(strategy_name, useMPI, num_procs): import numpy as np from pySDC.helpers.stats_helper import get_sorted import pySDC.projects.Resilience.strategies as strategies @@ -103,18 +103,18 @@ def single_test_vdp(strategy_name, useMPI, num_procs): @pytest.mark.mpi4py @pytest.mark.parametrize('strategy_name', STRATEGY_NAMES + STRATEGY_NAMES_MPIONLY) -def test_strategy_with_vdp_MPI(strategy_name, num_procs=1): - single_test_vdp(strategy_name=strategy_name, useMPI=True, num_procs=num_procs) +def test_strategy_MPI(strategy_name, num_procs=1): + single_test(strategy_name=strategy_name, useMPI=True, num_procs=num_procs) @pytest.mark.base @pytest.mark.parametrize('strategy_name', STRATEGY_NAMES + STRATEGY_NAMES_NONMPIONLY) -def test_strategy_with_vdp_nonMPI(strategy_name, num_procs=1): - single_test_vdp(strategy_name=strategy_name, useMPI=False, num_procs=num_procs) +def test_strategy_nonMPI(strategy_name, num_procs=1): + single_test(strategy_name=strategy_name, useMPI=False, num_procs=num_procs) if __name__ == '__main__': for name in STRATEGY_NAMES + STRATEGY_NAMES_NONMPIONLY: - test_strategy_with_vdp_nonMPI(name) + test_strategy_nonMPI(name) for name in STRATEGY_NAMES: - test_strategy_with_vdp_MPI(name) + test_strategy_MPI(name)