Coverage for pySDC/projects/Resilience/work_precision.py: 0%
310 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1from mpi4py import MPI
2import numpy as np
3import matplotlib.pyplot as plt
4import pickle
5import logging
6from time import perf_counter
7import copy
9from pySDC.projects.Resilience.strategies import merge_descriptions
10from pySDC.projects.Resilience.Lorenz import run_Lorenz
11from pySDC.projects.Resilience.vdp import run_vdp
12from pySDC.projects.Resilience.Schroedinger import run_Schroedinger
13from pySDC.projects.Resilience.quench import run_quench
14from pySDC.projects.Resilience.AC import run_AC
16from pySDC.helpers.stats_helper import get_sorted, filter_stats
17from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal
19setup_mpl(reset=True)
20LOGGER_LEVEL = 25
21LOG_TO_FILE = False
23logging.getLogger('matplotlib.texmanager').setLevel(90)
25MAPPINGS = {
26 'e_global': ('e_global_post_run', max, False),
27 'e_global_rel': ('e_global_rel_post_run', max, False),
28 't': ('timing_run', max, False),
29 # 'e_local_max': ('e_local_post_step', max, False),
30 'k_SDC': ('k', sum, None),
31 'k_SDC_no_restart': ('k', sum, False),
32 'k_Newton': ('work_newton', sum, None),
33 'k_linear': ('work_linear', sum, None),
34 'k_Newton_no_restart': ('work_newton', sum, False),
35 'k_rhs': ('work_rhs', sum, None),
36 'num_steps': ('dt', len, None),
37 'restart': ('restart', sum, None),
38 'dt_mean': ('dt', np.mean, False),
39 'dt_max': ('dt', max, False),
40 'dt_min': ('dt', min, False),
41 'dt_sigma': ('dt', np.std, False),
42 'e_embedded_max': ('error_embedded_estimate', max, False),
43 'u0_increment_max': ('u0_increment', max, None),
44 'u0_increment_mean': ('u0_increment', np.mean, None),
45 'u0_increment_max_no_restart': ('u0_increment', max, False),
46 'u0_increment_mean_no_restart': ('u0_increment', np.mean, False),
47}
49logger = logging.getLogger('WorkPrecision')
50logger.setLevel(LOGGER_LEVEL)
53def get_forbidden_combinations(problem, strategy, **kwargs):
54 """
55 Check if the combination of strategy and problem is forbidden
57 Args:
58 problem (function): A problem to run
59 strategy (Strategy): SDC strategy
60 """
61 if strategy.name == 'ERK':
62 if problem.__name__ in ['run_quench', 'run_Schroedinger', 'run_AC']:
63 return True
65 return False
68def single_run(
69 problem,
70 strategy,
71 data,
72 custom_description,
73 num_procs=1,
74 comm_world=None,
75 problem_args=None,
76 hooks=None,
77 Tend=None,
78 num_procs_sweeper=1,
79):
80 """
81 Make a single run of a particular problem with a certain strategy.
83 Args:
84 problem (function): A problem to run
85 strategy (Strategy): SDC strategy
86 data (dict): Put the results in here
87 custom_description (dict): Overwrite presets
88 num_procs (int): Number of processes for the time communicator
89 comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script
90 hooks (list): List of additional hooks
91 num_procs_sweeper (int): Number of processes for the sweeper
93 Returns:
94 dict: Stats generated by the run
95 """
96 from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun
97 from pySDC.implementations.hooks.log_work import LogWork
98 from pySDC.projects.Resilience.hook import LogData
100 hooks = hooks if hooks else []
102 t_last = perf_counter()
104 num_procs_tot = num_procs * num_procs_sweeper
105 comm = comm_world.Split(comm_world.rank < num_procs_tot)
106 if comm_world.rank >= num_procs_tot:
107 comm.Free()
108 return None
110 # make communicators for time and sweepers
111 comm_time = comm.Split(comm.rank // num_procs)
112 comm_sweep = comm.Split(comm_time.rank)
114 strategy_description = strategy.get_custom_description(problem, num_procs)
115 description = merge_descriptions(strategy_description, custom_description)
116 if comm_sweep.size > 1:
117 description['sweeper_params']['comm'] = comm_sweep
119 controller_params = {
120 'logger_level': LOGGER_LEVEL,
121 'log_to_file': LOG_TO_FILE,
122 'fname': 'out.txt',
123 **strategy.get_controller_params(),
124 }
125 problem_args = {} if problem_args is None else problem_args
127 stats, controller, crash = problem(
128 custom_description=description,
129 Tend=strategy.get_Tend(problem, num_procs) if Tend is None else Tend,
130 hook_class=[LogData, LogWork, LogGlobalErrorPostRun] + hooks,
131 custom_controller_params=controller_params,
132 use_MPI=True,
133 comm=comm_time,
134 **problem_args,
135 )
137 t_now = perf_counter()
138 logger.debug(f'Finished run in {t_now - t_last:.2e} s')
139 t_last = perf_counter()
141 # record all the metrics
142 stats_all = filter_stats(stats, comm=comm_sweep)
143 comm_sweep.Free()
145 for key, mapping in MAPPINGS.items():
146 if crash:
147 data[key] += [np.nan]
148 continue
149 me = get_sorted(stats_all, comm=comm_time, type=mapping[0], recomputed=mapping[2])
150 if len(me) == 0:
151 data[key] += [np.nan]
152 else:
153 data[key] += [mapping[1]([you[1] for you in me])]
155 t_now = perf_counter()
156 logger.debug(f'Recorded all data after {t_now - t_last:.2e} s')
157 t_last = perf_counter()
159 comm_time.Free()
160 comm.Free()
161 return stats
164def get_parameter(dictionary, where):
165 """
166 Get a parameter at a certain position in a dictionary of dictionaries.
168 Args:
169 dictionary (dict): The dictionary
170 where (list): The list of keys leading to the value you want
172 Returns:
173 The value of the dictionary
174 """
175 if len(where) == 1:
176 return dictionary[where[0]]
177 else:
178 return get_parameter(dictionary[where[0]], where[1:])
181def set_parameter(dictionary, where, parameter):
182 """
183 Set a parameter at a certain position in a dictionary of dictionaries
185 Args:
186 dictionary (dict): The dictionary
187 where (list): The list of keys leading to the value you want to set
188 parameter: Whatever you want to set the parameter to
190 Returns:
191 None
192 """
193 if len(where) == 1:
194 dictionary[where[0]] = parameter
195 else:
196 set_parameter(dictionary[where[0]], where[1:], parameter)
199def get_path(problem, strategy, num_procs, handle='', base_path='data/work_precision', num_procs_sweeper=1, mode=''):
200 """
201 Get the path to a certain data.
203 Args:
204 problem (function): A problem to run
205 strategy (Strategy): SDC strategy
206 num_procs (int): Number of processes for the time communicator
207 handle (str): The name of the configuration
208 base_path (str): Some path where all the files are stored
209 num_procs_sweeper (int): Number of processes for the sweeper
210 mode (str): The mode this was generated for
212 Returns:
213 str: The path to the data you are looking for
214 """
215 return f'{base_path}/{problem.__name__}-{strategy.__class__.__name__}-{handle}{"-wp" if handle else "wp"}-{num_procs}-{num_procs_sweeper}procs-{mode}.pickle'
218def record_work_precision(
219 problem,
220 strategy,
221 num_procs=1,
222 custom_description=None,
223 handle='',
224 runs=1,
225 comm_world=None,
226 problem_args=None,
227 param_range=None,
228 Tend=None,
229 hooks=None,
230 num_procs_sweeper=1,
231 mode='',
232):
233 """
234 Run problem with strategy and record the cost parameters.
236 Args:
237 problem (function): A problem to run
238 strategy (Strategy): SDC strategy
239 num_procs (int): Number of processes for the time communicator
240 custom_description (dict): Overwrite presets
241 handle (str): The name of the configuration
242 runs (int): Number of runs you want to do
243 comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script
244 num_procs_sweeper (int): Number of processes for the sweeper
246 Returns:
247 None
248 """
249 if get_forbidden_combinations(problem, strategy):
250 return None
252 data = {}
254 # prepare precision parameters
255 param = strategy.precision_parameter
256 description = merge_descriptions(
257 strategy.get_custom_description(problem, num_procs),
258 {} if custom_description is None else custom_description,
259 )
260 if param == 'e_tol':
261 power = 10.0
262 set_parameter(description, strategy.precision_parameter_loc[:-1] + ['dt_min'], 0)
263 exponents = [-3, -2, -1, 0, 1, 2, 3][::-1]
264 if problem.__name__ == 'run_vdp':
265 exponents = [-4, -3, -2, -1, 0, 1, 2]
266 elif param == 'dt':
267 power = 2.0
268 exponents = [-1, 0, 1, 2, 3][::-1]
269 elif param == 'restol':
270 power = 10.0
271 exponents = [-2, -1, 0, 1, 2, 3]
272 if problem.__name__ == 'run_vdp':
273 exponents = [-4, -3, -2, -1, 0, 1]
274 else:
275 raise NotImplementedError(f"I don't know how to get default value for parameter \"{param}\"")
277 where = strategy.precision_parameter_loc
278 default = get_parameter(description, where)
279 param_range = [default * power**i for i in exponents] if param_range is None else param_range
281 if problem.__name__ == 'run_quench':
282 if param == 'restol':
283 param_range = [1e-5, 1e-6, 1e-7, 1e-8, 1e-9]
284 elif param == 'dt':
285 param_range = [1.25, 2.5, 5.0, 10.0, 20.0][::-1]
287 # run multiple times with different parameters
288 for i in range(len(param_range)):
289 set_parameter(description, where, param_range[i])
291 data[param_range[i]] = {key: [] for key in MAPPINGS.keys()}
292 data[param_range[i]]['param'] = [param_range[i]]
293 data[param_range[i]][param] = [param_range[i]]
295 description = merge_descriptions(
296 descA=description, descB=strategy.get_description_for_tolerance(problem=problem, param=param_range[i])
297 )
298 for _j in range(runs):
299 if comm_world.rank == 0:
300 logger.log(
301 24,
302 f'Starting: {problem.__name__}: {strategy} {handle} {num_procs}-{num_procs_sweeper} procs, {param}={param_range[i]:.2e}',
303 )
304 single_run(
305 problem,
306 strategy,
307 data[param_range[i]],
308 custom_description=description,
309 comm_world=comm_world,
310 problem_args=problem_args,
311 num_procs=num_procs,
312 hooks=hooks,
313 Tend=Tend,
314 num_procs_sweeper=num_procs_sweeper,
315 )
317 comm_world.Barrier()
319 if comm_world.rank == 0:
320 if np.isfinite(data[param_range[i]]["k_linear"][-1]):
321 k_type = "k_linear"
322 elif np.isfinite(data[param_range[i]]["k_Newton"][-1]):
323 k_type = 'k_Newton'
324 else:
325 k_type = "k_SDC"
326 logger.log(
327 25,
328 f'{problem.__name__}: {strategy} {handle} {num_procs}-{num_procs_sweeper} procs, {param}={param_range[i]:.2e}: e={data[param_range[i]]["e_global"][-1]}, t={data[param_range[i]]["t"][-1]}, {k_type}={data[param_range[i]][k_type][-1]}',
329 )
331 if comm_world.rank == 0:
332 import socket
333 import time
335 data['meta'] = {
336 'hostname': socket.gethostname(),
337 'time': time.time,
338 'runs': runs,
339 }
340 path = get_path(problem, strategy, num_procs, handle, num_procs_sweeper=num_procs_sweeper, mode=mode)
341 with open(path, 'wb') as f:
342 logger.debug(f'Dumping file \"{path}\"')
343 pickle.dump(data, f)
344 return data
347def load(**kwargs):
348 """
349 Load stored data. Arguments are passed on to `get_path`
351 Returns:
352 dict: The data
353 """
354 path = get_path(**kwargs)
355 with open(path, 'rb') as f:
356 logger.debug(f'Loading file \"{path}\"')
357 data = pickle.load(f)
358 return data
361def extract_data(data, work_key, precision_key):
362 """
363 Get the work and precision from a data object.
365 Args:
366 data (dict): Data from work-precision measurements
367 work_key (str): Name of variable on x-axis
368 precision_key (str): Name of variable on y-axis
370 Returns:
371 list: Work
372 list: Precision
373 """
374 keys = [key for key in data.keys() if key not in ['meta']]
375 work = [np.nanmean(data[key][work_key]) for key in keys]
376 precision = [np.nanmean(data[key][precision_key]) for key in keys]
377 return work, precision
380def get_order(work_key='e_global', precision_key='param', strategy=None, handle=None, **kwargs):
381 data = load(**kwargs, strategy=strategy, handle=handle)
382 work, precision = extract_data(data, work_key, precision_key)
384 order = [np.log(precision[i + 1] / precision[i]) / np.log(work[i + 1] / work[i]) for i in range(len(work) - 1)]
386 print(f'Order for {strategy} {handle}: {np.mean(order):.2f}')
389def plot_work_precision(
390 problem,
391 strategy,
392 num_procs,
393 ax,
394 work_key='k_SDC',
395 precision_key='e_global',
396 handle='',
397 plotting_params=None,
398 comm_world=None,
399 num_procs_sweeper=1,
400 mode='',
401): # pragma: no cover
402 """
403 Plot data from running a problem with a strategy.
405 Args:
406 problem (function): A problem to run
407 strategy (Strategy): SDC strategy
408 num_procs (int): Number of processes for the time communicator
409 ax (matplotlib.pyplot.axes): Somewhere to plot
410 work_key (str): The key in the recorded data you want on the x-axis
411 precision_key (str): The key in the recorded data you want on the y-axis
412 handle (str): The name of the configuration
413 plotting_params (dict): Will be passed when plotting
414 comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script
415 num_procs_sweeper (int): Number of processes for the sweeper
416 mode (str): The of the configurations you want to retrieve
418 Returns:
419 None
420 """
421 if comm_world.rank > 0 or get_forbidden_combinations(problem, strategy):
422 return None
424 data = load(
425 problem=problem,
426 strategy=strategy,
427 num_procs=num_procs,
428 handle=handle,
429 num_procs_sweeper=num_procs_sweeper,
430 mode=mode,
431 )
433 work, precision = extract_data(data, work_key, precision_key)
434 keys = [key for key in data.keys() if key not in ['meta']]
436 for key in [work_key, precision_key]:
437 rel_variance = [np.std(data[me][key]) / max([np.nanmean(data[me][key]), 1.0]) for me in keys]
438 if not all(me < 1e-1 or not np.isfinite(me) for me in rel_variance):
439 logger.warning(
440 f"Variance in \"{key}\" for {get_path(problem, strategy, num_procs, handle, num_procs_sweeper=num_procs_sweeper, mode=mode)} too large! Got {rel_variance}"
441 )
443 style = merge_descriptions(
444 {**strategy.style, 'label': f'{strategy.style["label"]}{f" {handle}" if handle else ""}'},
445 plotting_params if plotting_params else {},
446 )
448 ax.loglog(work, precision, **style)
450 # get_order(
451 # problem=problem,
452 # strategy=strategy,
453 # num_procs=num_procs,
454 # handle=handle,
455 # work_key=work_key,
456 # precision_key=precision_key,
457 # )
459 if 't' in [work_key, precision_key]:
460 meta = data.get('meta', {})
462 if meta.get('hostname', None) in ['thomas-work']:
463 ax.text(0.1, 0.1, "Laptop timings!", transform=ax.transAxes)
464 if meta.get('runs', None) == 1:
465 ax.text(0.1, 0.2, "No sampling!", transform=ax.transAxes)
468def plot_parallel_efficiency_diagonalSDC(
469 ax, work_key, precision_key, num_procs_sweeper, num_procs=1, **kwargs
470): # pragma: no cover
471 serial_data = load(
472 num_procs=num_procs,
473 num_procs_sweeper=1,
474 **kwargs,
475 )
476 parallel_data = load(
477 num_procs=num_procs,
478 num_procs_sweeper=num_procs_sweeper,
479 **kwargs,
480 )
481 serial_work, serial_precision = extract_data(serial_data, work_key, precision_key)
482 parallel_work, parallel_precision = extract_data(parallel_data, work_key, precision_key)
483 # assert np.allclose(serial_precision, parallel_precision)
485 serial_work = np.asarray(serial_work)
486 parallel_work = np.asarray(parallel_work)
488 # ax.loglog(serial_work, serial_precision)
489 # ax.loglog(parallel_work, parallel_precision)
491 speedup = serial_work / parallel_work
492 parallel_efficiency = np.median(speedup) / num_procs_sweeper
493 ax.plot(serial_precision, speedup)
494 ax.set_xscale('log')
495 ax.set_ylabel('speedup')
497 if 't' in [work_key, precision_key]:
498 meta = parallel_data.get('meta', {})
500 if meta.get('hostname', None) in ['thomas-work']:
501 ax.text(0.1, 0.1, "Laptop timings!", transform=ax.transAxes)
502 if meta.get('runs', None) == 1:
503 ax.text(0.1, 0.2, "No sampling!", transform=ax.transAxes)
505 return np.median(speedup), parallel_efficiency
508def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only=False): # pragma: no cover
509 """
510 Decorate a plot
512 Args:
513 ax (matplotlib.pyplot.axes): Somewhere to plot
514 problem (function): A problem to run
515 work_key (str): The key in the recorded data you want on the x-axis
516 precision_key (str): The key in the recorded data you want on the y-axis
517 num_procs (int): Number of processes for the time communicator
518 title_only (bool): Put only the title on top, or do the whole shebang
520 Returns:
521 None
522 """
523 labels = {
524 'k_SDC': 'SDC iterations',
525 'k_SDC_no_restart': 'SDC iterations (restarts excluded)',
526 'k_Newton': 'Newton iterations',
527 'k_Newton_no_restart': 'Newton iterations (restarts excluded)',
528 'k_rhs': 'right hand side evaluations',
529 't': 'wall clock time / s',
530 'e_global': 'global error',
531 'e_global_rel': 'relative global error',
532 'e_local_max': 'max. local error',
533 'restart': 'restarts',
534 'dt_max': r'$\Delta t_\mathrm{max}$',
535 'dt_min': r'$\Delta t_\mathrm{min}$',
536 'dt_sigma': r'$\sigma(\Delta t)$',
537 'dt_mean': r'$\bar{\Delta t}$',
538 'param': 'parameter',
539 'u0_increment_max': r'$\| \Delta u_0 \|_{\infty} $',
540 'u0_increment_mean': r'$\bar{\Delta u_0}$',
541 'u0_increment_max_no_restart': r'$\| \Delta u_0 \|_{\infty} $ (restarts excluded)',
542 'u0_increment_mean_no_restart': r'$\bar{\Delta u_0}$ (restarts excluded)',
543 'k_linear': 'Linear solver iterations',
544 'speedup': 'Speedup',
545 'nprocs': r'$N_\mathrm{procs}$',
546 '': '',
547 }
549 if not title_only:
550 ax.set_xlabel(labels.get(work_key, 'work'))
551 ax.set_ylabel(labels.get(precision_key, 'precision'))
552 # ax.legend(frameon=False)
554 titles = {
555 'run_vdp': 'Van der Pol',
556 'run_Lorenz': 'Lorenz attractor',
557 'run_Schroedinger': r'Schr\"odinger',
558 'run_quench': 'Quench',
559 'run_AC': 'Allen-Cahn',
560 }
561 ax.set_title(titles.get(problem.__name__, ''))
564def execute_configurations(
565 problem,
566 configurations,
567 work_key,
568 precision_key,
569 num_procs,
570 ax,
571 decorate,
572 record,
573 runs,
574 comm_world,
575 plotting,
576 Tend=None,
577 num_procs_sweeper=1,
578 mode='',
579):
580 """
581 Run for multiple configurations.
583 Args:
584 problem (function): A problem to run
585 configurations (dict): The configurations you want to run with
586 work_key (str): The key in the recorded data you want on the x-axis
587 precision_key (str): The key in the recorded data you want on the y-axis
588 num_procs (int): Number of processes for the time communicator
589 ax (matplotlib.pyplot.axes): Somewhere to plot
590 decorate (bool): Whether to decorate fully or only put the title
591 record (bool): Whether to only plot or also record the data first
592 runs (int): Number of runs you want to do
593 comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script
594 plotting (bool): Whether to plot something
595 num_procs_sweeper (int): Number of processes for the sweeper
596 mode (str): What you want to look at
598 Returns:
599 None
600 """
601 for _, config in configurations.items():
602 for strategy in config['strategies']:
603 shared_args = {
604 'problem': problem,
605 'strategy': strategy,
606 'handle': config.get('handle', ''),
607 'num_procs': config.get('num_procs', num_procs),
608 'num_procs_sweeper': config.get('num_procs_sweeper', num_procs_sweeper),
609 }
610 if record:
611 logger.debug('Recording work precision')
612 record_work_precision(
613 **shared_args,
614 custom_description=config.get('custom_description', {}),
615 runs=runs,
616 comm_world=comm_world,
617 problem_args=config.get('problem_args', {}),
618 param_range=config.get('param_range', None),
619 hooks=config.get('hooks', None),
620 Tend=Tend,
621 mode=mode,
622 )
623 if plotting and comm_world.rank == 0:
624 logger.debug('Plotting')
625 plot_work_precision(
626 **shared_args,
627 work_key=work_key,
628 precision_key=precision_key,
629 ax=ax,
630 plotting_params=config.get('plotting_params', {}),
631 comm_world=comm_world,
632 mode=mode,
633 )
635 if comm_world.rank == 0:
636 decorate_panel(
637 ax=ax,
638 problem=problem,
639 work_key=work_key,
640 precision_key=precision_key,
641 num_procs=num_procs,
642 title_only=not decorate,
643 )
646def get_configs(mode, problem):
647 """
648 Get configurations for work-precision plots. These are dictionaries containing strategies and handles and so on.
650 Args:
651 mode (str): The of the configurations you want to retrieve
652 problem (function): A problem to run
654 Returns:
655 dict: Configurations
656 """
657 configurations = {}
658 if mode == 'regular':
659 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, BaseStrategy, IterateStrategy
661 handle = 'regular'
662 configurations[0] = {
663 'handle': handle,
664 'strategies': [AdaptivityStrategy(useMPI=True), BaseStrategy(useMPI=True), IterateStrategy(useMPI=True)],
665 }
666 elif mode == 'step_size_limiting':
667 from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter
668 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, ESDIRKStrategy
670 limits = [
671 25.0,
672 50.0,
673 ]
674 colors = ['teal', 'magenta']
675 markers = ['v', 'x']
676 markersize = 9
677 for i in range(len(limits)):
678 configurations[i] = {
679 'custom_description': {'convergence_controllers': {StepSizeLimiter: {'dt_max': limits[i]}}},
680 'handle': f'steplimiter{limits[i]:.0f}',
681 'strategies': [AdaptivityStrategy(useMPI=True)],
682 'plotting_params': {
683 'color': colors[i],
684 'marker': markers[i],
685 'label': rf'$\Delta t \leq { {limits[i]:.0f}} $',
686 # 'ls': '',
687 'markersize': markersize,
688 },
689 'num_procs': 1,
690 }
691 configurations[99] = {
692 'custom_description': {},
693 'handle': 'no limits',
694 'plotting_params': {
695 'label': 'no limiter',
696 # 'ls': '',
697 'markersize': markersize,
698 },
699 'strategies': [AdaptivityStrategy(useMPI=True)],
700 'num_procs': 1,
701 }
702 elif mode == 'dynamic_restarts':
703 """
704 Compare Block Gauss-Seidel SDC with restarting the first step in the block or the first step that exceeded the error threshold.
705 """
706 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, AdaptivityRestartFirstStep
708 desc = {}
709 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'}
710 desc['step_params'] = {'maxiter': 5}
712 ls = {
713 1: '-',
714 2: '--',
715 3: '-.',
716 4: ':',
717 5: ':',
718 }
720 configurations[-1] = {
721 'strategies': [AdaptivityStrategy(useMPI=True)],
722 'num_procs': 1,
723 }
725 for num_procs in [4, 2]:
726 plotting_params = {'ls': ls[num_procs], 'label': f'adaptivity {num_procs} procs'}
727 configurations[num_procs] = {
728 'strategies': [AdaptivityStrategy(useMPI=True), AdaptivityRestartFirstStep(useMPI=True)],
729 'custom_description': desc,
730 'num_procs': num_procs,
731 'plotting_params': plotting_params,
732 }
734 elif mode == 'compare_strategies':
735 """
736 Compare the different SDC strategies.
737 """
738 from pySDC.projects.Resilience.strategies import (
739 AdaptivityStrategy,
740 kAdaptivityStrategy,
741 AdaptivityPolynomialError,
742 BaseStrategy,
743 )
745 configurations[2] = {
746 'strategies': [kAdaptivityStrategy(useMPI=True)],
747 }
748 configurations[1] = {
749 'strategies': [AdaptivityPolynomialError(useMPI=True)],
750 }
751 configurations[0] = {
752 'custom_description': {
753 'step_params': {'maxiter': 5},
754 'sweeper_params': {'num_nodes': 3, 'quad_type': 'RADAU-RIGHT'},
755 },
756 'strategies': [
757 AdaptivityStrategy(useMPI=True),
758 BaseStrategy(useMPI=True),
759 ],
760 }
762 elif mode == 'interpolate_between_restarts':
763 """
764 Compare adaptivity with interpolation between restarts and without
765 """
766 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
768 i = 0
769 for interpolate_between_restarts, handle, ls in zip(
770 [True, False], ['Interpolation between restarts', 'regular'], ['--', '-']
771 ):
772 configurations[i] = {
773 'strategies': [
774 AdaptivityPolynomialError(interpolate_between_restarts=interpolate_between_restarts, useMPI=True)
775 ],
776 'plotting_params': {'ls': ls},
777 'handle': handle,
778 }
779 i += 1
780 elif mode == 'diagonal_SDC':
781 """
782 Run diagonal SDC with different number of nodes and ranks. You can use this to compute a speedup, but it's not strong scaling!
783 """
784 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
786 if problem.__name__ in ['run_Schroedinger']:
787 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
788 else:
789 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
790 generic_implicit_MPI as parallel_sweeper,
791 )
793 for parallel in [False, True]:
794 desc = {'sweeper_class': parallel_sweeper} if parallel else {}
795 for num_nodes, ls in zip([3, 4, 2], ['-', '--', ':', '-.']):
796 configurations[num_nodes + (99 if parallel else 0)] = {
797 'custom_description': {**desc, 'sweeper_params': {'num_nodes': num_nodes}},
798 'strategies': [
799 AdaptivityPolynomialError(useMPI=True, newton_inexactness=True, linear_inexactness=True)
800 ],
801 'num_procs_sweeper': num_nodes if parallel else 1,
802 'num_procs': 1,
803 'handle': f'{num_nodes} nodes',
804 'plotting_params': {
805 'ls': ls,
806 'label': f'{num_nodes} procs',
807 # **{'color': 'grey' if parallel else None},
808 },
809 }
811 elif mode == 'parallel_efficiency':
812 """
813 Compare parallel runs of the step size adaptive SDC
814 """
815 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, AdaptivityPolynomialError
817 if problem.__name__ in ['run_Schroedinger', 'run_AC']:
818 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
819 else:
820 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
821 generic_implicit_MPI as parallel_sweeper,
822 )
824 desc = {}
825 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': 'EE'}
826 desc['step_params'] = {'maxiter': 5}
828 ls = {
829 1: '-',
830 2: '--',
831 3: '-.',
832 4: '--',
833 5: ':',
834 12: ':',
835 }
837 for num_procs in [4, 1]:
838 plotting_params = (
839 {'ls': ls[num_procs], 'label': fr'$\Delta t$-adaptivity $N$={num_procs}x1'} if num_procs > 1 else {}
840 )
841 configurations[num_procs] = {
842 'strategies': [AdaptivityStrategy(useMPI=True)],
843 'custom_description': desc.copy(),
844 'num_procs': num_procs,
845 'plotting_params': plotting_params.copy(),
846 }
847 configurations[num_procs * 100 + 79] = {
848 'custom_description': {'sweeper_class': parallel_sweeper},
849 'strategies': [
850 AdaptivityPolynomialError(useMPI=True, newton_inexactness=True, linear_inexactness=True)
851 ],
852 'num_procs_sweeper': 3,
853 'num_procs': num_procs,
854 'plotting_params': {
855 'ls': ls.get(num_procs * 3, '-'),
856 'label': rf'$\Delta t$-$k$-adaptivity $N$={num_procs}x3',
857 },
858 }
860 configurations[num_procs * 200 + 79] = {
861 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=True, linear_inexactness=True)],
862 'num_procs': 1,
863 }
865 elif mode[:13] == 'vdp_stiffness':
866 """
867 Run van der Pol with different parameter for the nonlinear term, which controls the stiffness.
868 """
869 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, ERKStrategy, ESDIRKStrategy
871 mu = float(mode[14:])
873 problem_desc = {'problem_params': {'mu': mu}}
875 desc = {}
876 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'}
877 desc['step_params'] = {'maxiter': 5}
878 desc['problem_params'] = problem_desc['problem_params']
880 ls = {
881 1: '-',
882 2: '--',
883 3: '-.',
884 4: ':',
885 5: ':',
886 }
888 for num_procs in [5]:
889 plotting_params = {'ls': ls[num_procs], 'label': f'GSSDC {num_procs} procs'}
890 configurations[num_procs] = {
891 'strategies': [AdaptivityStrategy(True)],
892 'custom_description': desc,
893 'num_procs': num_procs,
894 'plotting_params': plotting_params,
895 'handle': mode,
896 }
898 configurations[1] = {
899 'strategies': [AdaptivityStrategy(True)],
900 'custom_description': desc,
901 'num_procs': 1,
902 'plotting_params': {'ls': ls[1], 'label': 'SDC'},
903 'handle': mode,
904 }
906 configurations[2] = {
907 'strategies': [ERKStrategy(useMPI=True)],
908 'num_procs': 1,
909 'handle': mode,
910 'plotting_params': {'label': 'CP5(4)'},
911 'custom_description': problem_desc,
912 }
913 configurations[4] = {
914 'strategies': [ESDIRKStrategy(useMPI=True)],
915 'num_procs': 1,
916 'handle': mode,
917 'plotting_params': {'label': 'ESDIRK5(3)'},
918 'custom_description': problem_desc,
919 }
920 elif mode == 'inexactness':
921 """
922 Compare inexact SDC to exact SDC
923 """
924 from pySDC.projects.Resilience.strategies import (
925 AdaptivityPolynomialError,
926 )
928 if problem.__name__ in ['run_Schroedinger']:
929 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
930 else:
931 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
932 generic_implicit_MPI as parallel_sweeper,
933 )
935 strategies = [
936 AdaptivityPolynomialError,
937 ]
939 inexactness = {
940 'newton_inexactness': True,
941 'linear_inexactness': True,
942 }
943 no_inexactness = {
944 'newton_inexactness': False,
945 'linear_inexactness': False,
946 'SDC_maxiter': 99,
947 'use_restol_rel': False,
948 }
950 configurations[1] = {
951 'custom_description': {'sweeper_class': parallel_sweeper},
952 'strategies': [me(useMPI=True, **no_inexactness) for me in strategies],
953 'num_procs_sweeper': 3,
954 'handle': 'exact',
955 'plotting_params': {'ls': '--'},
956 }
957 configurations[0] = {
958 'custom_description': {'sweeper_class': parallel_sweeper},
959 'strategies': [me(useMPI=True, **inexactness) for me in strategies],
960 'handle': 'inexact',
961 'num_procs_sweeper': 3,
962 }
963 elif mode == 'compare_adaptivity':
964 """
965 Compare various modes of adaptivity
966 """
967 # TODO: configurations not final!
968 from pySDC.projects.Resilience.strategies import (
969 # AdaptivityCollocationTypeStrategy,
970 # AdaptivityCollocationRefinementStrategy,
971 AdaptivityStrategy,
972 # AdaptivityExtrapolationWithinQStrategy,
973 ESDIRKStrategy,
974 ARKStrategy,
975 AdaptivityPolynomialError,
976 )
978 if problem.__name__ in ['run_Schroedinger']:
979 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
980 else:
981 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
982 generic_implicit_MPI as parallel_sweeper,
983 )
985 inexactness_params = {
986 # 'double_adaptivity': True,
987 'newton_inexactness': True,
988 'linear_inexactness': True,
989 }
991 strategies = [
992 AdaptivityPolynomialError,
993 # AdaptivityCollocationTypeStrategy,
994 # AdaptivityExtrapolationWithinQStrategy,
995 ]
997 # restol = None
998 # for strategy in strategies:
999 # strategy.restol = restol
1001 configurations[1] = {
1002 'custom_description': {'sweeper_class': parallel_sweeper},
1003 'strategies': [me(useMPI=True, **inexactness_params) for me in strategies],
1004 'handle': 'parallel',
1005 'num_procs_sweeper': 3,
1006 'plotting_params': {'ls': '-', 'label': '3 procs'},
1007 }
1008 configurations[2] = {
1009 'strategies': [me(useMPI=True, **inexactness_params) for me in strategies],
1010 'plotting_params': {'ls': '--'},
1011 }
1012 configurations[4] = {
1013 'custom_description': {'step_params': {'maxiter': 5}},
1014 'strategies': [AdaptivityStrategy(useMPI=True)],
1015 }
1017 desc_RK = {}
1018 configurations[-1] = {
1019 'strategies': [
1020 ARKStrategy(useMPI=True) if problem.__name__ == 'run_Schroedinger' else ESDIRKStrategy(useMPI=True),
1021 ],
1022 'num_procs': 1,
1023 'custom_description': desc_RK,
1024 }
1026 elif mode == 'preconditioners':
1027 """
1028 Compare different preconditioners
1029 """
1030 from pySDC.projects.Resilience.strategies import (
1031 AdaptivityStrategy,
1032 IterateStrategy,
1033 BaseStrategy,
1034 ESDIRKStrategy,
1035 ERKStrategy,
1036 AdaptivityPolynomialError,
1037 )
1039 inexacness = True
1040 strategies = [
1041 AdaptivityPolynomialError(
1042 useMPI=True, SDC_maxiter=29, newton_inexactness=inexacness, linear_inexactness=inexacness
1043 ),
1044 BaseStrategy(useMPI=True),
1045 ]
1047 desc = {}
1048 desc['sweeper_params'] = {
1049 'num_nodes': 3,
1050 }
1051 # desc['step_params'] = {'maxiter': 5}
1053 precons = ['IE', 'LU']
1054 ls = ['-.', '--', '-', ':']
1055 for i in range(len(precons) + 1):
1056 if i < len(precons):
1057 desc['sweeper_params']['QI'] = precons[i]
1058 handle = precons[i]
1059 else:
1060 handle = None
1061 configurations[i] = {
1062 'strategies': strategies,
1063 'custom_description': copy.deepcopy(desc),
1064 'handle': handle,
1065 'plotting_params': {'ls': ls[i]},
1066 }
1068 elif mode == 'RK_comp':
1069 """
1070 Compare parallel adaptive SDC to Runge-Kutta
1071 """
1072 from pySDC.projects.Resilience.strategies import (
1073 AdaptivityStrategy,
1074 ERKStrategy,
1075 ESDIRKStrategy,
1076 ARKStrategy,
1077 AdaptivityPolynomialError,
1078 )
1080 if problem.__name__ in ['run_Schroedinger', 'run_AC']:
1081 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1082 else:
1083 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1084 generic_implicit_MPI as parallel_sweeper,
1085 )
1087 desc = {}
1088 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': "EE"}
1089 desc['step_params'] = {'maxiter': 5}
1091 desc_poly = {}
1092 desc_poly['sweeper_class'] = parallel_sweeper
1094 ls = {
1095 1: '-',
1096 2: '--',
1097 3: '-.',
1098 4: ':',
1099 5: ':',
1100 }
1102 configurations[3] = {
1103 'custom_description': desc_poly,
1104 'strategies': [AdaptivityPolynomialError(useMPI=True)],
1105 'num_procs': 1,
1106 'num_procs_sweeper': 3,
1107 'plotting_params': {'label': r'$\Delta t$-$k$-adaptivity $N$=1x3'},
1108 }
1109 configurations[-1] = {
1110 'strategies': [
1111 ERKStrategy(useMPI=True),
1112 (
1113 ARKStrategy(useMPI=True)
1114 if problem.__name__ in ['run_Schroedinger', 'run_AC']
1115 else ESDIRKStrategy(useMPI=True)
1116 ),
1117 ],
1118 'num_procs': 1,
1119 }
1120 configurations[2] = {
1121 'strategies': [AdaptivityStrategy(useMPI=True)],
1122 'custom_description': desc,
1123 'num_procs': 4,
1124 'plotting_params': {'label': r'$\Delta t$-adaptivity $N$=4x1'},
1125 }
1127 elif mode == 'RK_comp_high_order':
1128 """
1129 Compare higher order SDC than we can get with RKM to RKM
1130 """
1131 from pySDC.projects.Resilience.strategies import (
1132 AdaptivityStrategy,
1133 ERKStrategy,
1134 ESDIRKStrategy,
1135 ARKStrategy,
1136 AdaptivityPolynomialError,
1137 )
1139 if problem.__name__ in ['run_Schroedinger']:
1140 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1141 else:
1142 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1143 generic_implicit_MPI as parallel_sweeper,
1144 )
1146 desc = {}
1147 desc['sweeper_params'] = {'num_nodes': 4, 'QI': 'IE', 'QE': "EE"}
1148 desc['step_params'] = {'maxiter': 7}
1150 desc_poly = {}
1151 desc_poly['sweeper_class'] = parallel_sweeper
1152 desc_poly['sweeper_params'] = {'num_nodes': 4}
1154 ls = {
1155 1: '-',
1156 2: '--',
1157 3: '-.',
1158 4: ':',
1159 5: ':',
1160 }
1162 desc_RK = {}
1163 if problem.__name__ in ['run_Schroedinger']:
1164 desc_RK['problem_params'] = {'imex': True}
1166 configurations[3] = {
1167 'custom_description': desc_poly,
1168 'strategies': [AdaptivityPolynomialError(useMPI=True)],
1169 'num_procs': 1,
1170 'num_procs_sweeper': 4,
1171 }
1172 configurations[-1] = {
1173 'strategies': [
1174 ERKStrategy(useMPI=True),
1175 ARKStrategy(useMPI=True) if problem.__name__ in ['run_Schroedinger'] else ESDIRKStrategy(useMPI=True),
1176 ],
1177 'num_procs': 1,
1178 'custom_description': desc_RK,
1179 }
1181 configurations[2] = {
1182 'strategies': [AdaptivityStrategy(useMPI=True)],
1183 'custom_description': desc,
1184 'num_procs': 4,
1185 }
1186 elif mode == 'avoid_restarts':
1187 """
1188 Test how well avoiding restarts works.
1189 """
1190 from pySDC.projects.Resilience.strategies import (
1191 AdaptivityStrategy,
1192 AdaptivityAvoidRestartsStrategy,
1193 AdaptivityPolynomialStrategy,
1194 )
1196 desc = {'sweeper_params': {'QI': 'IE'}, 'step_params': {'maxiter': 3}}
1197 param_range = [1e-3, 1e-5]
1198 configurations[0] = {
1199 'strategies': [AdaptivityPolynomialStrategy(useMPI=True)],
1200 'plotting_params': {'ls': '--'},
1201 'custom_description': desc,
1202 'param_range': param_range,
1203 }
1204 configurations[1] = {
1205 'strategies': [AdaptivityAvoidRestartsStrategy(useMPI=True)],
1206 'plotting_params': {'ls': '-.'},
1207 'custom_description': desc,
1208 'param_range': param_range,
1209 }
1210 configurations[2] = {
1211 'strategies': [AdaptivityStrategy(useMPI=True)],
1212 'custom_description': desc,
1213 'param_range': param_range,
1214 }
1215 else:
1216 raise NotImplementedError(f'Don\'t know the mode "{mode}"!')
1218 return configurations
1221def get_fig(x=1, y=1, **kwargs): # pragma: no cover
1222 """
1223 Get a figure to plot in.
1225 Args:
1226 x (int): How many panels in horizontal direction you want
1227 y (int): How many panels in vertical direction you want
1229 Returns:
1230 matplotlib.pyplot.Figure
1231 """
1232 width = 1.0
1233 ratio = 1.0 if y == 2 else 0.5
1234 keyword_arguments = {
1235 'figsize': figsize_by_journal('Springer_Numerical_Algorithms', width, ratio),
1236 'layout': 'constrained',
1237 **kwargs,
1238 }
1239 return plt.subplots(y, x, **keyword_arguments)
1242def save_fig(
1243 fig, name, work_key, precision_key, legend=True, format='pdf', base_path='data', squares=True, ncols=None, **kwargs
1244): # pragma: no cover
1245 """
1246 Save a figure with a legend on the bottom.
1248 Args:
1249 fig (matplotlib.pyplot.Figure): Figure you want to save
1250 name (str): Name of the plot to put in the path
1251 work_key (str): The key in the recorded data you want on the x-axis
1252 precision_key (str): The key in the recorded data you want on the y-axis
1253 legend (bool): Put a legend or not
1254 format (str): Format to store the figure with
1255 base_path (str): Some path where all the files are stored
1256 squares (bool): Adjust aspect ratio to squares if true
1258 Returns:
1259 None
1260 """
1261 handles = []
1262 labels = []
1263 for ax in fig.get_axes():
1264 h, l = ax.get_legend_handles_labels()
1265 handles += [h[i] for i in range(len(h)) if l[i] not in labels]
1266 labels += [me for me in l if me not in labels]
1267 if squares:
1268 ax.set_box_aspect(1)
1269 order = np.argsort([me[0] for me in labels])
1270 fig.legend(
1271 [handles[i] for i in order],
1272 [labels[i] for i in order],
1273 loc='outside lower center',
1274 ncols=ncols if ncols else 3 if len(handles) % 3 == 0 else 4,
1275 frameon=False,
1276 fancybox=True,
1277 handlelength=2.2,
1278 )
1280 path = f'{base_path}/wp-{name}-{work_key}-{precision_key}.{format}'
1281 fig.savefig(path, bbox_inches='tight', **kwargs)
1282 print(f'Stored figure \"{path}\"')
1285def all_problems(mode='compare_strategies', plotting=True, base_path='data', **kwargs): # pragma: no cover
1286 """
1287 Make a plot comparing various strategies for all problems.
1289 Args:
1290 work_key (str): The key in the recorded data you want on the x-axis
1291 precision_key (str): The key in the recorded data you want on the y-axis
1293 Returns:
1294 None
1295 """
1297 fig, axs = get_fig(2, 2)
1299 shared_params = {
1300 'work_key': 'k_SDC',
1301 'precision_key': 'e_global',
1302 'num_procs': 1,
1303 'runs': 1,
1304 'comm_world': MPI.COMM_WORLD,
1305 'record': False,
1306 'plotting': plotting,
1307 **kwargs,
1308 }
1310 problems = [run_vdp, run_quench, run_Schroedinger, run_AC]
1312 logger.log(26, f"Doing for all problems {mode}")
1313 for i in range(len(problems)):
1314 execute_configurations(
1315 **shared_params,
1316 problem=problems[i],
1317 ax=axs.flatten()[i],
1318 decorate=True,
1319 configurations=get_configs(mode, problems[i]),
1320 mode=mode,
1321 )
1323 if plotting and shared_params['comm_world'].rank == 0:
1324 save_fig(
1325 fig=fig,
1326 name=mode,
1327 work_key=shared_params['work_key'],
1328 precision_key=shared_params['precision_key'],
1329 legend=True,
1330 base_path=base_path,
1331 ncols=3 if mode in ['parallel_efficiency'] else None,
1332 )
1335def ODEs(mode='compare_strategies', plotting=True, base_path='data', **kwargs): # pragma: no cover
1336 """
1337 Make a plot comparing various strategies for the two ODEs.
1339 Args:
1340 work_key (str): The key in the recorded data you want on the x-axis
1341 precision_key (str): The key in the recorded data you want on the y-axis
1343 Returns:
1344 None
1345 """
1347 fig, axs = get_fig(x=2, y=1)
1349 shared_params = {
1350 'work_key': 'k_SDC',
1351 'precision_key': 'e_global',
1352 'num_procs': 1,
1353 'runs': 1,
1354 'comm_world': MPI.COMM_WORLD,
1355 'record': False,
1356 'plotting': plotting,
1357 **kwargs,
1358 }
1360 problems = [run_vdp, run_Lorenz]
1362 for i in range(len(problems)):
1363 execute_configurations(
1364 **shared_params,
1365 problem=problems[i],
1366 ax=axs.flatten()[i],
1367 decorate=i == 0,
1368 configurations=get_configs(mode, problems[i]),
1369 )
1371 if plotting and shared_params['comm_world'].rank == 0:
1372 save_fig(
1373 fig=fig,
1374 name=f'ODEs-{mode}',
1375 work_key=shared_params['work_key'],
1376 precision_key=shared_params['precision_key'],
1377 legend=True,
1378 base_path=base_path,
1379 )
1382def single_problem(mode, problem, plotting=True, base_path='data', **kwargs): # pragma: no cover
1383 """
1384 Make a plot for a single problem
1386 Args:
1387 mode (str): What you want to look at
1388 problem (function): A problem to run
1389 """
1390 fig, ax = get_fig(1, 1, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1, 0.8))
1392 params = {
1393 'work_key': 'k_SDC',
1394 'precision_key': 'e_global',
1395 'num_procs': 1,
1396 'runs': 1,
1397 'comm_world': MPI.COMM_WORLD,
1398 'record': False,
1399 'plotting': plotting,
1400 **kwargs,
1401 }
1403 logger.log(26, f"Doing single problem {mode}")
1404 execute_configurations(
1405 **params, problem=problem, ax=ax, decorate=True, configurations=get_configs(mode, problem), mode=mode
1406 )
1408 if plotting:
1409 save_fig(
1410 fig=fig,
1411 name=f'{problem.__name__}-{mode}',
1412 work_key=params['work_key'],
1413 precision_key=params['precision_key'],
1414 legend=False,
1415 base_path=base_path,
1416 )
1419def vdp_stiffness_plot(base_path='data', format='pdf', **kwargs): # pragma: no cover
1420 fig, axs = get_fig(2, 2, sharex=True, sharey=True)
1422 # mus = [0, 5, 10, 15]
1423 mus = [0, 10, 20, 40]
1425 for i in range(len(mus)):
1426 params = {
1427 'runs': 1,
1428 'problem': run_vdp,
1429 'record': False,
1430 'work_key': 't',
1431 'precision_key': 'e_global',
1432 'comm_world': MPI.COMM_WORLD,
1433 **kwargs,
1434 }
1435 params['num_procs'] = min(params['comm_world'].size, 5)
1436 params['plotting'] = params['comm_world'].rank == 0
1438 configurations = get_configs(mode=f'vdp_stiffness-{mus[i]}', problem=run_vdp)
1439 execute_configurations(**params, ax=axs.flatten()[i], decorate=True, configurations=configurations, Tend=100)
1440 axs.flatten()[i].set_title(rf'$\mu={ {mus[i]}} $')
1442 fig.suptitle('Van der Pol')
1443 if params['comm_world'].rank == 0:
1444 save_fig(
1445 fig=fig,
1446 name='vdp-stiffness',
1447 work_key=params['work_key'],
1448 precision_key=params['precision_key'],
1449 legend=False,
1450 base_path=base_path,
1451 format=format,
1452 )
1455def aggregate_parallel_efficiency_plot(): # pragma: no cover
1456 """
1457 Make a "weak scaling" plot for diagonal SDC
1458 """
1459 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
1461 fig, axs = plt.subplots(2, 2)
1463 _fig, _ax = plt.subplots(1, 1)
1464 num_procs = 1
1465 num_procs_sweeper = 2
1466 problem = run_quench
1468 num_procs_sweeper_list = [2, 3, 4]
1470 for problem, ax in zip([run_vdp, run_Lorenz, run_quench], axs.flatten()):
1471 speedup = []
1472 for num_procs_sweeper in num_procs_sweeper_list:
1473 s, e = plot_parallel_efficiency_diagonalSDC(
1474 ax=_ax,
1475 work_key='t',
1476 precision_key='e_global_rel',
1477 num_procs=num_procs,
1478 num_procs_sweeper=num_procs_sweeper,
1479 problem=problem,
1480 strategy=AdaptivityPolynomialError(),
1481 mode='diagonal_SDC',
1482 handle=f'{num_procs_sweeper} nodes',
1483 )
1484 speedup += [s]
1485 decorate_panel(ax, problem, work_key='nprocs', precision_key='')
1487 ax.plot(num_procs_sweeper_list, speedup, label='speedup')
1488 ax.plot(
1489 num_procs_sweeper_list,
1490 [speedup[i] / num_procs_sweeper_list[i] for i in range(len(speedup))],
1491 label='parallel efficiency',
1492 )
1494 fig.tight_layout()
1495 save_fig(fig, 'parallel_efficiency', 'nprocs', 'speedup')
1498if __name__ == "__main__":
1499 comm_world = MPI.COMM_WORLD
1501 # record = False
1502 # for mode in [
1503 # 'compare_strategies',
1504 # 'RK_comp',
1505 # 'parallel_efficiency',
1506 # ]:
1507 # params = {
1508 # 'mode': mode,
1509 # 'runs': 5,
1510 # 'plotting': comm_world.rank == 0,
1511 # }
1512 # params_single = {
1513 # **params,
1514 # 'problem': run_AC,
1515 # }
1516 # single_problem(**params_single, work_key='t', precision_key='e_global_rel', record=record)
1518 all_params = {
1519 'record': True,
1520 'runs': 5,
1521 'work_key': 't',
1522 'precision_key': 'e_global_rel',
1523 'plotting': comm_world.rank == 0,
1524 }
1526 for mode in [
1527 'RK_comp',
1528 'parallel_efficiency',
1529 'compare_strategies',
1530 ]:
1531 all_problems(**all_params, mode=mode)
1532 comm_world.Barrier()
1534 if comm_world.rank == 0:
1535 plt.show()