diff --git a/pySDC/helpers/plot_helper.py b/pySDC/helpers/plot_helper.py index 51e8783b96..8ddea7cd0c 100644 --- a/pySDC/helpers/plot_helper.py +++ b/pySDC/helpers/plot_helper.py @@ -43,11 +43,13 @@ def figsize_by_journal(journal, scale, ratio): # pragma: no cover 'JSC_beamer': 426.79135, 'Springer_Numerical_Algorithms': 338.58778, 'JSC_thesis': 434.26027, + 'TUHH_thesis': 426.79135, } # store text height in points here, get this from LaTeX using \the\textheight textheights = { 'JSC_beamer': 214.43411, 'JSC_thesis': 635.5, + 'TUHH_thesis': 631.65118, } assert ( journal in textwidths.keys() diff --git a/pySDC/helpers/pysdc_helper.py b/pySDC/helpers/pysdc_helper.py index 4acd969932..666d3a4826 100644 --- a/pySDC/helpers/pysdc_helper.py +++ b/pySDC/helpers/pysdc_helper.py @@ -15,7 +15,7 @@ class FrozenClass(object): def __setattr__(self, key, value): """ - Function called when setting arttributes + Function called when setting attributes Args: key: the attribute @@ -35,7 +35,7 @@ def __getattr__(self, key): if key in self.attrs: return None else: - super().__getattr__(key) + super().__getattribute__(key) @classmethod def add_attr(cls, key, raise_error_if_exists=False): diff --git a/pySDC/implementations/convergence_controller_classes/step_size_limiter.py b/pySDC/implementations/convergence_controller_classes/step_size_limiter.py index cb090eebd3..c3c7ede145 100644 --- a/pySDC/implementations/convergence_controller_classes/step_size_limiter.py +++ b/pySDC/implementations/convergence_controller_classes/step_size_limiter.py @@ -146,7 +146,7 @@ def get_new_step_size(self, controller, S, **kwargs): S, ) L.status.dt_new = dt_new - elif abs(L.status.dt_new / L.params.dt - 1) < self.params.dt_rel_min_slope: + elif abs(L.status.dt_new / L.params.dt - 1) < self.params.dt_rel_min_slope and not S.status.restart: L.status.dt_new = L.params.dt self.log( f"Step size did not change sufficiently to warrant step size change, keeping {L.status.dt_new:.2e}", @@ -154,3 +154,67 @@ def get_new_step_size(self, controller, S, **kwargs): ) return None + + +class StepSizeRounding(ConvergenceController): + """ + Class to round step size when using adaptive step size selection. + """ + + def setup(self, controller, params, description, **kwargs): + """ + Define parameters here + + Args: + controller (pySDC.Controller): The controller + params (dict): The params passed for this specific convergence controller + description (dict): The description object used to instantiate the controller + + Returns: + (dict): The updated params dictionary + """ + defaults = { + "control_order": +93, + "digits": 1, + "fac": 5, + } + return {**defaults, **super().setup(controller, params, description, **kwargs)} + + @staticmethod + def _round_step_size(dt, fac, digits): + dt_rounded = None + exponent = np.log10(dt) // 1 + + dt_norm = dt / 10 ** (exponent - digits) + dt_norm_round = (dt_norm // fac) * fac + dt_rounded = dt_norm_round * 10 ** (exponent - digits) + return dt_rounded + + def get_new_step_size(self, controller, S, **kwargs): + """ + Enforce an upper and lower limit to the step size here. + Be aware that this is only tested when a new step size has been determined. That means if you set an initial + value for the step size outside of the limits, and you don't do any further step size control, that value will + go through. + Also, the final step is adjusted such that we reach Tend as best as possible, which might give step sizes below + the lower limit set here. + + Args: + controller (pySDC.Controller): The controller + S (pySDC.Step): The current step + + Returns: + None + """ + for L in S.levels: + if L.status.dt_new is not None: + dt_rounded = self._round_step_size(L.status.dt_new, self.params.fac, self.params.digits) + + if L.status.dt_new != dt_rounded: + self.log( + f"Step size rounded from {L.status.dt_new:.6e} to {dt_rounded:.6e}", + S, + ) + L.status.dt_new = dt_rounded + + return None diff --git a/pySDC/implementations/problem_classes/generic_spectral.py b/pySDC/implementations/problem_classes/generic_spectral.py index 04f88823fc..556ff46760 100644 --- a/pySDC/implementations/problem_classes/generic_spectral.py +++ b/pySDC/implementations/problem_classes/generic_spectral.py @@ -387,3 +387,25 @@ def compute_residual_DAE_MPI(self, stage=None): L.status.updated = False return None + + +def get_extrapolated_error_DAE(self, S, **kwargs): + """ + The extrapolation estimate combines values of u and f from multiple steps to extrapolate and compare to the + solution obtained by the time marching scheme. This function can be used in `EstimateExtrapolationError`. + + Args: + S (pySDC.Step): The current step + + Returns: + None + """ + u_ex = self.get_extrapolated_solution(S) + diff_mask = S.levels[0].prob.diff_mask + if u_ex is not None: + S.levels[0].status.error_extrapolation_estimate = ( + abs((u_ex - S.levels[0].u[-1])[diff_mask]) * self.coeff.prefactor + ) + # print([abs(me) for me in (u_ex - S.levels[0].u[-1]) * self.coeff.prefactor]) + else: + S.levels[0].status.error_extrapolation_estimate = None diff --git a/pySDC/implementations/sweeper_classes/Runge_Kutta.py b/pySDC/implementations/sweeper_classes/Runge_Kutta.py index 801beb9620..a0805f1f2f 100644 --- a/pySDC/implementations/sweeper_classes/Runge_Kutta.py +++ b/pySDC/implementations/sweeper_classes/Runge_Kutta.py @@ -441,9 +441,12 @@ def update_nodes(self): rhs += lvl.dt * (self.QI[m + 1, j] * lvl.f[j].impl + self.QE[m + 1, j] * lvl.f[j].expl) # implicit solve with prefactor stemming from the diagonal of Qd, use previous stage as initial guess - lvl.u[m + 1][:] = prob.solve_system( - rhs, lvl.dt * self.QI[m + 1, m + 1], lvl.u[m], lvl.time + lvl.dt * self.coll.nodes[m + 1] - ) + if self.QI[m + 1, m + 1] != 0: + lvl.u[m + 1][:] = prob.solve_system( + rhs, lvl.dt * self.QI[m + 1, m + 1], lvl.u[m], lvl.time + lvl.dt * self.coll.nodes[m + 1] + ) + else: + lvl.u[m + 1][:] = rhs[:] # update function values (we don't usually need to evaluate the RHS at the solution of the step) if m < M - self.coll.num_solution_stages or self.params.eval_rhs_at_right_boundary: diff --git a/pySDC/projects/Resilience/RBC.py b/pySDC/projects/Resilience/RBC.py new file mode 100644 index 0000000000..46366c5d56 --- /dev/null +++ b/pySDC/projects/Resilience/RBC.py @@ -0,0 +1,334 @@ +# script to run a Rayleigh-Benard Convection problem +from pySDC.implementations.problem_classes.generic_spectral import compute_residual_DAE, get_extrapolated_error_DAE +from pySDC.implementations.problem_classes.RayleighBenard import RayleighBenard +from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI +from pySDC.projects.Resilience.hook import hook_collection, LogData +from pySDC.projects.Resilience.strategies import merge_descriptions +from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order +from pySDC.core.convergence_controller import ConvergenceController +from pySDC.implementations.convergence_controller_classes.estimate_extrapolation_error import ( + EstimateExtrapolationErrorNonMPI, +) +from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence + +from pySDC.core.errors import ConvergenceError + +import numpy as np + + +def u_exact(self, t, u_init=None, t_init=None, recompute=False, _t0=None): + import pickle + import os + + path = f'data/stats/RBC-u_init-{t:.8f}.pickle' + if os.path.exists(path) and not recompute and t_init is None: + with open(path, 'rb') as file: + data = pickle.load(file) + else: + from pySDC.helpers.stats_helper import get_sorted + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + + convergence_controllers = { + Adaptivity: {'e_tol': 1e-8, 'dt_rel_min_slope': 0.25, 'dt_min': 1e-5}, + } + desc = {'convergence_controllers': convergence_controllers} + + u0 = self._u_exact(0) if u_init is None else u_init + t0 = 0 if t_init is None else t_init + + if t == t0: + return u0 + else: + u0 = u0 if _t0 is None else self.u_exact(_t0) + _t0 = t0 if _t0 is None else _t0 + + stats, _, _ = run_RBC(Tend=t, u0=u0, t0=_t0, custom_description=desc) + + u = get_sorted(stats, type='u', recomputed=False)[-1] + data = u[1] + + if t0 == 0: + with open(path, 'wb') as file: + pickle.dump(data, file) + + return data + + +RayleighBenard._u_exact = RayleighBenard.u_exact +RayleighBenard.u_exact = u_exact + +PROBLEM_PARAMS = {'Rayleigh': 2e4, 'nx': 256, 'nz': 128} + + +class ReachTendExactly(ConvergenceController): + """ + This convergence controller will adapt the step size of (hopefully) the last step such that `Tend` is reached very closely. + Please pass the same `Tend` that you pass to the controller to the params for this to work. + """ + + def setup(self, controller, params, description, **kwargs): + defaults = { + "control_order": +94, + "Tend": None, + 'min_step_size': 1e-10, + } + return {**defaults, **super().setup(controller, params, description, **kwargs)} + + def get_new_step_size(self, controller, step, **kwargs): + L = step.levels[0] + + if not CheckConvergence.check_convergence(step): + return None + + dt = L.status.dt_new if L.status.dt_new else L.params.dt + time_left = self.params.Tend - L.time - L.dt + if time_left <= dt + self.params.min_step_size and not step.status.restart and time_left > 0: + dt_new = ( + min([dt + self.params.min_step_size, max([time_left, self.params.min_step_size])]) + + np.finfo('float').eps * 10 + ) + if dt_new != L.status.dt_new: + L.status.dt_new = dt_new + self.log( + f'Reducing step size from {dt:12e} to {L.status.dt_new:.12e} because there is only {time_left:.12e} left.', + step, + ) + + +def run_RBC( + custom_description=None, + num_procs=1, + Tend=21.0, + hook_class=LogData, + fault_stuff=None, + custom_controller_params=None, + u0=None, + t0=20.0, + use_MPI=False, + **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 + 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 + """ + EstimateExtrapolationErrorNonMPI.get_extrapolated_error = get_extrapolated_error_DAE + + level_params = {} + level_params['dt'] = 1e-3 + level_params['restol'] = -1 + + sweeper_params = {} + sweeper_params['quad_type'] = 'RADAU-RIGHT' + sweeper_params['num_nodes'] = 3 + sweeper_params['QI'] = 'LU' + sweeper_params['QE'] = 'PIC' + + from mpi4py import MPI + + problem_params = {'comm': MPI.COMM_SELF, **PROBLEM_PARAMS} + + step_params = {} + step_params['maxiter'] = 5 + + convergence_controllers = {} + convergence_controllers[ReachTendExactly] = {'Tend': Tend} + from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeRounding + + convergence_controllers[StepSizeRounding] = {} + + 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} + + imex_1st_order.compute_residual = compute_residual_DAE + + description = {} + description['problem_class'] = RayleighBenard + description['problem_params'] = problem_params + description['sweeper_class'] = imex_1st_order + description['sweeper_params'] = sweeper_params + description['level_params'] = level_params + description['step_params'] = step_params + description['convergence_controllers'] = convergence_controllers + + if custom_description is not None: + description = merge_descriptions(description, custom_description) + + 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 + + t0 = 0.0 if t0 is None else t0 + 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 generate_data_for_fault_stats(): + prob = RayleighBenard(**PROBLEM_PARAMS) + _ts = np.linspace(0, 22, 220, dtype=float) + for i in range(len(_ts) - 1): + prob.u_exact(_ts[i + 1], _t0=_ts[i]) + + +def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recompute=False): + from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun + from pySDC.implementations.hooks.log_work import LogSDCIterations, LogWork + from pySDC.implementations.convergence_controller_classes.crash import StopAtNan + from pySDC.helpers.stats_helper import get_sorted + import pickle + import os + + sweeper_params = {'num_nodes': num_nodes} + level_params = {'e_tol': e_tol, 'restol': restol} + step_params = {'maxiter': 99} + convergence_controllers = {StopAtNan: {}} + + errors = [] + dts = [] + for i in range(steps): + dts += [dt / 2**i] + + path = f'data/stats/RBC-u_order-{t:.8f}-{dts[-1]:.8e}-{num_nodes}-{e_tol:.2e}-{restol:.2e}.pickle' + + if os.path.exists(path) and not recompute: + with open(path, 'rb') as file: + stats = pickle.load(file) + else: + + level_params['dt'] = dts[-1] + + desc = { + 'sweeper_params': sweeper_params, + 'level_params': level_params, + 'step_params': step_params, + 'convergence_controllers': convergence_controllers, + } + + stats, _, _ = run_RBC( + Tend=t + dt, + t0=t, + custom_description=desc, + hook_class=[LogGlobalErrorPostRun, LogSDCIterations, LogWork], + ) + + with open(path, 'wb') as file: + pickle.dump(stats, file) + + e = get_sorted(stats, type='e_global_post_run') + # k = get_sorted(stats, type='k') + + if len(e) > 0: + errors += [e[-1][1]] + else: + errors += [np.nan] + + errors = np.array(errors) + dts = np.array(dts) + mask = np.isfinite(errors) + max_error = np.nanmax(errors) + + errors = errors[mask] + dts = dts[mask] + ax.loglog(dts, errors, label=f'{num_nodes} nodes', marker='x') + ax.loglog( + dts, [max_error * (me / dts[0]) ** (2 * num_nodes - 1) for me in dts], ls='--', label=f'order {2*num_nodes-1}' + ) + ax.set_xlabel(r'$\Delta t$') + ax.set_ylabel(r'global error') + ax.legend(frameon=False) + + +def check_order(t=14, dt=1e-1, steps=6): + prob = RayleighBenard(**PROBLEM_PARAMS) + _ts = [0, t, t + dt] + for i in range(len(_ts) - 1): + prob.u_exact(_ts[i + 1], _t0=_ts[i]) + + import matplotlib.pyplot as plt + + fig, ax = plt.subplots(1, 1) + for num_nodes in [1, 2, 3, 4]: + plot_order(t=t, dt=dt, steps=steps, num_nodes=num_nodes, ax=ax, restol=1e-9) + ax.set_title(f't={t:.2f}, dt={dt:.2f}') + plt.show() + + +def plot_step_size(t0=0, Tend=30, e_tol=1e-3, recompute=False): + import matplotlib.pyplot as plt + import pickle + import os + from pySDC.helpers.stats_helper import get_sorted + + path = f'data/stats/RBC-u-{t0:.8f}-{Tend:.8f}-{e_tol:.2e}.pickle' + if os.path.exists(path) and not recompute: + with open(path, 'rb') as file: + stats = pickle.load(file) + else: + from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity + + convergence_controllers = { + Adaptivity: {'e_tol': e_tol, 'dt_rel_min_slope': 0.25, 'dt_min': 1e-5}, + } + desc = {'convergence_controllers': convergence_controllers} + + stats, _, _ = run_RBC(Tend=Tend, t0=t0, custom_description=desc) + + with open(path, 'wb') as file: + pickle.dump(stats, file) + + fig, ax = plt.subplots(1, 1) + + dt = get_sorted(stats, type='dt', recomputed=False) + ax.plot([me[0] for me in dt], [me[1] for me in dt]) + ax.set_ylabel(r'$\Delta t$') + ax.set_xlabel(r'$t$') + ax.set_yscale('log') + plt.show() + + +if __name__ == '__main__': + # plot_step_size(0, 3) + generate_data_for_fault_stats() + # check_order(t=20, dt=1., steps=7) + # stats, _, _ = run_RBC() diff --git a/pySDC/projects/Resilience/README.rst b/pySDC/projects/Resilience/README.rst index becbd2742f..3092c023b2 100644 --- a/pySDC/projects/Resilience/README.rst +++ b/pySDC/projects/Resilience/README.rst @@ -50,10 +50,21 @@ Then, navigate to this directory, `pySDC/projects/Resilience/` and run the follo .. code-block:: 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 + python paper_plots.py --target=adaptivity Possibly, you need to create some directories in this one to store and load things, if path errors occur. + +Reproduction of the plots in the resilience paper +------------------------------------------------- +To reproduce the plots you need to install pySDC using this project's `environment.yml` file, which is in the same directory as this README. + +.. code-block:: bash + + mpirun -np 4 python work_precision.py + mpirun -np 4 python fault_stats.py prob run_Lorenz + mpirun -np 4 python fault_stats.py prob run_Schroedinger + mpirun -np 4 python fault_stats.py prob run_AC + mpirun -np 4 python fault_stats.py prob run_RBC + python paper_plots.py --target=resilience + +Please be aware that generating the fault data for Rayleigh-Benard requires generating reference solutions, which may take several hours. diff --git a/pySDC/projects/Resilience/fault_stats.py b/pySDC/projects/Resilience/fault_stats.py index f28aa11746..6bb22c9827 100644 --- a/pySDC/projects/Resilience/fault_stats.py +++ b/pySDC/projects/Resilience/fault_stats.py @@ -21,6 +21,7 @@ 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.RBC import run_RBC from pySDC.projects.Resilience.strategies import BaseStrategy, AdaptivityStrategy, IterateStrategy, HotRodStrategy import logging @@ -632,6 +633,8 @@ def get_name(self, strategy=None, faults=True, mode=None): prob_name = 'Quench' elif self.prob.__name__ == 'run_AC': prob_name = 'Allen-Cahn' + elif self.prob.__name__ == 'run_RBC': + prob_name = 'Rayleigh-Benard' else: raise NotImplementedError(f'Name not implemented for problem {self.prob}') @@ -1565,6 +1568,10 @@ def parse_args(): kwargs['prob'] = run_Schroedinger elif sys.argv[i + 1] == 'run_quench': kwargs['prob'] = run_quench + elif sys.argv[i + 1] == 'run_AC': + kwargs['prob'] = run_AC + elif sys.argv[i + 1] == 'run_RBC': + kwargs['prob'] = run_RBC else: raise NotImplementedError elif 'num_procs' in sys.argv[i]: @@ -1654,11 +1661,11 @@ def compare_adaptivity_modes(): def main(): kwargs = { - 'prob': run_AC, + 'prob': run_RBC, 'num_procs': 1, 'mode': 'default', 'runs': 2000, - 'reload': True, + 'reload': False, **parse_args(), } diff --git a/pySDC/projects/Resilience/paper_plots.py b/pySDC/projects/Resilience/paper_plots.py index e3ad6d435e..04bea2d954 100644 --- a/pySDC/projects/Resilience/paper_plots.py +++ b/pySDC/projects/Resilience/paper_plots.py @@ -9,6 +9,7 @@ run_vdp, run_quench, run_AC, + run_RBC, RECOVERY_THRESH_ABS, ) from pySDC.projects.Resilience.strategies import ( @@ -187,21 +188,24 @@ def plot_recovery_rate_recoverable_only(stats_analyser, fig, ax, **kwargs): # p ) -def compare_recovery_rate_problems(**kwargs): # pragma: no cover +def compare_recovery_rate_problems(target='resilience', **kwargs): # pragma: no cover """ - Compare the recovery rate for vdP, Lorenz and Schroedinger problems. + Compare the recovery rate for various problems. Only faults that can be recovered are shown. Returns: None """ - stats = [ - get_stats(run_vdp, **kwargs), - get_stats(run_quench, **kwargs), - get_stats(run_Schroedinger, **kwargs), - get_stats(run_AC, **kwargs), - ] - titles = ['Van der Pol', 'Quench', r'Schr\"odinger', 'Allen-Cahn'] + if target == 'resilience': + problems = [run_Lorenz, run_Schroedinger, run_AC, run_RBC] + titles = ['Lorenz', r'Schr\"odinger', 'Allen-Cahn', 'Rayleigh-Benard'] + elif target == 'thesis': + problems = [run_vdp, run_Lorenz, run_AC, run_RBC] # TODO: swap in Gray-Scott + titles = ['Van der Pol', 'Lorenz', 'Allen-Cahn', 'Rayleigh-Benard'] + else: + raise NotImplementedError() + + stats = [get_stats(problem, **kwargs) for problem in problems] my_setup_mpl() fig, axs = plt.subplots(2, 2, figsize=figsize_by_journal(JOURNAL, 1, 0.8), sharey=True) @@ -420,6 +424,42 @@ def plot_quench_solution(): # pragma: no cover savefig(fig, 'quench_sol') +def plot_RBC_solution(): # pragma: no cover + """ + Plot solution of Rayleigh-Benard convection + """ + my_setup_mpl() + + from mpl_toolkits.axes_grid1 import make_axes_locatable + + plt.rcParams['figure.constrained_layout.use'] = True + fig, axs = plt.subplots(2, 1, sharex=True, sharey=True, figsize=figsize_by_journal(JOURNAL, 1.0, 0.45)) + caxs = [] + divider = make_axes_locatable(axs[0]) + caxs += [divider.append_axes('right', size='3%', pad=0.03)] + divider2 = make_axes_locatable(axs[1]) + caxs += [divider2.append_axes('right', size='3%', pad=0.03)] + + from pySDC.projects.Resilience.RBC import RayleighBenard, PROBLEM_PARAMS + + prob = RayleighBenard(**PROBLEM_PARAMS) + + def _plot(t, ax, cax): + u_hat = prob.u_exact(t) + u = prob.itransform(u_hat) + im = ax.pcolormesh(prob.X, prob.Z, u[prob.index('T')], rasterized=True, cmap='plasma') + fig.colorbar(im, cax, label=f'$T(t={{{t}}})$') + + _plot(0, axs[0], caxs[0]) + _plot(21, axs[1], caxs[1]) + + axs[1].set_xlabel('$x$') + axs[0].set_ylabel('$z$') + axs[1].set_ylabel('$z$') + + savefig(fig, 'RBC_sol', tight_layout=False) + + def plot_Schroedinger_solution(): # pragma: no cover from pySDC.implementations.problem_classes.NonlinearSchroedinger_MPIFFT import nonlinearschroedinger_imex @@ -596,10 +636,11 @@ def make_plots_for_adaptivity_paper(): # pragma: no cover def make_plots_for_resilience_paper(): # pragma: no cover + compare_recovery_rate_problems(target='resilience', num_procs=1, strategy_type='SDC') + plot_RBC_solution() plot_recovery_rate(get_stats(run_vdp)) plot_fault_vdp(0) plot_fault_vdp(13) - compare_recovery_rate_problems(num_procs=1, strategy_type='SDC') def make_plots_for_notes(): # pragma: no cover @@ -614,8 +655,37 @@ def make_plots_for_notes(): # pragma: no cover analyse_resilience(run_quench, format='png') +def make_plots_for_thesis(): # pragma: no cover + global JOURNAL + JOURNAL = 'TUHH_thesis' + + plot_RBC_solution() + # plot_vdp_solution() + + # plot_adaptivity_stuff() + compare_recovery_rate_problems(target='thesis', num_procs=1, strategy_type='SDC') + + if __name__ == "__main__": - # make_plots_for_notes() - # make_plots_for_SIAM_CSE23() - # make_plots_for_TIME_X_website() - make_plots_for_adaptivity_paper() + import argparse + + parser = argparse.ArgumentParser() + parser.add_argument( + '--target', choices=['adaptivity', 'resilience', 'thesis', 'notes', 'SIAM_CSE23', 'TIME_X_website'], type=str + ) + args = parser.parse_args() + + if args.target == 'adaptivity': + make_plots_for_adaptivity_paper() + elif args.target == 'resilience': + make_plots_for_resilience_paper() + elif args.target == 'thesis': + make_plots_for_thesis() + elif args.target == 'notes': + make_plots_for_notes() + elif args.target == 'SIAM_CSE23': + make_plots_for_SIAM_CSE23() + elif args.target == 'TIME_X_website': + make_plots_for_TIME_X_website() + else: + raise NotImplementedError(f'Don\'t know how to make plots for target {args.target}') diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py index 2d4edaa689..09bcbff62d 100644 --- a/pySDC/projects/Resilience/strategies.py +++ b/pySDC/projects/Resilience/strategies.py @@ -130,6 +130,8 @@ def get_fault_args(self, problem, num_procs): args['time'] = 0.3 elif problem.__name__ == "run_AC": args['time'] = 1e-2 + elif problem.__name__ == "run_RBC": + args['time'] = 20.19 return args @@ -150,13 +152,16 @@ 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', 'run_AC']: + if problem.__name__ in ['run_Schroedinger', 'run_quench', 'run_AC', 'run_RBC']: rnd_params['min_node'] = 1 if problem.__name__ == "run_quench": rnd_params['iteration'] = 5 elif problem.__name__ == 'run_Lorenz': rnd_params['iteration'] = 5 + elif problem.__name__ == 'run_RBC': + rnd_params['problem_pos'] = [3, 16, 16] + return rnd_params @property @@ -206,6 +211,8 @@ def get_Tend(cls, problem, num_procs=1): return 500.0 elif problem.__name__ == "run_AC": return 0.025 + elif problem.__name__ == "run_RBC": + return 21 else: raise NotImplementedError('I don\'t have a final time for your problem!') @@ -269,6 +276,9 @@ def get_base_parameters(self, problem, num_procs=1): 'nu': 2, } custom_description['level_params'] = {'restol': -1, 'dt': 0.1 * eps**2} + elif problem.__name__ == 'run_RBC': + custom_description['level_params']['dt'] = 5e-2 + custom_description['step_params'] = {'maxiter': 5} custom_description['convergence_controllers'] = { # StepSizeLimiter: {'dt_min': self.get_Tend(problem=problem, num_procs=num_procs) / self.max_steps} @@ -287,6 +297,7 @@ def get_base_parameters(self, problem, num_procs=1): 'run_Schroedinger': 150, 'run_quench': 150, 'run_AC': 150, + 'run_RBC': 1000, } custom_description['convergence_controllers'][StopAtMaxRuntime] = { @@ -373,7 +384,7 @@ def get_custom_description(self, problem, num_procs=1): 'maxiter': 15, } - if self.newton_inexactness and problem.__name__ not in ['run_Schroedinger', 'run_AC']: + if self.newton_inexactness and problem.__name__ not in ['run_Schroedinger', 'run_AC', 'run_RBC']: if problem.__name__ == 'run_quench': inexactness_params['ratio'] = 1e-1 inexactness_params['min_tol'] = 1e-11 @@ -521,6 +532,7 @@ def get_custom_description(self, problem, num_procs): dt_max = np.inf dt_slope_max = np.inf + dt_slope_min = 0 if problem.__name__ == "run_piline": e_tol = 1e-7 @@ -545,6 +557,9 @@ def get_custom_description(self, problem, num_procs): elif problem.__name__ == "run_AC": e_tol = 1e-7 # dt_max = 0.1 * base_params['problem_params']['eps'] ** 2 + elif problem.__name__ == 'run_RBC': + e_tol = 1e-4 + dt_slope_min = 1 else: raise NotImplementedError( @@ -555,6 +570,7 @@ def get_custom_description(self, problem, num_procs): custom_description['convergence_controllers'][Adaptivity] = { 'e_tol': e_tol, 'dt_slope_max': dt_slope_max, + 'dt_rel_min_slope': dt_slope_min, } custom_description['convergence_controllers'][StepSizeLimiter] = { 'dt_max': dt_max, @@ -756,6 +772,8 @@ def get_custom_description(self, problem, num_procs): restol = 1e-7 elif problem.__name__ == "run_AC": restol = 1e-11 + elif problem.__name__ == "run_RBC": + restol = 1e-4 else: raise NotImplementedError( 'I don\'t have a residual tolerance for your problem. Please add one to the \ @@ -827,6 +845,9 @@ def get_custom_description(self, problem, num_procs, *args, **kwargs): elif problem.__name__ == "run_AC": desc['level_params']['restol'] = 1e-11 desc['level_params']['dt'] = 0.4 * desc['problem_params']['eps'] ** 2 / 8.0 + elif problem.__name__ == "run_RBC": + desc['level_params']['dt'] = 7e-2 + desc['level_params']['restol'] = 1e-9 return desc def get_custom_description_for_faults(self, problem, *args, **kwargs): @@ -835,6 +856,8 @@ def get_custom_description_for_faults(self, problem, *args, **kwargs): desc['level_params']['dt'] = 5.0 elif problem.__name__ == 'run_AC': desc['level_params']['dt'] = 0.6 * desc['problem_params']['eps'] ** 2 + elif problem.__name__ == 'run_RBC': + desc['level_params']['restol'] = 1e-6 return desc def get_reference_value(self, problem, key, op, num_procs=1): @@ -928,6 +951,9 @@ def get_custom_description(self, problem, num_procs): elif problem.__name__ == 'run_AC': HotRod_tol = 9.564437e-06 maxiter = 6 + elif problem.__name__ == 'run_RBC': + HotRod_tol = 6.34e-6 + maxiter = 6 else: raise NotImplementedError( 'I don\'t have a tolerance for Hot Rod for your problem. Please add one to the\ @@ -1302,6 +1328,7 @@ def get_custom_description(self, problem, num_procs): The custom descriptions you can supply to the problem when running it ''' from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityRK, Adaptivity + from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeSlopeLimiter from pySDC.implementations.convergence_controller_classes.basic_restarting import BasicRestarting from pySDC.implementations.sweeper_classes.Runge_Kutta import ARK548L2SA @@ -1323,6 +1350,9 @@ def get_custom_description(self, problem, num_procs): }, } + if problem.__name__ == "run_RBC": + rk_params['convergence_controllers'][StepSizeSlopeLimiter] = {'dt_rel_min_slope': 0.25} + custom_description = merge_descriptions(adaptivity_description, rk_params) return custom_description @@ -1857,6 +1887,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 + from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence base_params = super().get_custom_description(problem, num_procs) custom_description = {} @@ -1864,7 +1895,10 @@ def get_custom_description(self, problem, num_procs): dt_max = np.inf restol_rel = 1e-4 restol_min = 1e-12 + restol_max = 1e-5 + dt_slope_min = 0 dt_min = 0 + abort_at_growing_residual = True level_params = {} problem_params = {} @@ -1890,6 +1924,14 @@ def get_custom_description(self, problem, num_procs): e_tol = 1.0e-4 restol_rel = 1e-3 # dt_max = 0.1 * base_params['problem_params']['eps'] ** 2 + elif problem.__name__ == "run_RBC": + e_tol = 5e-3 + dt_slope_min = 1.0 + abort_at_growing_residual = False + restol_rel = 1e-3 + restol_max = 1e-1 + restol_min = 5e-7 + self.max_slope = 4 else: raise NotImplementedError( 'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\ @@ -1901,14 +1943,17 @@ def get_custom_description(self, problem, num_procs): 'e_tol': e_tol, 'restol_rel': restol_rel if self.use_restol_rel else 1e-11, 'restol_min': restol_min if self.use_restol_rel else 1e-12, + 'restol_max': restol_max if self.use_restol_rel else 1e-5, 'restart_at_maxiter': True, 'factor_if_not_converged': self.max_slope, 'interpolate_between_restarts': self.interpolate_between_restarts, + 'abort_at_growing_residual': abort_at_growing_residual, }, StepSizeLimiter: { 'dt_max': dt_max, 'dt_slope_max': self.max_slope, 'dt_min': dt_min, + 'dt_rel_min_slope': dt_slope_min, }, } custom_description['level_params'] = level_params diff --git a/pySDC/projects/Resilience/tests/test_RBC.py b/pySDC/projects/Resilience/tests/test_RBC.py new file mode 100644 index 0000000000..71dce01bb8 --- /dev/null +++ b/pySDC/projects/Resilience/tests/test_RBC.py @@ -0,0 +1,69 @@ +import pytest + + +@pytest.mark.mpi4py +@pytest.mark.parametrize('Tend', [1.2345, 1, 4.0 / 3.0]) +@pytest.mark.parametrize('dt', [0.1, 1.0 / 3.0, 0.999999]) +def test_ReachTendExactly(Tend, dt, num_procs=1): + from pySDC.projects.Resilience.RBC import ReachTendExactly + from pySDC.implementations.hooks.log_solution import LogSolution + from pySDC.implementations.problem_classes.Van_der_Pol_implicit import vanderpol + from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit + from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + from pySDC.helpers.stats_helper import get_sorted + import numpy as np + + level_params = {} + level_params['dt'] = dt + + sweeper_params = {} + sweeper_params['quad_type'] = 'RADAU-RIGHT' + sweeper_params['num_nodes'] = 2 + sweeper_params['QI'] = 'LU' + + problem_params = { + 'mu': 0.0, + 'newton_tol': 1e-1, + 'newton_maxiter': 9, + 'u0': np.array([2.0, 0.0]), + 'relative_tolerance': True, + } + + step_params = {} + step_params['maxiter'] = 2 + + convergence_controllers = {ReachTendExactly: {'Tend': Tend}} + + controller_params = {} + controller_params['logger_level'] = 15 + controller_params['hook_class'] = LogSolution + controller_params['mssdc_jac'] = False + + description = {} + description['problem_class'] = vanderpol + 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 + description['convergence_controllers'] = convergence_controllers + + t0 = 0.0 + + controller = controller_nonMPI(num_procs=num_procs, controller_params=controller_params, description=description) + + P = controller.MS[0].levels[0].prob + uinit = P.u_exact(t0) + + uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) + + u = get_sorted(stats, type='u') + t_last = u[-1][0] + + assert np.isclose( + t_last, Tend, atol=1e-10 + ), f'Expected {Tend=}, but got {t_last}, which is off by {t_last - Tend:.8e}' + + +if __name__ == '__main__': + test_ReachTendExactly(1.2345, 0.1) diff --git a/pySDC/projects/Resilience/work_precision.py b/pySDC/projects/Resilience/work_precision.py index 0a30ce47e0..86f8600ed8 100644 --- a/pySDC/projects/Resilience/work_precision.py +++ b/pySDC/projects/Resilience/work_precision.py @@ -12,6 +12,7 @@ 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.RBC import run_RBC from pySDC.helpers.stats_helper import get_sorted, filter_stats from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal @@ -38,6 +39,7 @@ def std_log(x): 'k_linear': ('work_linear', sum, None), 'k_Newton_no_restart': ('work_newton', sum, False), 'k_rhs': ('work_rhs', sum, None), + 'k_factorizations': ('work_factorizations', sum, None), 'num_steps': ('dt', len, None), 'restart': ('restart', sum, None), 'dt_mean': ('dt', np.mean, False), @@ -70,6 +72,14 @@ def get_forbidden_combinations(problem, strategy, **kwargs): return False +Tends = { + 'run_RBC': 16, +} +t0s = { + 'run_RBC': 0.2, +} + + def single_run( problem, strategy, @@ -132,12 +142,16 @@ def single_run( } problem_args = {} if problem_args is None else problem_args + Tend = Tends.get(problem.__name__, None) if Tend is None else Tend + t0 = t0s.get(problem.__name__, None) + stats, controller, crash = problem( custom_description=description, Tend=strategy.get_Tend(problem, num_procs) if Tend is None else Tend, hook_class=[LogData, LogWork, LogGlobalErrorPostRun] + hooks, custom_controller_params=controller_params, use_MPI=True, + t0=t0, comm=comm_time, **problem_args, ) @@ -274,6 +288,8 @@ def record_work_precision( exponents = [0, 1, 2, 3, 5][::-1] else: exponents = [-3, -2, -1, 0, 0.2, 0.8, 1][::-1] + if problem.__name__ == 'run_RBC': + exponents = [1, 0, -1, -2, -3, -4] elif param == 'dt': power = 2.0 exponents = [-1, 0, 1, 2, 3][::-1] @@ -294,6 +310,9 @@ def record_work_precision( param_range = [1e-5, 1e-6, 1e-7, 1e-8, 1e-9] elif param == 'dt': param_range = [1.25, 2.5, 5.0, 10.0, 20.0][::-1] + if problem.__name__ == 'run_RBC': + if param == 'dt': + param_range = [2e-1, 1e-1, 8e-2, 6e-2] # run multiple times with different parameters for i in range(len(param_range)): @@ -538,6 +557,7 @@ def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only 'k_Newton': 'Newton iterations', 'k_Newton_no_restart': 'Newton iterations (restarts excluded)', 'k_rhs': 'right hand side evaluations', + 'k_factorizations': 'matrix factorizations', 't': 'wall clock time / s', 'e_global': 'global error', 'e_global_rel': 'relative global error', @@ -785,14 +805,14 @@ def get_configs(mode, problem): AdaptivityPolynomialError, ) - if problem.__name__ in ['run_Schroedinger', 'run_AC']: + if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_RBC']: from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper else: from pySDC.implementations.sweeper_classes.generic_implicit_MPI import ( generic_implicit_MPI as parallel_sweeper, ) - newton_inexactness = False if problem.__name__ in ['run_vdp'] else True + newton_inexactness = False if problem.__name__ in ['run_vdp', 'run_RBC'] else True desc = {} desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': "EE"} @@ -811,7 +831,7 @@ def get_configs(mode, problem): RK_strategies = [] if problem.__name__ in ['run_Lorenz']: RK_strategies.append(ERKStrategy(useMPI=True)) - if problem.__name__ in ['run_Schroedinger', 'run_AC']: + if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_RBC']: RK_strategies.append(ARKStrategy(useMPI=True)) else: RK_strategies.append(ESDIRKStrategy(useMPI=True)) @@ -1591,12 +1611,12 @@ def aggregate_parallel_efficiency_plot(): # pragma: no cover ]: params = { 'mode': mode, - 'runs': 5, + 'runs': 1, 'plotting': comm_world.rank == 0, } params_single = { **params, - 'problem': run_AC, + 'problem': run_RBC, } single_problem(**params_single, work_key='t', precision_key='e_global_rel', record=record) diff --git a/pySDC/tests/test_convergence_controllers/test_step_size_limiter.py b/pySDC/tests/test_convergence_controllers/test_step_size_limiter.py index 3c988976aa..618fbea66e 100644 --- a/pySDC/tests/test_convergence_controllers/test_step_size_limiter.py +++ b/pySDC/tests/test_convergence_controllers/test_step_size_limiter.py @@ -108,5 +108,17 @@ def test_step_size_limiter(): assert L.status.dt_new == 0.5 +@pytest.mark.base +@pytest.mark.parametrize('dt', [1 / 3, 2 / 30]) +def test_step_size_rounding(dt): + from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeRounding + + expect = { + 1 / 3: 0.3, + 2 / 30: 0.065, + } + assert StepSizeRounding._round_step_size(dt, 5, 1) == expect[dt] + + if __name__ == '__main__': test_step_size_slope_limiter()