diff --git a/pySDC/projects/Resilience/GS.py b/pySDC/projects/Resilience/GS.py
new file mode 100644
index 0000000000000000000000000000000000000000..5fdc47d5414656ea06dc7d06b0d44c13bc2045db
--- /dev/null
+++ b/pySDC/projects/Resilience/GS.py
@@ -0,0 +1,176 @@
+# script to run a Gray-Scott problem
+from pySDC.implementations.problem_classes.GrayScott_MPIFFT import grayscott_imex_diffusion
+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.projects.Resilience.sweepers import imex_1st_order_efficient
+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.projects.Resilience.reachTendExactly import ReachTendExactly
+
+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/GS-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},
+ }
+ 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_GS(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
+
+
+if not hasattr(grayscott_imex_diffusion, '_u_exact'):
+ grayscott_imex_diffusion._u_exact = grayscott_imex_diffusion.u_exact
+ grayscott_imex_diffusion.u_exact = u_exact
+
+
+def run_GS(
+ custom_description=None,
+ num_procs=1,
+ Tend=1e2,
+ hook_class=LogData,
+ fault_stuff=None,
+ custom_controller_params=None,
+ u0=None,
+ t0=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
+ """
+
+ level_params = {}
+ level_params['dt'] = 1e0
+ level_params['restol'] = -1
+
+ sweeper_params = {}
+ sweeper_params['quad_type'] = 'RADAU-RIGHT'
+ sweeper_params['num_nodes'] = 3
+ sweeper_params['QI'] = 'MIN-SR-S'
+ sweeper_params['QE'] = 'PIC'
+ sweeper_params['initial_guess'] = 'copy'
+
+ from mpi4py import MPI
+
+ problem_params = {
+ 'comm': MPI.COMM_SELF,
+ 'num_blobs': -48,
+ 'L': 2,
+ 'nvars': (128,) * 2,
+ 'A': 0.062,
+ 'B': 0.1229,
+ 'Du': 2e-5,
+ 'Dv': 1e-5,
+ }
+
+ step_params = {}
+ step_params['maxiter'] = 5
+
+ convergence_controllers = {}
+ convergence_controllers[ReachTendExactly] = {'Tend': Tend}
+
+ controller_params = {}
+ controller_params['logger_level'] = 15
+ 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}
+
+ description = {}
+ description['problem_class'] = grayscott_imex_diffusion
+ description['problem_params'] = problem_params
+ description['sweeper_class'] = imex_1st_order_efficient
+ 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
+
+
+if __name__ == '__main__':
+ stats, _, _ = run_GS()
diff --git a/pySDC/projects/Resilience/Lorenz.py b/pySDC/projects/Resilience/Lorenz.py
index 817ad46be48ad2c9ff8317bc8f4631340df7f22a..c9cf9ebe3034a857fdef30f4f17da283a2b47cf6 100644
--- a/pySDC/projects/Resilience/Lorenz.py
+++ b/pySDC/projects/Resilience/Lorenz.py
@@ -49,6 +49,7 @@ def run_Lorenz(
sweeper_params['quad_type'] = 'RADAU-RIGHT'
sweeper_params['num_nodes'] = 3
sweeper_params['QI'] = 'IE'
+ sweeper_params['initial_guess'] = 'copy'
problem_params = {
'newton_tol': 1e-9,
diff --git a/pySDC/projects/Resilience/RBC.py b/pySDC/projects/Resilience/RBC.py
index 46366c5d562d44c360f07a18380c4827308135b8..94b898cdff2ecd64cc0628c0c4685e9d8e5af679 100644
--- a/pySDC/projects/Resilience/RBC.py
+++ b/pySDC/projects/Resilience/RBC.py
@@ -4,17 +4,18 @@ 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.projects.Resilience.sweepers import imex_1st_order_efficient
from pySDC.implementations.convergence_controller_classes.estimate_extrapolation_error import (
EstimateExtrapolationErrorNonMPI,
)
-from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence
+from pySDC.projects.Resilience.reachTendExactly import ReachTendExactly
from pySDC.core.errors import ConvergenceError
import numpy as np
+PROBLEM_PARAMS = {'Rayleigh': 3.2e5, 'nx': 256, 'nz': 128, 'max_cached_factorizations': 30}
+
def u_exact(self, t, u_init=None, t_init=None, recompute=False, _t0=None):
import pickle
@@ -54,45 +55,9 @@ def u_exact(self, t, u_init=None, t_init=None, recompute=False, _t0=None):
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,
- )
+if not hasattr(RayleighBenard, '_u_exact'):
+ RayleighBenard._u_exact = RayleighBenard.u_exact
+ RayleighBenard.u_exact = u_exact
def run_RBC(
@@ -105,6 +70,7 @@ def run_RBC(
u0=None,
t0=20.0,
use_MPI=False,
+ step_size_rounding=False,
**kwargs,
):
"""
@@ -135,6 +101,7 @@ def run_RBC(
sweeper_params['num_nodes'] = 3
sweeper_params['QI'] = 'LU'
sweeper_params['QE'] = 'PIC'
+ sweeper_params['initial_guess'] = 'copy'
from mpi4py import MPI
@@ -147,22 +114,23 @@ def run_RBC(
convergence_controllers[ReachTendExactly] = {'Tend': Tend}
from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeRounding
- convergence_controllers[StepSizeRounding] = {}
+ if step_size_rounding:
+ convergence_controllers[StepSizeRounding] = {}
controller_params = {}
- controller_params['logger_level'] = 30
+ controller_params['logger_level'] = 15
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
+ imex_1st_order_efficient.compute_residual = compute_residual_DAE
description = {}
description['problem_class'] = RayleighBenard
description['problem_params'] = problem_params
- description['sweeper_class'] = imex_1st_order
+ description['sweeper_class'] = imex_1st_order_efficient
description['sweeper_params'] = sweeper_params
description['level_params'] = level_params
description['step_params'] = step_params
@@ -204,14 +172,15 @@ def run_RBC(
return stats, controller, crash
-def generate_data_for_fault_stats():
+def generate_data_for_fault_stats(Tend):
prob = RayleighBenard(**PROBLEM_PARAMS)
- _ts = np.linspace(0, 22, 220, dtype=float)
+ _ts = np.linspace(0, Tend, Tend * 10 + 1, dtype=float)
for i in range(len(_ts) - 1):
- prob.u_exact(_ts[i + 1], _t0=_ts[i])
+ print(f'Generating reference solution from {_ts[i]:.4e} to {_ts[i+1]:.4e}')
+ prob.u_exact(_ts[i + 1], _t0=_ts[i], recompute=False)
-def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recompute=False):
+def plot_order(t, dt, steps, num_nodes, e_tol=1e-9, restol=1e-9, ax=None, recompute=False): # pragma: no cover
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
@@ -294,7 +263,7 @@ def check_order(t=14, dt=1e-1, steps=6):
plt.show()
-def plot_step_size(t0=0, Tend=30, e_tol=1e-3, recompute=False):
+def plot_step_size(t0=0, Tend=30, e_tol=1e-3, recompute=False): # pragma: no cover
import matplotlib.pyplot as plt
import pickle
import os
@@ -327,8 +296,86 @@ def plot_step_size(t0=0, Tend=30, e_tol=1e-3, recompute=False):
plt.show()
+def plot_factorizations_over_time(t0=0, Tend=50, e_tol=1e-3, recompute=False, adaptivity_mode='dt'): # pragma: no cover
+ import matplotlib.pyplot as plt
+ import pickle
+ import os
+ from pySDC.helpers.stats_helper import get_sorted
+ from pySDC.implementations.hooks.log_work import LogWork
+ from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity, AdaptivityPolynomialError
+ from pySDC.helpers.plot_helper import figsize_by_journal, setup_mpl
+ from pySDC.projects.Resilience.paper_plots import savefig
+
+ setup_mpl()
+
+ fig, axs = plt.subplots(1, 2, figsize=figsize_by_journal('TUHH_thesis', 1.0, 0.4))
+
+ if adaptivity_mode == 'dt':
+ adaptivity = Adaptivity
+ elif adaptivity_mode == 'dt_k':
+ adaptivity = AdaptivityPolynomialError
+
+ dt_controllers = {
+ 'basic': {
+ adaptivity: {
+ 'e_tol': e_tol,
+ }
+ },
+ 'min slope': {adaptivity: {'e_tol': e_tol, 'beta': 0.5, 'dt_rel_min_slope': 2}},
+ 'fixed': {},
+ # 'rounding': {adaptivity: {'e_tol': e_tol, 'beta': 0.5, 'dt_rel_min_slope': 2}, StepSizeRounding: {}},
+ }
+
+ for name, params in dt_controllers.items():
+ if adaptivity_mode == 'dt':
+ path = f'data/stats/RBC-u-{t0:.8f}-{Tend:.8f}-{e_tol:.2e}-{name}.pickle'
+ elif adaptivity_mode == 'dt_k':
+ path = f'data/stats/RBC-u-{t0:.8f}-{Tend:.8f}-{e_tol:.2e}-{name}-dtk.pickle'
+
+ if os.path.exists(path) and not recompute:
+ with open(path, 'rb') as file:
+ stats = pickle.load(file)
+ else:
+
+ convergence_controllers = {
+ **params,
+ }
+ desc = {'convergence_controllers': convergence_controllers}
+
+ if name == 'fixed':
+ if adaptivity_mode == 'dt':
+ desc['level_params'] = {'dt': 2e-2}
+ elif adaptivity_mode == 'dt_k':
+ desc['level_params'] = {'dt': 2e-3}
+ elif adaptivity_mode == 'dt_k':
+ desc['level_params'] = {'restol': 1e-7}
+
+ stats, _, _ = run_RBC(
+ Tend=Tend, t0=t0, custom_description=desc, hook_class=LogWork, step_size_rounding=False
+ )
+
+ with open(path, 'wb') as file:
+ pickle.dump(stats, file)
+
+ factorizations = get_sorted(stats, type='work_factorizations')
+ rhs_evals = get_sorted(stats, type='work_rhs')
+ axs[0].plot([me[0] for me in factorizations], np.cumsum([me[1] for me in factorizations]), label=name)
+ axs[1].plot([me[0] for me in rhs_evals], np.cumsum([me[1] for me in rhs_evals]), label=name)
+
+ axs[0].set_ylabel(r'matrix factorizations')
+ axs[1].set_ylabel(r'right hand side evaluations')
+ axs[0].set_xlabel(r'$t$')
+ axs[1].set_xlabel(r'$t$')
+ axs[0].set_yscale('log')
+ axs[1].set_yscale('log')
+ axs[0].legend(frameon=False)
+ savefig(fig, f'RBC_step_size_controller_{adaptivity_mode}')
+
+
if __name__ == '__main__':
- # plot_step_size(0, 3)
- generate_data_for_fault_stats()
+ # plot_step_size(0, 30)
+ generate_data_for_fault_stats(Tend=30)
+ # plot_factorizations_over_time(e_tol=1e-3, adaptivity_mode='dt')
+ # plot_factorizations_over_time(recompute=False, e_tol=1e-5, adaptivity_mode='dt_k')
# 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 3092c023b21dc46831ad3ee6159aad5093341b0d..75a5e44083a3878320dfeaa6ecc2ee6d4a79cdc1 100644
--- a/pySDC/projects/Resilience/README.rst
+++ b/pySDC/projects/Resilience/README.rst
@@ -43,13 +43,26 @@ You can see the results below, except for the solution, which looks the same as
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:
-
+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.
+Then, run the following commands in the same directory:
.. code-block:: bash
- mpirun -np 4 python work_precision.py
+ mpirun -np 4 python work_precision.py --mode=compare_strategies --problem=vdp
+ mpirun -np 4 python work_precision.py --mode=compare_strategies --problem=quench
+ mpirun -np 4 python work_precision.py --mode=compare_strategies --problem=Schroedinger
+ mpirun -np 4 python work_precision.py --mode=compare_strategies --problem=AC
+
+ mpirun -np 12 python work_precision.py --mode=parallel_efficiency --problem=vdp
+ mpirun -np 12 python work_precision.py --mode=parallel_efficiency --problem=quench
+ mpirun -np 12 python work_precision.py --mode=parallel_efficiency --problem=Schroedinger
+ mpirun -np 12 python work_precision.py --mode=parallel_efficiency --problem=AC
+
+ mpirun -np 4 python work_precision.py --mode=RK_comp --problem=vdp
+ mpirun -np 4 python work_precision.py --mode=RK_comp --problem=quench
+ mpirun -np 4 python work_precision.py --mode=RK_comp --problem=Schroedinger
+ mpirun -np 4 python work_precision.py --mode=RK_comp --problem=AC
+
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.
@@ -57,10 +70,10 @@ Possibly, you need to create some directories in this one to store and load thin
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.
+Then, run the following commands in the same directory:
.. 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
@@ -68,3 +81,34 @@ To reproduce the plots you need to install pySDC using this project's `environme
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.
+
+Reproduction of the plots in Thomas Baumann's thesis
+----------------------------------------------------
+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.
+
+To run the resilience experiments, run
+
+.. code-block:: bash
+
+ mpirun -np 4 python fault_stats.py prob run_Lorenz
+ mpirun -np 4 python fault_stats.py prob run_vdp
+ mpirun -np 4 python fault_stats.py prob run_GS
+ mpirun -np 4 python fault_stats.py prob run_RBC
+
+Please be aware that generating the fault data for Rayleigh-Benard requires generating reference solutions, which may take several hours.
+To run the experiments on computational efficiency, run for the choices of problems `vdp`, `Lorenz``, `GS`, `RBC` the following commands
+
+.. code-block:: bash
+
+ mpirun -np 4 python work_precision.py --mode=compare_strategies --problem=<problem>
+ mpirun -np 12 python work_precision.py --mode=parallel_efficiency_dt_k --problem=<problem>
+ mpirun -np 12 python work_precision.py --mode=parallel_efficiency_dt --problem=<problem>
+ mpirun -np 12 python work_precision.py --mode=RK_comp --problem=<problem>
+
+After you ran all experiments, you can generate the plots with the following command:
+
+.. code-block:: bash
+
+ python paper_plots.py --target=thesis
+
+As this plotting script also runs a few additional smaller experiments, do not be surprised if it takes a while.
diff --git a/pySDC/projects/Resilience/fault_stats.py b/pySDC/projects/Resilience/fault_stats.py
index b75e47ded873f9e3dbdd4f552f3e6ba8d151f5f4..2acdfb2243da58b7a10b508d5c740b0d97a89c3a 100644
--- a/pySDC/projects/Resilience/fault_stats.py
+++ b/pySDC/projects/Resilience/fault_stats.py
@@ -22,6 +22,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.GS import run_GS
from pySDC.projects.Resilience.strategies import BaseStrategy, AdaptivityStrategy, IterateStrategy, HotRodStrategy
import logging
@@ -108,7 +109,7 @@ class FaultStats:
Returns:
float: Tend to put into the run
'''
- return self.strategies[0].get_Tend(self.prob, self.num_procs)
+ return self.strategies[0].get_Tend(self.prob, self.num_procs, resilience_experiment=True)
def run_stats_generation(
self, runs=1000, step=None, comm=None, kwargs_range=None, faults=None, _reload=False, _runs_partial=0
@@ -637,6 +638,8 @@ class FaultStats:
prob_name = 'Allen-Cahn'
elif self.prob.__name__ == 'run_RBC':
prob_name = 'Rayleigh-Benard'
+ elif self.prob.__name__ == 'run_GS':
+ prob_name = 'Gray-Scott'
else:
raise NotImplementedError(f'Name not implemented for problem {self.prob}')
@@ -717,7 +720,7 @@ class FaultStats:
else:
try:
with_faults = self.load(faults=True, **kwargs)
- with_faults['recovered'] = with_faults['error'] < self.get_thresh(kwargs['strategy'])
+ with_faults['recovered'] = with_faults['error'] <= self.get_thresh(kwargs['strategy'])
self.store(faults=True, dat=with_faults, **kwargs)
except KeyError as error:
print(
@@ -836,10 +839,10 @@ class FaultStats:
me_recovered = np.zeros_like(me)
for i in range(len(me)):
- _mask = (dat[thingB] == admissable_thingB[i]) & mask
+ _mask = np.logical_and((dat[thingB] == admissable_thingB[i]), mask)
if _mask.any():
me[i] = op(dat, no_faults, thingA, _mask)
- me_recovered[i] = op(dat, no_faults, thingA, _mask & dat['recovered'])
+ me_recovered[i] = op(dat, no_faults, thingA, np.logical_and(_mask, dat['recovered']))
if recovered:
ax.plot(
@@ -965,7 +968,7 @@ class FaultStats:
if recoverable_only:
recoverable_mask = self.get_fixable_faults_only(strategy)
else:
- recoverable_mask = self.get_mask()
+ recoverable_mask = self.get_mask(strategy=strategy)
for thresh_idx in range(len(thresh_range)):
rec_mask = self.get_mask(
@@ -1514,7 +1517,7 @@ class FaultStats:
HR_strategy = HotRodStrategy(useMPI=self.use_MPI)
- description = HR_strategy.get_custom_description(self.prob, self.num_procs)
+ description = HR_strategy.get_custom_description_for_faults(self.prob, self.num_procs)
description['convergence_controllers'][HotRod]['HotRod_tol'] = 1e2
stats, _, _ = self.single_run(HR_strategy, force_params=description)
@@ -1578,6 +1581,8 @@ def parse_args():
kwargs['prob'] = run_AC
elif sys.argv[i + 1] == 'run_RBC':
kwargs['prob'] = run_RBC
+ elif sys.argv[i + 1] == 'run_GS':
+ kwargs['prob'] = run_GS
else:
raise NotImplementedError
elif 'num_procs' in sys.argv[i]:
@@ -1667,11 +1672,11 @@ def compare_adaptivity_modes():
def main():
kwargs = {
- 'prob': run_RBC,
+ 'prob': run_vdp,
'num_procs': 1,
'mode': 'default',
- 'runs': 2000,
- 'reload': False,
+ 'runs': 4000,
+ 'reload': True,
**parse_args(),
}
@@ -1690,12 +1695,14 @@ def main():
AdaptivityPolynomialError(**strategy_args),
],
faults=[False, True],
- recovery_thresh=1.15,
+ recovery_thresh=1.1,
recovery_thresh_abs=RECOVERY_THRESH_ABS.get(kwargs.get('prob', None), 0),
stats_path='data/stats-jusuf',
**kwargs,
)
- stats_analyser.run_stats_generation(runs=kwargs['runs'], step=25)
+ stats_analyser.get_recovered()
+
+ stats_analyser.run_stats_generation(runs=kwargs['runs'], step=8)
if MPI.COMM_WORLD.rank > 0: # make sure only one rank accesses the data
return None
diff --git a/pySDC/projects/Resilience/paper_plots.py b/pySDC/projects/Resilience/paper_plots.py
index 67f9aca486846800b8b44460bf08360a0887b6b9..69d176236567ba41acc860f3ea8a7805b4f38964 100644
--- a/pySDC/projects/Resilience/paper_plots.py
+++ b/pySDC/projects/Resilience/paper_plots.py
@@ -10,6 +10,7 @@ from pySDC.projects.Resilience.fault_stats import (
run_quench,
run_AC,
run_RBC,
+ run_GS,
RECOVERY_THRESH_ABS,
)
from pySDC.projects.Resilience.strategies import (
@@ -45,9 +46,9 @@ def get_stats(problem, path='data/stats-jusuf', num_procs=1, strategy_type='SDC'
FaultStats: Object to analyse resilience statistics from
"""
if strategy_type == 'SDC':
- strategies = [BaseStrategy(), AdaptivityStrategy(), IterateStrategy()]
+ strategies = [BaseStrategy(), AdaptivityStrategy(), IterateStrategy(), AdaptivityPolynomialError()]
if JOURNAL not in ['JSC_beamer']:
- strategies += [HotRodStrategy(), AdaptivityPolynomialError()]
+ strategies += [HotRodStrategy()]
elif strategy_type == 'RK':
strategies = [DIRKStrategy()]
if problem.__name__ in ['run_Lorenz', 'run_vdp']:
@@ -142,7 +143,8 @@ def plot_recovery_rate(stats_analyser, **kwargs): # pragma: no cover
None
"""
my_setup_mpl()
- fig, axs = plt.subplots(1, 2, figsize=(TEXTWIDTH, 5 * cm), sharex=True, sharey=True)
+ # fig, axs = plt.subplots(1, 2, figsize=(TEXTWIDTH, 5 * cm), sharex=True, sharey=True)
+ fig, axs = plt.subplots(1, 2, figsize=figsize_by_journal(JOURNAL, 1, 0.4), sharex=True)
stats_analyser.plot_things_per_things(
'recovered',
'bit',
@@ -185,7 +187,7 @@ def plot_recovery_rate_recoverable_only(stats_analyser, fig, ax, **kwargs): # p
ax=ax,
fig=fig,
strategies=[stats_analyser.strategies[i]],
- plotting_args={'markevery': 5},
+ plotting_args={'markevery': 10 if stats_analyser.prob.__name__ in ['run_RBC', 'run_Schroedinger'] else 5},
)
@@ -200,9 +202,9 @@ def compare_recovery_rate_problems(target='resilience', **kwargs): # pragma: no
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']
+ elif target in ['thesis', 'talk']:
+ problems = [run_vdp, run_Lorenz, run_GS, run_RBC] # TODO: swap in Gray-Scott
+ titles = ['Van der Pol', 'Lorenz', 'Gray-Scott', 'Rayleigh-Benard']
else:
raise NotImplementedError()
@@ -219,13 +221,17 @@ def compare_recovery_rate_problems(target='resilience', **kwargs): # pragma: no
ax.get_legend().remove()
if kwargs.get('strategy_type', 'SDC') == 'SDC':
- axs[1, 1].legend(frameon=False, loc="lower right")
+ axs[1, 0].legend(frameon=False, loc="lower right")
else:
axs[0, 1].legend(frameon=False, loc="lower right")
axs[0, 0].set_ylim((-0.05, 1.05))
axs[1, 0].set_ylabel('recovery rate')
axs[0, 0].set_ylabel('recovery rate')
+ if target == 'talk':
+ axs[0, 0].set_xlabel('')
+ axs[0, 1].set_xlabel('')
+
name = ''
for key, val in kwargs.items():
name = f'{name}_{key}-{val}'
@@ -233,6 +239,37 @@ def compare_recovery_rate_problems(target='resilience', **kwargs): # pragma: no
savefig(fig, f'compare_equations{name}.pdf')
+def plot_recovery_rate_detailed_Lorenz(target='resilience'): # pragma: no cover
+ stats = get_stats(run_Lorenz, num_procs=1, strategy_type='SDC')
+ stats.get_recovered()
+ mask = None
+
+ for x in ['node', 'iteration', 'bit']:
+ if target == 'talk':
+ fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.6, 0.6))
+ else:
+ fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.8, 0.5))
+
+ stats.plot_things_per_things(
+ 'recovered',
+ x,
+ False,
+ op=stats.rec_rate,
+ mask=mask,
+ args={'ylabel': 'recovery rate'},
+ ax=ax,
+ plotting_args={'markevery': 5 if x == 'bit' else 1},
+ )
+ ax.set_ylim((-0.05, 1.05))
+
+ if x == 'node':
+ ax.set_xticks([0, 1, 2, 3])
+ elif x == 'iteration':
+ ax.set_xticks([1, 2, 3, 4, 5])
+
+ savefig(fig, f'recovery_rate_Lorenz_{x}')
+
+
def plot_adaptivity_stuff(): # pragma: no cover
"""
Plot the solution for a van der Pol problem as well as the local error and cost associated with the base scheme and
@@ -393,7 +430,7 @@ def plot_fault_vdp(bit=0): # pragma: no cover
savefig(fig, f'fault_bit_{bit}')
-def plot_fault_Lorenz(bit=0): # pragma: no cover
+def plot_fault_Lorenz(bit=0, target='resilience'): # pragma: no cover
"""
Make a plot showing the impact of a fault on the Lorenz attractor without any resilience.
The faults are inserted in the last iteration in the last node in x such that you can best see the impact.
@@ -423,12 +460,15 @@ def plot_fault_Lorenz(bit=0): # pragma: no cover
strategy = BaseStrategy()
my_setup_mpl()
- fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.8, 0.5))
+ if target == 'resilience':
+ fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.4, 0.6))
+ else:
+ fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.8, 0.5))
colors = ['grey', strategy.color, 'magenta']
ls = ['--', '-']
markers = [None, strategy.marker]
do_faults = [False, True]
- superscripts = ['*', '']
+ superscripts = [r'\mathrm{no~faults}', '']
labels = ['x', 'x']
run = 19 + 20 * bit
@@ -447,7 +487,7 @@ def plot_fault_Lorenz(bit=0): # pragma: no cover
[me[1][0] for me in u],
ls=ls[i],
color=colors[i],
- label=rf'${{{labels[i]}}}^{{{superscripts[i]}}}$',
+ label=rf'${{{labels[i]}}}_{{{superscripts[i]}}}$',
marker=markers[i],
markevery=500,
)
@@ -458,8 +498,19 @@ def plot_fault_Lorenz(bit=0): # pragma: no cover
)
ax.set_title(f'Fault in bit {faults[idx][1][4]}')
- ax.legend(frameon=True, loc='lower left')
ax.set_xlabel(r'$t$')
+
+ h, l = ax.get_legend_handles_labels()
+ fig.legend(
+ h,
+ l,
+ loc='outside lower center',
+ ncols=3,
+ frameon=False,
+ fancybox=True,
+ borderaxespad=0.01,
+ )
+
savefig(fig, f'fault_bit_{bit}')
@@ -521,7 +572,7 @@ def plot_quench_solution(): # pragma: no cover
savefig(fig, 'quench_sol')
-def plot_RBC_solution(): # pragma: no cover
+def plot_RBC_solution(setup='resilience'): # pragma: no cover
"""
Plot solution of Rayleigh-Benard convection
"""
@@ -529,13 +580,14 @@ def plot_RBC_solution(): # pragma: no cover
from mpl_toolkits.axes_grid1 import make_axes_locatable
+ nplots = 3 if setup == 'thesis_intro' else 2
+ aspect = 0.8 if nplots == 3 else 0.5
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))
+ fig, axs = plt.subplots(nplots, 1, sharex=True, sharey=True, figsize=figsize_by_journal(JOURNAL, 1.0, aspect))
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)]
+ for ax in axs:
+ divider = make_axes_locatable(ax)
+ caxs += [divider.append_axes('right', size='3%', pad=0.03)]
from pySDC.projects.Resilience.RBC import RayleighBenard, PROBLEM_PARAMS
@@ -547,14 +599,67 @@ def plot_RBC_solution(): # pragma: no cover
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])
+ if setup == 'resilience':
+ _plot(0, axs[0], caxs[0])
+ _plot(21, axs[1], caxs[1])
+ elif setup == 'work-precision':
+ _plot(10, axs[0], caxs[0])
+ _plot(16, axs[1], caxs[1])
+ elif setup == 'resilience_thesis':
+ _plot(20, axs[0], caxs[0])
+ _plot(21, axs[1], caxs[1])
+ elif setup == 'thesis_intro':
+ _plot(0, axs[0], caxs[0])
+ _plot(18, axs[1], caxs[1])
+ _plot(30, axs[2], caxs[2])
- axs[1].set_xlabel('$x$')
- axs[0].set_ylabel('$z$')
- axs[1].set_ylabel('$z$')
+ for ax in axs:
+ ax.set_ylabel('$z$')
+ ax.set_aspect(1)
+ axs[-1].set_xlabel('$x$')
+
+ savefig(fig, f'RBC_sol_{setup}', tight_layout=False)
+
+
+def plot_GS_solution(tend=500): # pragma: no cover
+ my_setup_mpl()
+
+ fig, axs = plt.subplots(1, 2, figsize=figsize_by_journal(JOURNAL, 1.0, 0.45), sharex=True, sharey=True)
+
+ from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+ plt.rcParams['figure.constrained_layout.use'] = True
+ cax = []
+ divider = make_axes_locatable(axs[0])
+ cax += [divider.append_axes('right', size='5%', pad=0.05)]
+ divider2 = make_axes_locatable(axs[1])
+ cax += [divider2.append_axes('right', size='5%', pad=0.05)]
- savefig(fig, 'RBC_sol', tight_layout=False)
+ from pySDC.projects.Resilience.GS import grayscott_imex_diffusion
+
+ problem_params = {
+ 'num_blobs': -48,
+ 'L': 2,
+ 'nvars': (128,) * 2,
+ 'A': 0.062,
+ 'B': 0.1229,
+ 'Du': 2e-5,
+ 'Dv': 1e-5,
+ }
+ P = grayscott_imex_diffusion(**problem_params)
+ Tend = tend
+ im = axs[0].pcolormesh(*P.X, P.u_exact(0)[1], rasterized=True, cmap='binary')
+ im1 = axs[1].pcolormesh(*P.X, P.u_exact(Tend)[1], rasterized=True, cmap='binary')
+
+ fig.colorbar(im, cax=cax[0])
+ fig.colorbar(im1, cax=cax[1])
+ axs[0].set_title(r'$v(t=0)$')
+ axs[1].set_title(rf'$v(t={{{Tend}}})$')
+ for ax in axs:
+ ax.set_aspect(1)
+ ax.set_xlabel('$x$')
+ ax.set_ylabel('$y$')
+ savefig(fig, f'GrayScott_sol{f"_{tend}" if tend != 500 else ""}')
def plot_Schroedinger_solution(): # pragma: no cover
@@ -630,7 +735,7 @@ def plot_AC_solution(): # pragma: no cover
savefig(fig, 'AC_sol')
-def plot_vdp_solution(): # pragma: no cover
+def plot_vdp_solution(setup='adaptivity'): # pragma: no cover
"""
Plot the solution of van der Pol problem over time to illustrate the varying time scales.
@@ -645,24 +750,41 @@ def plot_vdp_solution(): # pragma: no cover
else:
fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 1.0, 0.33))
- custom_description = {
- 'convergence_controllers': {Adaptivity: {'e_tol': 1e-7, 'dt_max': 1e0}},
- 'problem_params': {'mu': 1000, 'crash_at_maxiter': False},
- 'level_params': {'dt': 1e-3},
- }
-
- stats, _, _ = run_vdp(custom_description=custom_description, Tend=2000)
+ if setup == 'adaptivity':
+ custom_description = {
+ 'convergence_controllers': {Adaptivity: {'e_tol': 1e-7, 'dt_max': 1e0}},
+ 'problem_params': {'mu': 1000, 'crash_at_maxiter': False},
+ 'level_params': {'dt': 1e-3},
+ }
+ Tend = 2000
+ elif setup == 'resilience':
+ custom_description = {
+ 'convergence_controllers': {Adaptivity: {'e_tol': 1e-7, 'dt_max': 1e0}},
+ 'problem_params': {'mu': 5, 'crash_at_maxiter': False},
+ 'level_params': {'dt': 1e-3},
+ }
+ Tend = 50
+
+ stats, _, _ = run_vdp(custom_description=custom_description, Tend=Tend)
u = get_sorted(stats, type='u', recomputed=False)
_u = np.array([me[1][0] for me in u])
_x = np.array([me[0] for me in u])
- x1 = _x[abs(_u - 1.1) < 1e-2][0]
- ax.plot(_x, _u, color='black')
- ax.axvspan(x1, x1 + 20, alpha=0.4)
+ name = ''
+ if setup == 'adaptivity':
+ x1 = _x[abs(_u - 1.1) < 1e-2][0]
+ ax.plot(_x, _u, color='black')
+ ax.axvspan(x1, x1 + 20, alpha=0.4)
+ elif setup == 'resilience':
+ x1 = _x[abs(_u - 2.0) < 1e-2][0]
+ ax.plot(_x, _u, color='black')
+ ax.axvspan(x1, x1 + 11.5, alpha=0.4)
+ name = '_resilience'
+
ax.set_ylabel(r'$u$')
ax.set_xlabel(r'$t$')
- savefig(fig, 'vdp_sol')
+ savefig(fig, f'vdp_sol{name}')
def work_precision(): # pragma: no cover
@@ -683,10 +805,19 @@ def work_precision(): # pragma: no cover
all_problems(**{**all_params, 'work_key': 'param'}, mode='compare_strategies')
-def plot_recovery_rate_per_acceptance_threshold(problem): # pragma no cover
+def plot_recovery_rate_per_acceptance_threshold(problem, target='resilience'): # pragma no cover
stats_analyser = get_stats(problem)
- stats_analyser.plot_recovery_thresholds(thresh_range=np.linspace(0.5, 1.5, 1000), recoverable_only=True)
+ if target == 'talk':
+ fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.4, 0.6))
+ else:
+ fig, ax = plt.subplots(figsize=figsize_by_journal(JOURNAL, 0.8, 0.5))
+
+ stats_analyser.plot_recovery_thresholds(thresh_range=np.logspace(-1, 4, 500), recoverable_only=False, ax=ax)
+ ax.set_xscale('log')
+ ax.set_ylim((-0.05, 1.05))
+ ax.set_xlabel('relative threshold')
+ savefig(fig, 'recovery_rate_per_thresh')
def make_plots_for_TIME_X_website(): # pragma: no cover
@@ -739,13 +870,20 @@ def make_plots_for_adaptivity_paper(): # pragma: no cover
def make_plots_for_resilience_paper(): # pragma: no cover
+ global JOURNAL
+ JOURNAL = 'Springer_proceedings'
+
plot_Lorenz_solution()
+ plot_RBC_solution()
+ plot_AC_solution()
+ plot_Schroedinger_solution()
+
plot_fault_Lorenz(0)
plot_fault_Lorenz(20)
- plot_RBC_solution()
compare_recovery_rate_problems(target='resilience', num_procs=1, strategy_type='SDC')
- # plot_recovery_rate(get_stats(run_Lorenz))
- # plot_recovery_rate_per_acceptance_threshold(run_Lorenz)
+ plot_recovery_rate(get_stats(run_Lorenz))
+ plot_recovery_rate_detailed_Lorenz()
+ plot_recovery_rate_per_acceptance_threshold(run_Lorenz)
plt.show()
@@ -764,12 +902,77 @@ def make_plots_for_notes(): # pragma: no cover
def make_plots_for_thesis(): # pragma: no cover
global JOURNAL
JOURNAL = 'TUHH_thesis'
+ for setup in ['thesis_intro', 'resilience_thesis', 'work_precision']:
+ plot_RBC_solution(setup)
- plot_RBC_solution()
- # plot_vdp_solution()
+ from pySDC.projects.Resilience.RBC import plot_factorizations_over_time
+
+ plot_factorizations_over_time(t0=0, Tend=50)
+
+ from pySDC.projects.Resilience.work_precision import all_problems, single_problem
- # plot_adaptivity_stuff()
+ all_params = {
+ 'record': False,
+ 'work_key': 't',
+ 'precision_key': 'e_global_rel',
+ 'plotting': True,
+ 'base_path': 'data/paper',
+ 'target': 'thesis',
+ }
+
+ for mode in ['compare_strategies', 'parallel_efficiency_dt_k', 'parallel_efficiency_dt', 'RK_comp']:
+ all_problems(**all_params, mode=mode)
+ all_problems(**{**all_params, 'work_key': 'param'}, mode='compare_strategies')
+ single_problem(**all_params, mode='RK_comp_high_order_RBC', problem=run_RBC)
+
+ for tend in [500, 2000]:
+ plot_GS_solution(tend=tend)
+ for setup in ['resilience', 'adaptivity']:
+ plot_vdp_solution(setup=setup)
+
+ plot_adaptivity_stuff()
+
+ plot_fault_Lorenz(0)
+ plot_fault_Lorenz(20)
compare_recovery_rate_problems(target='thesis', num_procs=1, strategy_type='SDC')
+ plot_recovery_rate_per_acceptance_threshold(run_Lorenz)
+ plot_recovery_rate(get_stats(run_Lorenz))
+ plot_recovery_rate_detailed_Lorenz()
+
+
+def make_plots_for_TUHH_seminar(): # pragma: no cover
+ global JOURNAL
+ JOURNAL = 'JSC_beamer'
+
+ from pySDC.projects.Resilience.work_precision import (
+ all_problems,
+ )
+
+ all_params = {
+ 'record': False,
+ 'work_key': 't',
+ 'precision_key': 'e_global_rel',
+ 'plotting': True,
+ 'base_path': 'data/paper',
+ 'target': 'talk',
+ }
+
+ for mode in ['compare_strategies', 'parallel_efficiency_dt_k', 'parallel_efficiency_dt', 'RK_comp']:
+ all_problems(**all_params, mode=mode)
+ all_problems(**{**all_params, 'work_key': 'param'}, mode='compare_strategies')
+
+ plot_GS_solution()
+ for setup in ['resilience_thesis', 'work_precision']:
+ plot_RBC_solution(setup)
+ for setup in ['resilience', 'adaptivity']:
+ plot_vdp_solution(setup=setup)
+
+ plot_adaptivity_stuff()
+
+ plot_fault_Lorenz(20, target='talk')
+ compare_recovery_rate_problems(target='talk', num_procs=1, strategy_type='SDC')
+ plot_recovery_rate_per_acceptance_threshold(run_Lorenz, target='talk')
+ plot_recovery_rate_detailed_Lorenz(target='talk')
if __name__ == "__main__":
@@ -777,7 +980,9 @@ if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument(
- '--target', choices=['adaptivity', 'resilience', 'thesis', 'notes', 'SIAM_CSE23', 'TIME_X_website'], type=str
+ '--target',
+ choices=['adaptivity', 'resilience', 'thesis', 'notes', 'SIAM_CSE23', 'TIME_X_website', 'TUHH_seminar'],
+ type=str,
)
args = parser.parse_args()
@@ -793,5 +998,7 @@ if __name__ == "__main__":
make_plots_for_SIAM_CSE23()
elif args.target == 'TIME_X_website':
make_plots_for_TIME_X_website()
+ elif args.target == 'TUHH_seminar':
+ make_plots_for_TUHH_seminar()
else:
raise NotImplementedError(f'Don\'t know how to make plots for target {args.target}')
diff --git a/pySDC/projects/Resilience/reachTendExactly.py b/pySDC/projects/Resilience/reachTendExactly.py
new file mode 100644
index 0000000000000000000000000000000000000000..4f4c394ed03df923dad2b4444cfae1c9be5071d1
--- /dev/null
+++ b/pySDC/projects/Resilience/reachTendExactly.py
@@ -0,0 +1,46 @@
+import numpy as np
+from pySDC.core.convergence_controller import ConvergenceController
+from pySDC.implementations.convergence_controller_classes.check_convergence import CheckConvergence
+
+
+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]
+ time_size = step.status.time_size
+
+ 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) * time_size
+ and not step.status.restart
+ and time_left > 0
+ and step.status.last
+ ):
+ dt_new = (
+ min([(dt + self.params.min_step_size) * time_size, max([time_left, self.params.min_step_size])])
+ + time_size * np.finfo('float').eps * 10
+ ) / time_size
+
+ if dt_new != L.status.dt_new:
+ L.status.dt_new = dt_new
+ self.log(
+ f'Changing step size from {dt:12e} to {L.status.dt_new:.12e} because there is only {time_left:.12e} left.',
+ step,
+ )
diff --git a/pySDC/projects/Resilience/strategies.py b/pySDC/projects/Resilience/strategies.py
index 406e63995d14da2822dd6dc2ad10e65b959fda42..134356331d951ca3ed96a14117c0c27dbf083659 100644
--- a/pySDC/projects/Resilience/strategies.py
+++ b/pySDC/projects/Resilience/strategies.py
@@ -121,7 +121,7 @@ class Strategy:
args['target'] = 0
if problem.__name__ == "run_vdp":
- args['time'] = 5.25
+ args['time'] = 5.5 # 25
elif problem.__name__ == "run_Schroedinger":
args['time'] = 0.3
elif problem.__name__ == "run_quench":
@@ -132,6 +132,8 @@ class Strategy:
args['time'] = 1e-2
elif problem.__name__ == "run_RBC":
args['time'] = 20.19
+ elif problem.__name__ == "run_GS":
+ args['time'] = 100.0
return args
@@ -152,7 +154,7 @@ class Strategy:
rnd_params['iteration'] = base_params['step_params']['maxiter']
rnd_params['rank'] = num_procs
- if problem.__name__ in ['run_Schroedinger', 'run_quench', 'run_AC', 'run_RBC']:
+ if problem.__name__ in ['run_Schroedinger', 'run_quench', 'run_AC', 'run_RBC', 'run_GS']:
rnd_params['min_node'] = 1
if problem.__name__ == "run_quench":
@@ -161,6 +163,8 @@ class Strategy:
rnd_params['iteration'] = 5
elif problem.__name__ == 'run_RBC':
rnd_params['problem_pos'] = [3, 16, 16]
+ elif problem.__name__ == 'run_vdp':
+ rnd_params['iteration'] = 5
return rnd_params
@@ -188,7 +192,7 @@ class Strategy:
return self.name
@classmethod
- def get_Tend(cls, problem, num_procs=1):
+ def get_Tend(cls, problem, num_procs=1, resilience_experiment=False):
'''
Get the final time of runs for fault stats based on the problem
@@ -200,7 +204,10 @@ class Strategy:
float: Tend to put into the run
'''
if problem.__name__ == "run_vdp":
- return 20 # 11.5
+ if resilience_experiment:
+ return 11.5
+ else:
+ return 20
elif problem.__name__ == "run_piline":
return 20.0
elif problem.__name__ == "run_Lorenz":
@@ -213,6 +220,8 @@ class Strategy:
return 0.025
elif problem.__name__ == "run_RBC":
return 21
+ elif problem.__name__ == "run_GS":
+ return 500
else:
raise NotImplementedError('I don\'t have a final time for your problem!')
@@ -277,7 +286,10 @@ class Strategy:
}
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['level_params']['dt'] = 2.5e-2 if num_procs == 4 else 5e-2
+ custom_description['step_params'] = {'maxiter': 5}
+ elif problem.__name__ == 'run_GS':
+ custom_description['level_params']['dt'] = 1.0
custom_description['step_params'] = {'maxiter': 5}
custom_description['convergence_controllers'] = {
@@ -298,6 +310,7 @@ class Strategy:
'run_quench': 150,
'run_AC': 150,
'run_RBC': 1000,
+ 'run_GS': 100,
}
custom_description['convergence_controllers'][StopAtMaxRuntime] = {
@@ -319,14 +332,25 @@ class Strategy:
custom_description = self.get_base_parameters(problem, num_procs)
return merge_descriptions(custom_description, self.custom_description)
- def get_custom_description_for_faults(self, *args, **kwargs):
+ def get_custom_description_for_faults(self, problem, *args, **kwargs):
'''
Get a custom description based on the problem to run the fault stuff
Returns:
dict: Custom description
'''
- return self.get_custom_description(*args, **kwargs)
+ custom_description = self.get_custom_description(problem, *args, **kwargs)
+ if problem.__name__ == "run_vdp":
+ custom_description['step_params'] = {'maxiter': 5}
+ custom_description['problem_params'] = {
+ 'u0': np.array([2.0, 0], dtype=np.float64),
+ 'mu': 5,
+ 'crash_at_maxiter': False,
+ 'newton_tol': 1e-11,
+ 'stop_at_nan': False,
+ }
+ custom_description['level_params'] = {'dt': 1e-2}
+ return custom_description
def get_reference_value(self, problem, key, op, num_procs=1):
"""
@@ -367,7 +391,9 @@ class InexactBaseStrategy(Strategy):
def get_custom_description(self, problem, num_procs=1):
from pySDC.implementations.convergence_controller_classes.inexactness import NewtonInexactness
- preconditioner = 'MIN-SR-NS' if problem.__name__ in ['run_Lorenz'] else 'MIN-SR-S'
+ preconditioner = {
+ 'run_Lorenz': 'MIN-SR-NS',
+ }.get(problem.__name__, 'MIN-SR-S')
desc = {}
desc['sweeper_params'] = {'QI': preconditioner}
@@ -384,7 +410,7 @@ class InexactBaseStrategy(Strategy):
'maxiter': 15,
}
- if self.newton_inexactness and problem.__name__ not in ['run_Schroedinger', 'run_AC', 'run_RBC']:
+ if self.newton_inexactness and problem.__name__ not in ['run_Schroedinger', 'run_AC', 'run_RBC', 'run_GS']:
if problem.__name__ == 'run_quench':
inexactness_params['ratio'] = 1e-1
inexactness_params['min_tol'] = 1e-11
@@ -441,12 +467,24 @@ class BaseStrategy(Strategy):
desc['level_params']['dt'] = 1e-2 * desc['problem_params']['eps'] ** 2
return desc
- def get_custom_description_for_faults(self, problem, *args, **kwargs):
- desc = self.get_custom_description(problem, *args, **kwargs)
+ def get_custom_description_for_faults(self, problem, num_procs, *args, **kwargs):
+ desc = self.get_custom_description(problem, num_procs, *args, **kwargs)
if problem.__name__ == "run_quench":
desc['level_params']['dt'] = 5.0
elif problem.__name__ == "run_AC":
- desc['level_params']['dt'] = 8e-5
+ desc['level_params']['dt'] = 4e-5 if num_procs == 4 else 8e-5
+ elif problem.__name__ == "run_GS":
+ desc['level_params']['dt'] = 4e-1
+ elif problem.__name__ == "run_vdp":
+ desc['step_params'] = {'maxiter': 5}
+ desc['problem_params'] = {
+ 'u0': np.array([2.0, 0], dtype=np.float64),
+ 'mu': 5,
+ 'crash_at_maxiter': False,
+ 'newton_tol': 1e-11,
+ 'stop_at_nan': False,
+ }
+ desc['level_params'] = {'dt': 4.5e-2}
return desc
def get_reference_value(self, problem, key, op, num_procs=1):
@@ -495,24 +533,24 @@ class AdaptivityStrategy(Strategy):
def label(self):
return r'$\Delta t$-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_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):
'''
@@ -535,13 +573,14 @@ class AdaptivityStrategy(Strategy):
dt_max = np.inf
dt_slope_max = np.inf
dt_slope_min = 0
+ beta = 0.9
if problem.__name__ == "run_piline":
e_tol = 1e-7
elif problem.__name__ == "run_vdp":
e_tol = 2e-5
elif problem.__name__ == "run_Lorenz":
- e_tol = 1e-7
+ e_tol = 1e-6 if num_procs == 4 else 1e-7
elif problem.__name__ == "run_Schroedinger":
e_tol = 4e-7
elif problem.__name__ == "run_quench":
@@ -560,8 +599,14 @@ class AdaptivityStrategy(Strategy):
e_tol = 1e-7
# dt_max = 0.1 * base_params['problem_params']['eps'] ** 2
elif problem.__name__ == 'run_RBC':
- e_tol = 1e-4
+ if num_procs == 4:
+ e_tol = 2e-2
+ else:
+ e_tol = 1e-4
dt_slope_min = 1
+ beta = 0.5
+ elif problem.__name__ == 'run_GS':
+ e_tol = 1e-5
else:
raise NotImplementedError(
@@ -573,6 +618,7 @@ class AdaptivityStrategy(Strategy):
'e_tol': e_tol,
'dt_slope_max': dt_slope_max,
'dt_rel_min_slope': dt_slope_min,
+ 'beta': beta,
}
custom_description['convergence_controllers'][StepSizeLimiter] = {
'dt_max': dt_max,
@@ -603,6 +649,8 @@ class AdaptivityStrategy(Strategy):
def get_custom_description_for_faults(self, problem, num_procs, *args, **kwargs):
from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter
from pySDC.implementations.convergence_controller_classes.adaptivity import Adaptivity
+ from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeSlopeLimiter
+ from pySDC.projects.Resilience.RBC import ReachTendExactly
desc = self.get_custom_description(problem, num_procs, *args, **kwargs)
if problem.__name__ == "run_quench":
@@ -610,6 +658,23 @@ class AdaptivityStrategy(Strategy):
desc['convergence_controllers'][Adaptivity]['e_tol'] = 1e-6
elif problem.__name__ == "run_AC":
desc['convergence_controllers'][Adaptivity]['e_tol'] = 1e-5
+ elif problem.__name__ == "run_GS":
+ desc['convergence_controllers'][Adaptivity]['e_tol'] = 2e-6
+ elif problem.__name__ == "run_vdp":
+ desc['step_params'] = {'maxiter': 5}
+ desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'LU'}
+ desc['problem_params'] = {
+ 'u0': np.array([2.0, 0], dtype=np.float64),
+ 'mu': 5,
+ 'crash_at_maxiter': True,
+ 'newton_tol': 1e-8,
+ 'stop_at_nan': True,
+ 'relative_tolerance': False,
+ }
+ desc['level_params'] = {'dt': 8e-3}
+ desc['convergence_controllers'][Adaptivity]['e_tol'] = 2e-7
+ desc['convergence_controllers'][ReachTendExactly] = {'Tend': 11.5}
+ # desc['convergence_controllers'][StepSizeSlopeLimiter] = {'dt_slope_min': 1/4, 'dt_slope_max': 4}
return desc
@@ -778,6 +843,8 @@ class IterateStrategy(Strategy):
restol = 1e-11
elif problem.__name__ == "run_RBC":
restol = 1e-4
+ elif problem.__name__ == "run_GS":
+ restol = 1e-4
else:
raise NotImplementedError(
'I don\'t have a residual tolerance for your problem. Please add one to the \
@@ -851,17 +918,34 @@ class kAdaptivityStrategy(IterateStrategy):
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-6
+ desc['level_params']['e_tol'] = 1e-7
+ elif problem.__name__ == "run_GS":
+ desc['level_params']['dt'] = 1.0
desc['level_params']['restol'] = 1e-9
return desc
- def get_custom_description_for_faults(self, problem, *args, **kwargs):
- desc = self.get_custom_description(problem, *args, **kwargs)
+ def get_custom_description_for_faults(self, problem, num_procs, *args, **kwargs):
+ desc = self.get_custom_description(problem, num_procs, *args, **kwargs)
if problem.__name__ == 'run_quench':
desc['level_params']['dt'] = 5.0
elif problem.__name__ == 'run_AC':
- desc['level_params']['dt'] = 5e-4
+ desc['level_params']['dt'] = 4e-4 if num_procs == 4 else 5e-4
+ desc['level_params']['restol'] = 1e-5 if num_procs == 4 else 1e-11
elif problem.__name__ == 'run_RBC':
- desc['level_params']['restol'] = 1e-6
+ desc['level_params']['restol'] = 1e-3 if num_procs == 4 else 1e-6
+ elif problem.__name__ == 'run_Lorenz':
+ desc['level_params']['dt'] = 8e-3
+ elif problem.__name__ == "run_vdp":
+ desc['sweeper_params'] = {'num_nodes': 3}
+ desc['problem_params'] = {
+ 'u0': np.array([2.0, 0], dtype=np.float64),
+ 'mu': 5,
+ 'crash_at_maxiter': False,
+ 'newton_tol': 1e-11,
+ 'stop_at_nan': False,
+ }
+ desc['level_params'] = {'dt': 4.0e-2, 'restol': 1e-7}
return desc
def get_reference_value(self, problem, key, op, num_procs=1):
@@ -920,14 +1004,15 @@ class HotRodStrategy(Strategy):
base_params = super().get_custom_description(problem, num_procs)
if problem.__name__ == "run_vdp":
- if num_procs == 4:
- HotRod_tol = 1.800804e-04
- elif num_procs == 5:
- HotRod_tol = 9.329361e-05
- else: # 1 process
- HotRod_tol = 1.347949e-06
- HotRod_tol = 7e-6 if num_procs > 1 else 5e-7
- maxiter = 4
+ # if num_procs == 4:
+ # HotRod_tol = 1.800804e-04
+ # elif num_procs == 5:
+ # HotRod_tol = 9.329361e-05
+ # else: # 1 process
+ # HotRod_tol = 1.347949e-06
+ # HotRod_tol = 7e-6 if num_procs > 1 else 5e-7
+ HotRod_tol = 7.2e-05
+ maxiter = 6
elif problem.__name__ == "run_Lorenz":
if num_procs == 5:
HotRod_tol = 9.539348e-06
@@ -956,7 +1041,10 @@ class HotRodStrategy(Strategy):
HotRod_tol = 9.564437e-06
maxiter = 6
elif problem.__name__ == 'run_RBC':
- HotRod_tol = 6.34e-6
+ HotRod_tol = 3e-4 if num_procs == 4 else 6.34e-6
+ maxiter = 6
+ elif problem.__name__ == 'run_GS':
+ HotRod_tol = 3.22e-5
maxiter = 6
else:
raise NotImplementedError(
@@ -986,6 +1074,20 @@ class HotRodStrategy(Strategy):
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'] = 8e-5
+ elif problem.__name__ == "run_GS":
+ desc['level_params']['dt'] = 4e-1
+ elif problem.__name__ == "run_vdp":
+ desc['step_params'] = {'maxiter': 6}
+ desc['problem_params'] = {
+ 'u0': np.array([2.0, 0], dtype=np.float64),
+ 'mu': 5,
+ 'crash_at_maxiter': False,
+ 'newton_tol': 1e-11,
+ 'stop_at_nan': False,
+ }
+ desc['level_params'] = {'dt': 4.5e-2}
return desc
def get_reference_value(self, problem, key, op, num_procs=1):
@@ -1383,6 +1485,63 @@ class ARKStrategy(AdaptivityStrategy):
super().get_reference_value(problem, key, op, num_procs)
+class ARK3_CFL_Strategy(BaseStrategy):
+ """
+ This is special for RBC with CFL number for accuracy
+ """
+
+ def __init__(self, **kwargs):
+ '''
+ Initialization routine
+ '''
+ from pySDC.implementations.problem_classes.RayleighBenard import CFLLimit
+
+ super().__init__(**kwargs)
+ self.color = 'maroon'
+ self.marker = '<'
+ self.name = 'ARK3'
+ self.bar_plot_x_label = 'ARK3'
+ self.precision_parameter = 'cfl'
+ self.precision_parameter_loc = ['convergence_controllers', CFLLimit, 'cfl']
+ self.max_steps = 1e5
+
+ @property
+ def label(self):
+ return 'ARK3'
+
+ def get_custom_description(self, problem, num_procs):
+ '''
+ 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.basic_restarting import BasicRestarting
+ from pySDC.implementations.sweeper_classes.Runge_Kutta import ARK3
+ from pySDC.implementations.problem_classes.RayleighBenard import CFLLimit
+ from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeSlopeLimiter
+
+ desc = super().get_custom_description(problem, num_procs)
+
+ rk_params = {
+ 'step_params': {'maxiter': 1},
+ 'sweeper_class': ARK3,
+ 'convergence_controllers': {
+ CFLLimit: {
+ 'cfl': 0.5,
+ 'dt_max': 1.0,
+ },
+ StepSizeSlopeLimiter: {'dt_rel_min_slope': 0.2},
+ },
+ }
+
+ custom_description = merge_descriptions(desc, rk_params)
+
+ return custom_description
+
+
class ESDIRKStrategy(AdaptivityStrategy):
'''
ESDIRK5(3)
@@ -1515,7 +1674,7 @@ class ERKStrategy(DIRKStrategy):
@property
def label(self):
- return 'CP5(4)'
+ return 'CK5(4)'
def get_random_params(self, problem, num_procs):
'''
@@ -1905,6 +2064,7 @@ class AdaptivityPolynomialError(InexactBaseStrategy):
abort_at_growing_residual = True
level_params = {}
problem_params = {}
+ beta = 0.9
if problem.__name__ == "run_vdp":
e_tol = 6e-4
@@ -1926,16 +2086,23 @@ class AdaptivityPolynomialError(InexactBaseStrategy):
restol_rel = 1e-1
elif problem.__name__ == "run_AC":
e_tol = 1.0e-4
- restol_rel = 1e-3
+ restol_rel = 1e-3 if num_procs == 4 else 1e-3
# dt_max = 0.1 * base_params['problem_params']['eps'] ** 2
elif problem.__name__ == "run_RBC":
- e_tol = 5e-3
+ e_tol = 5e-2 if num_procs == 4 else 5e-3
dt_slope_min = 1.0
abort_at_growing_residual = False
- restol_rel = 1e-3
+ restol_rel = 1e-2 if num_procs == 4 else 1e-4
restol_max = 1e-1
- restol_min = 5e-7
+ restol_min = 5e-8
self.max_slope = 4
+ beta = 0.5
+ level_params['e_tol'] = 1e-5
+ elif problem.__name__ == 'run_GS':
+ e_tol = 1e-4
+ restol_rel = 4e-3
+ restol_max = 1e-4
+ restol_min = 1e-9
else:
raise NotImplementedError(
'I don\'t have a tolerance for adaptivity for your problem. Please add one to the\
@@ -1952,6 +2119,7 @@ class AdaptivityPolynomialError(InexactBaseStrategy):
'factor_if_not_converged': self.max_slope,
'interpolate_between_restarts': self.interpolate_between_restarts,
'abort_at_growing_residual': abort_at_growing_residual,
+ 'beta': beta,
},
StepSizeLimiter: {
'dt_max': dt_max,
@@ -1962,10 +2130,11 @@ class AdaptivityPolynomialError(InexactBaseStrategy):
}
custom_description['level_params'] = level_params
custom_description['problem_params'] = problem_params
+
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)
+ def get_custom_description_for_faults(self, problem, num_procs, *args, **kwargs):
+ desc = self.get_custom_description(problem, num_procs, *args, **kwargs)
if problem.__name__ == "run_quench":
from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError
@@ -1974,7 +2143,31 @@ class AdaptivityPolynomialError(InexactBaseStrategy):
elif problem.__name__ == "run_AC":
from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError
- desc['convergence_controllers'][AdaptivityPolynomialError]['e_tol'] = 1e-3
+ desc['convergence_controllers'][AdaptivityPolynomialError]['e_tol'] = 6e-3 if num_procs == 4 else 1e-3
+ if num_procs == 4:
+ desc['step_params'] = {'maxiter': 50}
+ elif problem.__name__ == "run_Lorenz":
+ from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError
+
+ desc['convergence_controllers'][AdaptivityPolynomialError]['e_tol'] = 2e-4
+ desc['convergence_controllers'][AdaptivityPolynomialError]['restol_min'] = 1e-11
+ desc['convergence_controllers'][AdaptivityPolynomialError]['restol_rel'] = 1e-11
+ elif problem.__name__ == "run_vdp":
+ from pySDC.implementations.convergence_controller_classes.adaptivity import AdaptivityPolynomialError
+ from pySDC.implementations.convergence_controller_classes.inexactness import NewtonInexactness
+
+ desc['step_params'] = {'maxiter': 16}
+ desc['problem_params'] = {
+ 'u0': np.array([2.0, 0], dtype=np.float64),
+ 'mu': 5,
+ 'crash_at_maxiter': False,
+ 'newton_tol': 1e-11,
+ 'stop_at_nan': False,
+ }
+ desc['convergence_controllers'][AdaptivityPolynomialError]['e_tol'] = 5e-4
+ desc['convergence_controllers'][AdaptivityPolynomialError]['restol_rel'] = 8e-5
+ desc['convergence_controllers'].pop(NewtonInexactness)
+ desc['level_params'] = {'dt': 4.5e-2}
return desc
def get_random_params(self, problem, num_procs):
@@ -1993,7 +2186,7 @@ class AdaptivityPolynomialError(InexactBaseStrategy):
if problem.__name__ == "run_quench":
rnd_params['iteration'] = 1
elif problem.__name__ == 'run_Lorenz':
- rnd_params['iteration'] = 4
+ rnd_params['iteration'] = 5
return rnd_params
def get_reference_value(self, problem, key, op, num_procs=1):
diff --git a/pySDC/projects/Resilience/tests/test_ReachTendExactly.py b/pySDC/projects/Resilience/tests/test_ReachTendExactly.py
new file mode 100644
index 0000000000000000000000000000000000000000..1c60d574e12b3a5df9c549f9681091ed70056949
--- /dev/null
+++ b/pySDC/projects/Resilience/tests/test_ReachTendExactly.py
@@ -0,0 +1,73 @@
+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])
+@pytest.mark.parametrize('num_procs', [1, 2])
+def test_ReachTendExactly(Tend, dt, num_procs):
+ if dt * num_procs > Tend:
+ return None
+
+ from pySDC.projects.Resilience.reachTendExactly 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, 1.0 / 3.0, 2)
diff --git a/pySDC/projects/Resilience/work_precision.py b/pySDC/projects/Resilience/work_precision.py
index 86f8600ed849943924d1d9503fcbfb4119004020..2a1c448aa304669cc008ec14951ea7fc6e24cbab 100644
--- a/pySDC/projects/Resilience/work_precision.py
+++ b/pySDC/projects/Resilience/work_precision.py
@@ -13,6 +13,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.GS import run_GS
from pySDC.helpers.stats_helper import get_sorted, filter_stats
from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal
@@ -23,6 +24,11 @@ LOG_TO_FILE = False
logging.getLogger('matplotlib.texmanager').setLevel(90)
+Tends = {'run_RBC': 16.0, 'run_Lorenz': 2.0}
+t0s = {
+ 'run_RBC': 10.0,
+}
+
def std_log(x):
return np.std(np.log(x))
@@ -72,14 +78,6 @@ def get_forbidden_combinations(problem, strategy, **kwargs):
return False
-Tends = {
- 'run_RBC': 16,
-}
-t0s = {
- 'run_RBC': 0.2,
-}
-
-
def single_run(
problem,
strategy,
@@ -161,7 +159,16 @@ def single_run(
t_last = perf_counter()
# record all the metrics
- stats_all = filter_stats(stats, comm=comm_sweep)
+ if comm_sweep.size > 1:
+ try:
+ stats_all = filter_stats(stats, comm=comm_sweep)
+ except MPI.Exception:
+ for key in MAPPINGS.keys():
+ data[key] += [np.nan]
+ return stats
+
+ else:
+ stats_all = stats
comm_sweep.Free()
for key, mapping in MAPPINGS.items():
@@ -289,7 +296,13 @@ def record_work_precision(
else:
exponents = [-3, -2, -1, 0, 0.2, 0.8, 1][::-1]
if problem.__name__ == 'run_RBC':
- exponents = [1, 0, -1, -2, -3, -4]
+ exponents = [1, 0, -0.5, -1, -2]
+ if problem.__name__ == 'run_GS':
+ exponents = [-2, -1, 0, 1, 2, 3][::-1]
+ if problem.__name__ == 'run_Lorenz':
+ exponents = [0, 1, 2, 3][::-1]
+ if type(strategy).__name__ in ["AdaptivityStrategy"]:
+ exponents = [0, 1, 2, 3, 4, 5][::-1]
elif param == 'dt':
power = 2.0
exponents = [-1, 0, 1, 2, 3][::-1]
@@ -298,6 +311,9 @@ def record_work_precision(
exponents = [-2, -1, 0, 1, 2, 3]
if problem.__name__ == 'run_vdp':
exponents = [-4, -3, -2, -1, 0, 1]
+ elif param == 'cfl':
+ power = 2
+ exponents = [-3, -2, -1, 0, 1]
else:
raise NotImplementedError(f"I don't know how to get default value for parameter \"{param}\"")
@@ -312,7 +328,13 @@ def record_work_precision(
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]
+ param_range = [8e-2, 6e-2, 4e-2, 3e-2, 2e-2]
+ if problem.__name__ == 'run_GS':
+ if param == 'dt':
+ param_range = [2, 1, 0.5, 0.1]
+ if problem.__name__ == 'run_Lorenz':
+ if param == 'dt':
+ param_range = [5e-2, 2e-2, 1e-2, 5e-3]
# run multiple times with different parameters
for i in range(len(param_range)):
@@ -495,6 +517,69 @@ def plot_work_precision(
if meta.get('runs', None) == 1:
ax.text(0.1, 0.2, "No sampling!", transform=ax.transAxes)
+ if problem.__name__ == 'run_vdp':
+ if mode == 'parallel_efficiency':
+ # ax.set_xticks([6e-1, 2e0])
+ ax.set_xticks(
+ ticks=[
+ 0.4,
+ 5e-1,
+ 6e-1,
+ 7e-1,
+ 8e-1,
+ 9e-1,
+ 2e0,
+ ],
+ labels=['']
+ + [r'$5\times 10^{-1}$']
+ + [
+ '',
+ ]
+ * 4
+ + [r'$2\times 10^0$'],
+ minor=True,
+ )
+ elif mode == 'RK_comp':
+ ax.set_xticks(
+ ticks=[
+ 5e-1,
+ 6e-1,
+ 7e-1,
+ 8e-1,
+ 9e-1,
+ 2e0,
+ ],
+ labels=[r'$5\times 10^{-1}$']
+ + [
+ '',
+ ]
+ * 4
+ + [r'$2\times 10^0$'],
+ minor=True,
+ )
+ elif problem.__name__ == 'run_quench':
+ if mode == 'RK_comp':
+ ax.set_xticks(
+ ticks=[
+ 0.2,
+ 0.3,
+ 0.4,
+ 5e-1,
+ 6e-1,
+ 7e-1,
+ 8e-1,
+ 9e-1,
+ 2e0,
+ ],
+ labels=['']
+ + [r'$3\times 10^{-1}$']
+ + [
+ '',
+ ]
+ * 7,
+ minor=True,
+ )
+
def plot_parallel_efficiency_diagonalSDC(
ax, work_key, precision_key, num_procs_sweeper, num_procs=1, **kwargs
@@ -589,6 +674,8 @@ def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only
'run_Schroedinger': r'Schr\"odinger',
'run_quench': 'Quench',
'run_AC': 'Allen-Cahn',
+ 'run_RBC': 'Rayleigh-Benard',
+ 'run_GS': 'Gray-Scott',
}
ax.set_title(titles.get(problem.__name__, ''))
@@ -788,8 +875,8 @@ def get_configs(mode, problem):
'sweeper_params': {'num_nodes': 3, 'quad_type': 'RADAU-RIGHT'},
},
'strategies': [
- AdaptivityStrategy(useMPI=True),
BaseStrategy(useMPI=True),
+ AdaptivityStrategy(useMPI=True),
],
}
@@ -803,55 +890,166 @@ def get_configs(mode, problem):
ESDIRKStrategy,
ARKStrategy,
AdaptivityPolynomialError,
+ ARK3_CFL_Strategy,
)
- if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_RBC']:
+ if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_RBC', 'run_GS']:
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', 'run_RBC'] else True
+ newton_inexactness = False if problem.__name__ in ['run_vdp', 'run_RBC', 'run_GS'] else True
desc = {}
desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': "EE"}
desc['step_params'] = {'maxiter': 5}
+ num_procs_dt = {
+ 'run_RBC': 1,
+ }.get(problem.__name__, 4)
desc_poly = {}
desc_poly['sweeper_class'] = parallel_sweeper
+ num_procs_dt_k = 3
ls = {
- 1: '-',
+ 1: '--',
2: '--',
- 3: '-.',
- 4: ':',
- 5: ':',
+ 3: '-',
+ 4: '-',
+ 5: '-',
+ 12: ':',
}
RK_strategies = []
if problem.__name__ in ['run_Lorenz']:
RK_strategies.append(ERKStrategy(useMPI=True))
- if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_RBC']:
+ desc_poly['sweeper_params'] = {'QI': 'MIN-SR-S', 'QE': 'PIC'}
+ desc['sweeper_params']['QI'] = 'MIN-SR-S'
+ desc['sweeper_params']['QE'] = 'PIC'
+ if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS']:
RK_strategies.append(ARKStrategy(useMPI=True))
+ elif problem.__name__ == 'run_RBC':
+ RK_strategies.append(ARK3_CFL_Strategy(useMPI=True))
+ desc['sweeper_params']['num_nodes'] = 2
+ desc['sweeper_params']['QI'] = 'LU'
+ desc['sweeper_params']['QE'] = 'PIC'
+ desc['step_params']['maxiter'] = 3
+
+ desc_poly['sweeper_params'] = {'num_nodes': 2, 'QI': 'MIN-SR-S'}
+ num_procs_dt_k = 2
else:
RK_strategies.append(ESDIRKStrategy(useMPI=True))
- configurations[3] = {
- 'custom_description': desc_poly,
- 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
+ configurations[-1] = {
+ 'strategies': RK_strategies,
'num_procs': 1,
- 'num_procs_sweeper': 3,
- 'plotting_params': {'label': r'$\Delta t$-$k$-adaptivity $N$=1x3'},
}
+ if problem.__name__ == 'run_Lorenz':
+ configurations[3] = {
+ 'custom_description': desc_poly,
+ 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
+ 'num_procs': 4,
+ 'num_procs_sweeper': num_procs_dt_k,
+ 'plotting_params': {
+ 'label': rf'$\Delta t$-$k$-adaptivity $N$=4x{num_procs_dt_k}',
+ 'ls': ls[num_procs_dt_k * 4],
+ },
+ }
+ else:
+ configurations[3] = {
+ 'custom_description': desc_poly,
+ 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
+ 'num_procs': 1,
+ 'num_procs_sweeper': num_procs_dt_k,
+ 'plotting_params': {
+ 'label': rf'$\Delta t$-$k$-adaptivity $N$=1x{num_procs_dt_k}',
+ 'ls': ls[num_procs_dt_k],
+ },
+ }
+ if problem.__name__ in ['run_Lorenz']:
+ configurations[2] = {
+ 'strategies': [AdaptivityStrategy(useMPI=True)],
+ 'custom_description': {**desc, 'sweeper_class': parallel_sweeper},
+ 'num_procs': num_procs_dt,
+ 'num_procs_sweeper': num_procs_dt_k,
+ 'plotting_params': {
+ 'label': rf'$\Delta t$-adaptivity $N$={num_procs_dt}x3',
+ 'ls': ls[num_procs_dt * num_procs_dt_k],
+ },
+ }
+ else:
+ configurations[2] = {
+ 'strategies': [AdaptivityStrategy(useMPI=True)],
+ 'custom_description': desc,
+ 'num_procs': num_procs_dt,
+ 'plotting_params': {'label': rf'$\Delta t$-adaptivity $N$={num_procs_dt}x1', 'ls': ls[num_procs_dt]},
+ }
+
+ elif mode == 'RK_comp_high_order_RBC':
+ """
+ Compare parallel adaptive SDC to Runge-Kutta at order five for RBC problem
+ """
+ from pySDC.projects.Resilience.strategies import (
+ AdaptivityStrategy,
+ ERKStrategy,
+ ESDIRKStrategy,
+ ARKStrategy,
+ AdaptivityPolynomialError,
+ ARK3_CFL_Strategy,
+ )
+
+ assert problem.__name__ == 'run_RBC'
+
+ from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
+
+ newton_inexactness = False
+
+ desc = {}
+ desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': "EE"}
+ desc['step_params'] = {'maxiter': 5}
+ num_procs_dt = 1
+
+ desc_poly = {}
+ desc_poly['sweeper_class'] = parallel_sweeper
+ num_procs_dt_k = 3
+
+ ls = {
+ 1: '--',
+ 2: '--',
+ 3: '-',
+ 4: '-',
+ 5: '-',
+ 12: ':',
+ }
+ RK_strategies = [ARK3_CFL_Strategy(useMPI=True)]
+ desc['sweeper_params']['num_nodes'] = 3
+ desc['sweeper_params']['QI'] = 'LU'
+ desc['sweeper_params']['QE'] = 'PIC'
+ desc['step_params']['maxiter'] = 5
+
+ desc_poly['sweeper_params'] = {'num_nodes': 3, 'QI': 'MIN-SR-S'}
+ num_procs_dt_k = 3
+
configurations[-1] = {
'strategies': RK_strategies,
'num_procs': 1,
}
+ configurations[3] = {
+ 'custom_description': desc_poly,
+ 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
+ 'num_procs': 1,
+ 'num_procs_sweeper': num_procs_dt_k,
+ 'plotting_params': {
+ 'label': rf'$\Delta t$-$k$-adaptivity $N$=1x{num_procs_dt_k}',
+ 'ls': ls[num_procs_dt_k],
+ },
+ }
configurations[2] = {
'strategies': [AdaptivityStrategy(useMPI=True)],
'custom_description': desc,
- 'num_procs': 4,
- 'plotting_params': {'label': r'$\Delta t$-adaptivity $N$=4x1'},
+ 'num_procs': num_procs_dt,
+ 'plotting_params': {'label': rf'$\Delta t$-adaptivity $N$={num_procs_dt}x1', 'ls': ls[num_procs_dt]},
}
elif mode == 'parallel_efficiency':
@@ -860,7 +1058,7 @@ def get_configs(mode, problem):
"""
from pySDC.projects.Resilience.strategies import AdaptivityStrategy, AdaptivityPolynomialError
- if problem.__name__ in ['run_Schroedinger', 'run_AC']:
+ if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS', '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 (
@@ -871,6 +1069,10 @@ def get_configs(mode, problem):
desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': 'EE'}
desc['step_params'] = {'maxiter': 5}
+ if problem.__name__ in ['run_RBC']:
+ desc['sweeper_params']['QE'] = 'PIC'
+ desc['sweeper_params']['QI'] = 'LU'
+
ls = {
1: '-',
2: '--',
@@ -907,13 +1109,137 @@ def get_configs(mode, problem):
},
}
- configurations[num_procs * 200 + 79] = {
+ configurations[200 + 79] = {
'strategies': [
AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True)
],
'num_procs': 1,
}
+ elif mode == 'parallel_efficiency_dt':
+ """
+ Compare parallel runs of the step size adaptive SDC
+ """
+ from pySDC.projects.Resilience.strategies import AdaptivityStrategy
+ if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS', '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,
+ )
+
+ desc = {}
+ desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': 'EE'}
+ desc['step_params'] = {'maxiter': 5}
+
+ if problem.__name__ in ['run_RBC']:
+ desc['sweeper_params']['QE'] = 'PIC'
+ desc['sweeper_params']['QI'] = 'LU'
+
+ desc_serial = {
+ 'step_params': {'maxiter': 5},
+ 'sweeper_params': {'num_nodes': 3, 'quad_type': 'RADAU-RIGHT'},
+ }
+
+ ls = {
+ 1: '-',
+ 2: '--',
+ 3: '-.',
+ 4: '--',
+ 5: ':',
+ 12: ':',
+ }
+
+ newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
+
+ for num_procs in [4, 1]:
+ configurations[num_procs] = {
+ 'strategies': [AdaptivityStrategy(useMPI=True)],
+ 'custom_description': desc.copy() if num_procs > 1 else desc_serial,
+ 'num_procs': num_procs,
+ 'plotting_params': {
+ 'ls': ls.get(num_procs, '-'),
+ 'label': rf'$\Delta t$-adaptivity $N$={num_procs}x1',
+ },
+ }
+ configurations[num_procs * 200 + 79] = {
+ 'custom_description': {
+ 'sweeper_class': parallel_sweeper,
+ 'sweeper_params': {'QI': 'MIN-SR-S', 'QE': 'PIC'},
+ 'step_params': {'maxiter': 5},
+ },
+ 'strategies': [AdaptivityStrategy(useMPI=True)],
+ 'num_procs_sweeper': 3,
+ 'num_procs': num_procs,
+ 'plotting_params': {
+ 'ls': ls.get(num_procs * 3, '-'),
+ 'label': rf'$\Delta t$-adaptivity $N$={num_procs}x3',
+ },
+ }
+ elif mode == 'parallel_efficiency_dt_k':
+ """
+ Compare parallel runs of the step size adaptive SDC
+ """
+ from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
+
+ if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS', '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,
+ )
+
+ ls = {
+ 1: '-',
+ 2: '--',
+ 3: '-.',
+ 4: '--',
+ 5: ':',
+ 12: ':',
+ }
+
+ QI = {
+ (1, 3, 'run_Lorenz'): 'MIN-SR-NS',
+ (1, 1, 'run_Lorenz'): 'MIN-SR-NS',
+ (4, 1, 'run_Lorenz'): 'IE',
+ }
+
+ newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
+
+ for num_procs in [4, 1]:
+ configurations[num_procs * 100 + 79] = {
+ 'custom_description': {
+ 'sweeper_class': parallel_sweeper,
+ 'sweeper_params': {'QI': QI.get((num_procs, 3, problem.__name__), 'MIN-SR-S'), 'QE': 'PIC'},
+ },
+ 'strategies': [
+ AdaptivityPolynomialError(
+ useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True
+ )
+ ],
+ 'num_procs_sweeper': 3,
+ 'num_procs': num_procs,
+ 'plotting_params': {
+ 'ls': ls.get(num_procs * 3, '-'),
+ 'label': rf'$\Delta t$-$k$-adaptivity $N$={num_procs}x3',
+ },
+ }
+ configurations[num_procs * 200 + 79] = {
+ 'custom_description': {
+ 'sweeper_params': {'QI': QI.get((num_procs, 1, problem.__name__), 'MIN-SR-S'), 'QE': 'PIC'},
+ },
+ 'strategies': [
+ AdaptivityPolynomialError(
+ useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True
+ )
+ ],
+ 'num_procs_sweeper': 1,
+ 'num_procs': num_procs,
+ 'plotting_params': {
+ 'ls': ls.get(num_procs, '-'),
+ 'label': rf'$\Delta t$-$k$-adaptivity $N$={num_procs}x1',
+ },
+ }
elif mode == 'interpolate_between_restarts':
"""
Compare adaptivity with interpolation between restarts and without
@@ -1298,21 +1624,31 @@ def get_configs(mode, problem):
return configurations
-def get_fig(x=1, y=1, **kwargs): # pragma: no cover
+def get_fig(x=1, y=1, target='adaptivity', **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
+ target (str): Where the plot is supposed to end up
Returns:
matplotlib.pyplot.Figure
"""
width = 1.0
ratio = 1.0 if y == 2 else 0.5
+ if target == 'adaptivity':
+ journal = 'Springer_Numerical_Algorithms'
+ elif target == 'thesis':
+ journal = 'TUHH_thesis'
+ elif target == 'talk':
+ journal = 'JSC_beamer'
+ else:
+ raise NotImplementedError
+
keyword_arguments = {
- 'figsize': figsize_by_journal('Springer_Numerical_Algorithms', width, ratio),
+ 'figsize': figsize_by_journal(journal, width, ratio),
'layout': 'constrained',
**kwargs,
}
@@ -1363,7 +1699,9 @@ def save_fig(
print(f'Stored figure \"{path}\"')
-def all_problems(mode='compare_strategies', plotting=True, base_path='data', **kwargs): # pragma: no cover
+def all_problems(
+ mode='compare_strategies', plotting=True, base_path='data', target='adaptivity', **kwargs
+): # pragma: no cover
"""
Make a plot comparing various strategies for all problems.
@@ -1375,7 +1713,10 @@ def all_problems(mode='compare_strategies', plotting=True, base_path='data', **k
None
"""
- fig, axs = get_fig(2, 2)
+ if target == 'talk':
+ fig, axs = get_fig(4, 1, target=target)
+ else:
+ fig, axs = get_fig(2, 2, target=target)
shared_params = {
'work_key': 'k_SDC',
@@ -1388,7 +1729,12 @@ def all_problems(mode='compare_strategies', plotting=True, base_path='data', **k
**kwargs,
}
- problems = [run_vdp, run_quench, run_Schroedinger, run_AC]
+ if target == 'adaptivity':
+ problems = [run_vdp, run_quench, run_Schroedinger, run_AC]
+ elif target in ['thesis', 'talk']:
+ problems = [run_vdp, run_Lorenz, run_GS, run_RBC]
+ else:
+ raise NotImplementedError
logger.log(26, f"Doing for all problems {mode}")
for i in range(len(problems)):
@@ -1404,17 +1750,18 @@ def all_problems(mode='compare_strategies', plotting=True, base_path='data', **k
if plotting and shared_params['comm_world'].rank == 0:
ncols = {
'parallel_efficiency': 2,
+ 'parallel_efficiency_dt': 2,
+ 'parallel_efficiency_dt_k': 2,
'RK_comp': 2,
}
- y_right_dt_fixed = [1e18, 4e1, 5, 1e8]
- y_right_dt = [1e-1, 1e4, 1, 2e-2]
- y_right_dtk = [1e-4, 1e4, 1e-2, 1e-3]
+ if target == 'talk':
+ _ncols = 4
+ else:
+ _ncols = ncols.get(mode, None)
if shared_params['work_key'] == 'param':
- for ax, yRfixed, yRdt, yRdtk in zip(fig.get_axes(), y_right_dt_fixed, y_right_dt, y_right_dtk):
- add_order_line(ax, 1, '--', yRdt, marker=None)
- add_order_line(ax, 5 / 4, ':', yRdtk, marker=None)
- add_order_line(ax, 5, '-.', yRfixed, marker=None)
+ for ax, prob in zip(fig.get_axes(), problems):
+ add_param_order_lines(ax, prob)
save_fig(
fig=fig,
name=mode,
@@ -1422,10 +1769,46 @@ def all_problems(mode='compare_strategies', plotting=True, base_path='data', **k
precision_key=shared_params['precision_key'],
legend=True,
base_path=base_path,
- ncols=ncols.get(mode, None),
+ ncols=_ncols,
)
+def add_param_order_lines(ax, problem):
+ if problem.__name__ == 'run_vdp':
+ yRfixed = 1e18
+ yRdt = 1e-1
+ yRdtk = 1e-4
+ elif problem.__name__ == 'run_quench':
+ yRfixed = 4e1
+ yRdt = 1e4
+ yRdtk = 1e4
+ elif problem.__name__ == 'run_Schroedinger':
+ yRfixed = 5
+ yRdt = 1
+ yRdtk = 1e-2
+ elif problem.__name__ == 'run_AC':
+ yRfixed = 1e8
+ yRdt = 2e-2
+ yRdtk = 1e-3
+ elif problem.__name__ == 'run_Lorenz':
+ yRfixed = 1e1
+ yRdt = 2e-2
+ yRdtk = 7e-4
+ elif problem.__name__ == 'run_RBC':
+ yRfixed = 1e-6
+ yRdt = 4e-5
+ yRdtk = 8e-6
+ elif problem.__name__ == 'run_GS':
+ yRfixed = 4e-3
+ yRdt = 5e0
+ yRdtk = 8e-1
+ else:
+ return None
+ add_order_line(ax, 1, '--', yRdt, marker=None)
+ add_order_line(ax, 5 / 4, ':', yRdtk, marker=None)
+ add_order_line(ax, 5, '-.', yRfixed, marker=None)
+
+
def ODEs(mode='compare_strategies', plotting=True, base_path='data', **kwargs): # pragma: no cover
"""
Make a plot comparing various strategies for the two ODEs.
@@ -1473,7 +1856,7 @@ def ODEs(mode='compare_strategies', plotting=True, base_path='data', **kwargs):
)
-def single_problem(mode, problem, plotting=True, base_path='data', **kwargs): # pragma: no cover
+def single_problem(mode, problem, plotting=True, base_path='data', target='thesis', **kwargs): # pragma: no cover
"""
Make a plot for a single problem
@@ -1481,7 +1864,10 @@ def single_problem(mode, problem, plotting=True, base_path='data', **kwargs): #
mode (str): What you want to look at
problem (function): A problem to run
"""
- fig, ax = get_fig(1, 1, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1, 0.8))
+ if target == 'thesis':
+ fig, ax = get_fig(1, 1, figsize=figsize_by_journal('TUHH_thesis', 0.7, 0.6))
+ else:
+ fig, ax = get_fig(1, 1, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1, 0.8))
params = {
'work_key': 'k_SDC',
@@ -1507,6 +1893,7 @@ def single_problem(mode, problem, plotting=True, base_path='data', **kwargs): #
precision_key=params['precision_key'],
legend=False,
base_path=base_path,
+ squares=target != 'thesis',
)
@@ -1603,38 +1990,42 @@ def aggregate_parallel_efficiency_plot(): # pragma: no cover
if __name__ == "__main__":
comm_world = MPI.COMM_WORLD
- record = comm_world.size > 1
- for mode in [
- # 'compare_strategies',
- # 'RK_comp',
- # 'parallel_efficiency',
- ]:
- params = {
- 'mode': mode,
- 'runs': 1,
- 'plotting': comm_world.rank == 0,
- }
- params_single = {
- **params,
- 'problem': run_RBC,
- }
- single_problem(**params_single, work_key='t', precision_key='e_global_rel', record=record)
-
- all_params = {
- 'record': True,
- 'runs': 5,
- 'work_key': 't',
- 'precision_key': 'e_global_rel',
- 'plotting': comm_world.rank == 0,
+ import argparse
+
+ parser = argparse.ArgumentParser()
+ parser.add_argument('--mode', type=str, default='compare_strategies')
+ parser.add_argument('--record', type=str, choices=['True', 'False'], default='True')
+ parser.add_argument('--plotting', type=str, choices=['True', 'False'], default='True')
+ parser.add_argument('--runs', type=int, default=5)
+ parser.add_argument(
+ '--problem', type=str, choices=['vdp', 'RBC', 'AC', 'quench', 'Lorenz', 'Schroedinger', 'GS'], default='vdp'
+ )
+ parser.add_argument('--work_key', type=str, default='t')
+ parser.add_argument('--precision_key', type=str, default='e_global_rel')
+ parser.add_argument('--logger_level', type=int, default='25')
+
+ args = parser.parse_args()
+
+ problems = {
+ 'Lorenz': run_Lorenz,
+ 'vdp': run_vdp,
+ 'Schroedinger': run_Schroedinger,
+ 'quench': run_quench,
+ 'AC': run_AC,
+ 'RBC': run_RBC,
+ 'GS': run_GS,
+ }
+
+ params = {
+ **vars(args),
+ 'record': args.record == 'True',
+ 'plotting': args.plotting == 'True' and comm_world.rank == 0,
+ 'problem': problems[args.problem],
}
- for mode in [
- 'RK_comp',
- 'parallel_efficiency',
- 'compare_strategies',
- ]:
- all_problems(**all_params, mode=mode)
- comm_world.Barrier()
+ LOGGER_LEVEL = params.pop('logger_level')
+
+ single_problem(**params)
if comm_world.rank == 0:
plt.show()