Coverage for pySDC/projects/Resilience/work_precision.py: 0%
488 statements
« prev ^ index » next coverage.py v7.6.12, created at 2025-03-04 07:15 +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
15from pySDC.projects.Resilience.RBC import run_RBC
16from pySDC.projects.Resilience.GS import run_GS
18from pySDC.helpers.stats_helper import get_sorted, filter_stats
19from pySDC.helpers.plot_helper import setup_mpl, figsize_by_journal
21setup_mpl(reset=True)
22LOGGER_LEVEL = 25
23LOG_TO_FILE = False
25logging.getLogger('matplotlib.texmanager').setLevel(90)
27Tends = {'run_RBC': 16.0, 'run_Lorenz': 2.0}
28t0s = {
29 'run_RBC': 10.0,
30}
33def std_log(x):
34 return np.std(np.log(x))
37MAPPINGS = {
38 'e_global': ('e_global_post_run', max, False),
39 'e_global_rel': ('e_global_rel_post_run', max, False),
40 't': ('timing_run', max, False),
41 # 'e_local_max': ('e_local_post_step', max, False),
42 'k_SDC': ('k', sum, None),
43 'k_SDC_no_restart': ('k', sum, False),
44 'k_Newton': ('work_newton', sum, None),
45 'k_linear': ('work_linear', sum, None),
46 'k_Newton_no_restart': ('work_newton', sum, False),
47 'k_rhs': ('work_rhs', sum, None),
48 'k_factorizations': ('work_factorizations', sum, None),
49 'num_steps': ('dt', len, None),
50 'restart': ('restart', sum, None),
51 'dt_mean': ('dt', np.mean, False),
52 'dt_max': ('dt', max, False),
53 'dt_min': ('dt', min, False),
54 'dt_sigma': ('dt', std_log, False),
55 'e_embedded_max': ('error_embedded_estimate', max, False),
56 'u0_increment_max': ('u0_increment', max, None),
57 'u0_increment_mean': ('u0_increment', np.mean, None),
58 'u0_increment_max_no_restart': ('u0_increment', max, False),
59 'u0_increment_mean_no_restart': ('u0_increment', np.mean, False),
60}
62logger = logging.getLogger('WorkPrecision')
63logger.setLevel(LOGGER_LEVEL)
66def get_forbidden_combinations(problem, strategy, **kwargs):
67 """
68 Check if the combination of strategy and problem is forbidden
70 Args:
71 problem (function): A problem to run
72 strategy (Strategy): SDC strategy
73 """
74 if strategy.name == 'ERK':
75 if problem.__name__ in ['run_quench', 'run_Schroedinger', 'run_AC']:
76 return True
78 return False
81def single_run(
82 problem,
83 strategy,
84 data,
85 custom_description,
86 num_procs=1,
87 comm_world=None,
88 problem_args=None,
89 hooks=None,
90 Tend=None,
91 num_procs_sweeper=1,
92):
93 """
94 Make a single run of a particular problem with a certain strategy.
96 Args:
97 problem (function): A problem to run
98 strategy (Strategy): SDC strategy
99 data (dict): Put the results in here
100 custom_description (dict): Overwrite presets
101 num_procs (int): Number of processes for the time communicator
102 comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script
103 hooks (list): List of additional hooks
104 num_procs_sweeper (int): Number of processes for the sweeper
106 Returns:
107 dict: Stats generated by the run
108 """
109 from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun
110 from pySDC.implementations.hooks.log_work import LogWork
111 from pySDC.projects.Resilience.hook import LogData
113 hooks = hooks if hooks else []
115 t_last = perf_counter()
117 num_procs_tot = num_procs * num_procs_sweeper
118 comm = comm_world.Split(comm_world.rank < num_procs_tot)
119 if comm_world.rank >= num_procs_tot:
120 comm.Free()
121 return None
123 # make communicators for time and sweepers
124 comm_time = comm.Split(comm.rank // num_procs)
125 comm_sweep = comm.Split(comm_time.rank)
127 if comm_time.size < num_procs:
128 raise Exception(f'Need at least {num_procs*num_procs_sweeper} processes, got only {comm.size}')
130 strategy_description = strategy.get_custom_description(problem, num_procs)
131 description = merge_descriptions(strategy_description, custom_description)
132 if comm_sweep.size > 1:
133 description['sweeper_params']['comm'] = comm_sweep
135 controller_params = {
136 'logger_level': LOGGER_LEVEL,
137 'log_to_file': LOG_TO_FILE,
138 'fname': 'out.txt',
139 **strategy.get_controller_params(),
140 }
141 problem_args = {} if problem_args is None else problem_args
143 Tend = Tends.get(problem.__name__, None) if Tend is None else Tend
144 t0 = t0s.get(problem.__name__, None)
146 stats, controller, crash = problem(
147 custom_description=description,
148 Tend=strategy.get_Tend(problem, num_procs) if Tend is None else Tend,
149 hook_class=[LogData, LogWork, LogGlobalErrorPostRun] + hooks,
150 custom_controller_params=controller_params,
151 use_MPI=True,
152 t0=t0,
153 comm=comm_time,
154 **problem_args,
155 )
157 t_now = perf_counter()
158 logger.debug(f'Finished run in {t_now - t_last:.2e} s')
159 t_last = perf_counter()
161 # record all the metrics
162 if comm_sweep.size > 1:
163 try:
164 stats_all = filter_stats(stats, comm=comm_sweep)
165 except MPI.Exception:
166 for key in MAPPINGS.keys():
167 data[key] += [np.nan]
168 return stats
170 else:
171 stats_all = stats
172 comm_sweep.Free()
174 for key, mapping in MAPPINGS.items():
175 if crash:
176 data[key] += [np.nan]
177 continue
178 me = get_sorted(stats_all, comm=comm_time, type=mapping[0], recomputed=mapping[2])
179 if len(me) == 0:
180 data[key] += [np.nan]
181 else:
182 data[key] += [mapping[1]([you[1] for you in me])]
184 t_now = perf_counter()
185 logger.debug(f'Recorded all data after {t_now - t_last:.2e} s')
186 t_last = perf_counter()
188 comm_time.Free()
189 comm.Free()
190 return stats
193def get_parameter(dictionary, where):
194 """
195 Get a parameter at a certain position in a dictionary of dictionaries.
197 Args:
198 dictionary (dict): The dictionary
199 where (list): The list of keys leading to the value you want
201 Returns:
202 The value of the dictionary
203 """
204 if len(where) == 1:
205 return dictionary[where[0]]
206 else:
207 return get_parameter(dictionary[where[0]], where[1:])
210def set_parameter(dictionary, where, parameter):
211 """
212 Set a parameter at a certain position in a dictionary of dictionaries
214 Args:
215 dictionary (dict): The dictionary
216 where (list): The list of keys leading to the value you want to set
217 parameter: Whatever you want to set the parameter to
219 Returns:
220 None
221 """
222 if len(where) == 1:
223 dictionary[where[0]] = parameter
224 else:
225 set_parameter(dictionary[where[0]], where[1:], parameter)
228def get_path(problem, strategy, num_procs, handle='', base_path='data/work_precision', num_procs_sweeper=1, mode=''):
229 """
230 Get the path to a certain data.
232 Args:
233 problem (function): A problem to run
234 strategy (Strategy): SDC strategy
235 num_procs (int): Number of processes for the time communicator
236 handle (str): The name of the configuration
237 base_path (str): Some path where all the files are stored
238 num_procs_sweeper (int): Number of processes for the sweeper
239 mode (str): The mode this was generated for
241 Returns:
242 str: The path to the data you are looking for
243 """
244 return f'{base_path}/{problem.__name__}-{strategy.__class__.__name__}-{handle}{"-wp" if handle else "wp"}-{num_procs}-{num_procs_sweeper}procs-{mode}.pickle'
247def record_work_precision(
248 problem,
249 strategy,
250 num_procs=1,
251 custom_description=None,
252 handle='',
253 runs=1,
254 comm_world=None,
255 problem_args=None,
256 param_range=None,
257 Tend=None,
258 hooks=None,
259 num_procs_sweeper=1,
260 mode='',
261):
262 """
263 Run problem with strategy and record the cost parameters.
265 Args:
266 problem (function): A problem to run
267 strategy (Strategy): SDC strategy
268 num_procs (int): Number of processes for the time communicator
269 custom_description (dict): Overwrite presets
270 handle (str): The name of the configuration
271 runs (int): Number of runs you want to do
272 comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script
273 num_procs_sweeper (int): Number of processes for the sweeper
275 Returns:
276 None
277 """
278 if get_forbidden_combinations(problem, strategy):
279 return None
281 data = {}
283 # prepare precision parameters
284 param = strategy.precision_parameter
285 description = merge_descriptions(
286 strategy.get_custom_description(problem, num_procs),
287 {} if custom_description is None else custom_description,
288 )
289 if param == 'e_tol':
290 power = 10.0
291 set_parameter(description, strategy.precision_parameter_loc[:-1] + ['dt_min'], 0)
292 exponents = [-3, -2, -1, 0, 1, 2, 3][::-1]
293 if problem.__name__ == 'run_vdp':
294 if type(strategy).__name__ in ["AdaptivityPolynomialError"]:
295 exponents = [0, 1, 2, 3, 5][::-1]
296 else:
297 exponents = [-3, -2, -1, 0, 0.2, 0.8, 1][::-1]
298 if problem.__name__ == 'run_RBC':
299 exponents = [1, 0, -0.5, -1, -2]
300 if problem.__name__ == 'run_GS':
301 exponents = [-2, -1, 0, 1, 2, 3][::-1]
302 if problem.__name__ == 'run_Lorenz':
303 exponents = [0, 1, 2, 3][::-1]
304 if type(strategy).__name__ in ["AdaptivityStrategy"]:
305 exponents = [0, 1, 2, 3, 4, 5][::-1]
306 elif param == 'dt':
307 power = 2.0
308 exponents = [-1, 0, 1, 2, 3][::-1]
309 elif param == 'restol':
310 power = 10.0
311 exponents = [-2, -1, 0, 1, 2, 3]
312 if problem.__name__ == 'run_vdp':
313 exponents = [-4, -3, -2, -1, 0, 1]
314 elif param == 'cfl':
315 power = 2
316 exponents = [-3, -2, -1, 0, 1]
317 else:
318 raise NotImplementedError(f"I don't know how to get default value for parameter \"{param}\"")
320 where = strategy.precision_parameter_loc
321 default = get_parameter(description, where)
322 param_range = [default * power**i for i in exponents] if param_range is None else param_range
324 if problem.__name__ == 'run_quench':
325 if param == 'restol':
326 param_range = [1e-5, 1e-6, 1e-7, 1e-8, 1e-9]
327 elif param == 'dt':
328 param_range = [1.25, 2.5, 5.0, 10.0, 20.0][::-1]
329 if problem.__name__ == 'run_RBC':
330 if param == 'dt':
331 param_range = [8e-2, 6e-2, 4e-2, 3e-2, 2e-2]
332 if problem.__name__ == 'run_GS':
333 if param == 'dt':
334 param_range = [2, 1, 0.5, 0.1]
335 if problem.__name__ == 'run_Lorenz':
336 if param == 'dt':
337 param_range = [5e-2, 2e-2, 1e-2, 5e-3]
339 # run multiple times with different parameters
340 for i in range(len(param_range)):
341 set_parameter(description, where, param_range[i])
343 data[param_range[i]] = {key: [] for key in MAPPINGS.keys()}
344 data[param_range[i]]['param'] = [param_range[i]]
345 data[param_range[i]][param] = [param_range[i]]
347 description = merge_descriptions(
348 descA=description, descB=strategy.get_description_for_tolerance(problem=problem, param=param_range[i])
349 )
350 for _j in range(runs):
351 if comm_world.rank == 0:
352 logger.log(
353 24,
354 f'Starting: {problem.__name__}: {strategy} {handle} {num_procs}-{num_procs_sweeper} procs, {param}={param_range[i]:.2e}',
355 )
356 single_run(
357 problem,
358 strategy,
359 data[param_range[i]],
360 custom_description=description,
361 comm_world=comm_world,
362 problem_args=problem_args,
363 num_procs=num_procs,
364 hooks=hooks,
365 Tend=Tend,
366 num_procs_sweeper=num_procs_sweeper,
367 )
369 comm_world.Barrier()
371 if comm_world.rank == 0:
372 if np.isfinite(data[param_range[i]]["k_linear"][-1]):
373 k_type = "k_linear"
374 elif np.isfinite(data[param_range[i]]["k_Newton"][-1]):
375 k_type = 'k_Newton'
376 else:
377 k_type = "k_SDC"
378 logger.log(
379 25,
380 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]}',
381 )
383 if comm_world.rank == 0:
384 import socket
385 import time
387 data['meta'] = {
388 'hostname': socket.gethostname(),
389 'time': time.time,
390 'runs': runs,
391 }
392 path = get_path(problem, strategy, num_procs, handle, num_procs_sweeper=num_procs_sweeper, mode=mode)
393 with open(path, 'wb') as f:
394 logger.debug(f'Dumping file \"{path}\"')
395 pickle.dump(data, f)
396 return data
399def load(**kwargs):
400 """
401 Load stored data. Arguments are passed on to `get_path`
403 Returns:
404 dict: The data
405 """
406 path = get_path(**kwargs)
407 with open(path, 'rb') as f:
408 logger.debug(f'Loading file \"{path}\"')
409 data = pickle.load(f)
410 return data
413def extract_data(data, work_key, precision_key):
414 """
415 Get the work and precision from a data object.
417 Args:
418 data (dict): Data from work-precision measurements
419 work_key (str): Name of variable on x-axis
420 precision_key (str): Name of variable on y-axis
422 Returns:
423 numpy array: Work
424 numpy array: Precision
425 """
426 keys = [key for key in data.keys() if key not in ['meta']]
427 work = [np.nanmean(data[key][work_key]) for key in keys]
428 precision = [np.nanmean(data[key][precision_key]) for key in keys]
429 return np.array(work), np.array(precision)
432def get_order(work_key='e_global', precision_key='param', strategy=None, handle=None, **kwargs):
433 data = load(**kwargs, strategy=strategy, handle=handle)
434 work, precision = extract_data(data, work_key, precision_key)
436 order = [np.log(precision[i + 1] / precision[i]) / np.log(work[i + 1] / work[i]) for i in range(len(work) - 1)]
438 print(f'Order for {strategy} {handle}: {np.mean(order):.2f}')
441def plot_work_precision(
442 problem,
443 strategy,
444 num_procs,
445 ax,
446 work_key='k_SDC',
447 precision_key='e_global',
448 handle='',
449 plotting_params=None,
450 comm_world=None,
451 num_procs_sweeper=1,
452 mode='',
453): # pragma: no cover
454 """
455 Plot data from running a problem with a strategy.
457 Args:
458 problem (function): A problem to run
459 strategy (Strategy): SDC strategy
460 num_procs (int): Number of processes for the time communicator
461 ax (matplotlib.pyplot.axes): Somewhere to plot
462 work_key (str): The key in the recorded data you want on the x-axis
463 precision_key (str): The key in the recorded data you want on the y-axis
464 handle (str): The name of the configuration
465 plotting_params (dict): Will be passed when plotting
466 comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script
467 num_procs_sweeper (int): Number of processes for the sweeper
468 mode (str): The of the configurations you want to retrieve
470 Returns:
471 None
472 """
473 if comm_world.rank > 0 or get_forbidden_combinations(problem, strategy):
474 return None
476 data = load(
477 problem=problem,
478 strategy=strategy,
479 num_procs=num_procs,
480 handle=handle,
481 num_procs_sweeper=num_procs_sweeper,
482 mode=mode,
483 )
485 work, precision = extract_data(data, work_key, precision_key)
486 keys = [key for key in data.keys() if key not in ['meta']]
488 for key in [work_key, precision_key]:
489 rel_variance = [np.std(data[me][key]) / max([np.nanmean(data[me][key]), 1.0]) for me in keys]
490 if not all(me < 1e-1 or not np.isfinite(me) for me in rel_variance):
491 logger.warning(
492 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}"
493 )
495 style = merge_descriptions(
496 {**strategy.style, 'label': f'{strategy.style["label"]}{f" {handle}" if handle else ""}'},
497 plotting_params if plotting_params else {},
498 )
500 mask = np.logical_and(np.isfinite(work), np.isfinite(precision))
501 ax.loglog(work[mask], precision[mask], **style)
503 # get_order(
504 # problem=problem,
505 # strategy=strategy,
506 # num_procs=num_procs,
507 # handle=handle,
508 # work_key=work_key,
509 # precision_key=precision_key,
510 # )
512 if 't' in [work_key, precision_key]:
513 meta = data.get('meta', {})
515 if meta.get('hostname', None) in ['thomas-work']:
516 ax.text(0.1, 0.1, "Laptop timings!", transform=ax.transAxes)
517 if meta.get('runs', None) == 1:
518 ax.text(0.1, 0.2, "No sampling!", transform=ax.transAxes)
520 if problem.__name__ == 'run_vdp':
521 if mode == 'parallel_efficiency':
522 # ax.set_xticks([6e-1, 2e0])
523 ax.set_xticks(
524 ticks=[
525 0.4,
526 5e-1,
527 6e-1,
528 7e-1,
529 8e-1,
530 9e-1,
531 2e0,
532 ],
533 labels=['']
534 + [r'$5\times 10^{-1}$']
535 + [
536 '',
537 ]
538 * 4
539 + [r'$2\times 10^0$'],
540 minor=True,
541 )
542 elif mode == 'RK_comp':
543 ax.set_xticks(
544 ticks=[
545 5e-1,
546 6e-1,
547 7e-1,
548 8e-1,
549 9e-1,
550 2e0,
551 ],
552 labels=[r'$5\times 10^{-1}$']
553 + [
554 '',
555 ]
556 * 4
557 + [r'$2\times 10^0$'],
558 minor=True,
559 )
560 elif problem.__name__ == 'run_quench':
561 if mode == 'RK_comp':
562 ax.set_xticks(
563 ticks=[
564 0.2,
565 0.3,
566 0.4,
567 5e-1,
568 6e-1,
569 7e-1,
570 8e-1,
571 9e-1,
572 2e0,
573 ],
574 labels=['']
575 + [r'$3\times 10^{-1}$']
576 + [
577 '',
578 ]
579 * 7,
580 minor=True,
581 )
584def plot_parallel_efficiency_diagonalSDC(
585 ax, work_key, precision_key, num_procs_sweeper, num_procs=1, **kwargs
586): # pragma: no cover
587 serial_data = load(
588 num_procs=num_procs,
589 num_procs_sweeper=1,
590 **kwargs,
591 )
592 parallel_data = load(
593 num_procs=num_procs,
594 num_procs_sweeper=num_procs_sweeper,
595 **kwargs,
596 )
597 serial_work, serial_precision = extract_data(serial_data, work_key, precision_key)
598 parallel_work, parallel_precision = extract_data(parallel_data, work_key, precision_key)
599 # assert np.allclose(serial_precision, parallel_precision)
601 serial_work = np.asarray(serial_work)
602 parallel_work = np.asarray(parallel_work)
604 # ax.loglog(serial_work, serial_precision)
605 # ax.loglog(parallel_work, parallel_precision)
607 speedup = serial_work / parallel_work
608 parallel_efficiency = np.median(speedup) / num_procs_sweeper
609 ax.plot(serial_precision, speedup)
610 ax.set_xscale('log')
611 ax.set_ylabel('speedup')
613 if 't' in [work_key, precision_key]:
614 meta = parallel_data.get('meta', {})
616 if meta.get('hostname', None) in ['thomas-work']:
617 ax.text(0.1, 0.1, "Laptop timings!", transform=ax.transAxes)
618 if meta.get('runs', None) == 1:
619 ax.text(0.1, 0.2, "No sampling!", transform=ax.transAxes)
621 return np.median(speedup), parallel_efficiency
624def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only=False): # pragma: no cover
625 """
626 Decorate a plot
628 Args:
629 ax (matplotlib.pyplot.axes): Somewhere to plot
630 problem (function): A problem to run
631 work_key (str): The key in the recorded data you want on the x-axis
632 precision_key (str): The key in the recorded data you want on the y-axis
633 num_procs (int): Number of processes for the time communicator
634 title_only (bool): Put only the title on top, or do the whole shebang
636 Returns:
637 None
638 """
639 labels = {
640 'k_SDC': 'SDC iterations',
641 'k_SDC_no_restart': 'SDC iterations (restarts excluded)',
642 'k_Newton': 'Newton iterations',
643 'k_Newton_no_restart': 'Newton iterations (restarts excluded)',
644 'k_rhs': 'right hand side evaluations',
645 'k_factorizations': 'matrix factorizations',
646 't': 'wall clock time / s',
647 'e_global': 'global error',
648 'e_global_rel': 'relative global error',
649 'e_local_max': 'max. local error',
650 'restart': 'restarts',
651 'dt_max': r'$\Delta t_\mathrm{max}$',
652 'dt_min': r'$\Delta t_\mathrm{min}$',
653 'dt_sigma': r'$\sigma(\Delta t)$',
654 'dt_mean': r'$\bar{\Delta t}$',
655 'param': 'accuracy parameter',
656 'u0_increment_max': r'$\| \Delta u_0 \|_{\infty} $',
657 'u0_increment_mean': r'$\bar{\Delta u_0}$',
658 'u0_increment_max_no_restart': r'$\| \Delta u_0 \|_{\infty} $ (restarts excluded)',
659 'u0_increment_mean_no_restart': r'$\bar{\Delta u_0}$ (restarts excluded)',
660 'k_linear': 'Linear solver iterations',
661 'speedup': 'Speedup',
662 'nprocs': r'$N_\mathrm{procs}$',
663 '': '',
664 }
666 if not title_only:
667 ax.set_xlabel(labels.get(work_key, 'work'))
668 ax.set_ylabel(labels.get(precision_key, 'precision'))
669 # ax.legend(frameon=False)
671 titles = {
672 'run_vdp': 'Van der Pol',
673 'run_Lorenz': 'Lorenz attractor',
674 'run_Schroedinger': r'Schr\"odinger',
675 'run_quench': 'Quench',
676 'run_AC': 'Allen-Cahn',
677 'run_RBC': 'Rayleigh-Benard',
678 'run_GS': 'Gray-Scott',
679 }
680 ax.set_title(titles.get(problem.__name__, ''))
683def execute_configurations(
684 problem,
685 configurations,
686 work_key,
687 precision_key,
688 num_procs,
689 ax,
690 decorate,
691 record,
692 runs,
693 comm_world,
694 plotting,
695 Tend=None,
696 num_procs_sweeper=1,
697 mode='',
698):
699 """
700 Run for multiple configurations.
702 Args:
703 problem (function): A problem to run
704 configurations (dict): The configurations you want to run with
705 work_key (str): The key in the recorded data you want on the x-axis
706 precision_key (str): The key in the recorded data you want on the y-axis
707 num_procs (int): Number of processes for the time communicator
708 ax (matplotlib.pyplot.axes): Somewhere to plot
709 decorate (bool): Whether to decorate fully or only put the title
710 record (bool): Whether to only plot or also record the data first
711 runs (int): Number of runs you want to do
712 comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script
713 plotting (bool): Whether to plot something
714 num_procs_sweeper (int): Number of processes for the sweeper
715 mode (str): What you want to look at
717 Returns:
718 None
719 """
720 for _, config in configurations.items():
721 for strategy in config['strategies']:
722 shared_args = {
723 'problem': problem,
724 'strategy': strategy,
725 'handle': config.get('handle', ''),
726 'num_procs': config.get('num_procs', num_procs),
727 'num_procs_sweeper': config.get('num_procs_sweeper', num_procs_sweeper),
728 }
729 if record:
730 logger.debug('Recording work precision')
731 record_work_precision(
732 **shared_args,
733 custom_description=config.get('custom_description', {}),
734 runs=runs,
735 comm_world=comm_world,
736 problem_args=config.get('problem_args', {}),
737 param_range=config.get('param_range', None),
738 hooks=config.get('hooks', None),
739 Tend=config.get('Tend') if Tend is None else Tend,
740 mode=mode,
741 )
742 if plotting and comm_world.rank == 0:
743 logger.debug('Plotting')
744 plot_work_precision(
745 **shared_args,
746 work_key=work_key,
747 precision_key=precision_key,
748 ax=ax,
749 plotting_params=config.get('plotting_params', {}),
750 comm_world=comm_world,
751 mode=mode,
752 )
754 if comm_world.rank == 0:
755 decorate_panel(
756 ax=ax,
757 problem=problem,
758 work_key=work_key,
759 precision_key=precision_key,
760 num_procs=num_procs,
761 title_only=not decorate,
762 )
765def get_configs(mode, problem):
766 """
767 Get configurations for work-precision plots. These are dictionaries containing strategies and handles and so on.
769 Args:
770 mode (str): The of the configurations you want to retrieve
771 problem (function): A problem to run
773 Returns:
774 dict: Configurations
775 """
776 configurations = {}
777 if mode == 'regular':
778 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, BaseStrategy, IterateStrategy
780 handle = 'regular'
781 configurations[0] = {
782 'handle': handle,
783 'strategies': [AdaptivityStrategy(useMPI=True), BaseStrategy(useMPI=True), IterateStrategy(useMPI=True)],
784 }
785 elif mode == 'step_size_limiting':
786 from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter
787 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, ESDIRKStrategy
789 limits = [
790 25.0,
791 50.0,
792 ]
793 colors = ['teal', 'magenta']
794 markers = ['v', 'x']
795 markersize = 9
796 for i in range(len(limits)):
797 configurations[i] = {
798 'custom_description': {'convergence_controllers': {StepSizeLimiter: {'dt_max': limits[i]}}},
799 'handle': f'steplimiter{limits[i]:.0f}',
800 'strategies': [AdaptivityStrategy(useMPI=True)],
801 'plotting_params': {
802 'color': colors[i],
803 'marker': markers[i],
804 'label': rf'$\Delta t \leq { {limits[i]:.0f}} $',
805 # 'ls': '',
806 'markersize': markersize,
807 },
808 'num_procs': 1,
809 }
810 configurations[99] = {
811 'custom_description': {},
812 'handle': 'no limits',
813 'plotting_params': {
814 'label': 'no limiter',
815 # 'ls': '',
816 'markersize': markersize,
817 },
818 'strategies': [AdaptivityStrategy(useMPI=True)],
819 'num_procs': 1,
820 }
821 elif mode == 'dynamic_restarts':
822 """
823 Compare Block Gauss-Seidel SDC with restarting the first step in the block or the first step that exceeded the error threshold.
824 """
825 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, AdaptivityRestartFirstStep
827 desc = {}
828 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'}
829 desc['step_params'] = {'maxiter': 5}
831 ls = {
832 1: '-',
833 2: '--',
834 3: '-.',
835 4: ':',
836 5: ':',
837 }
839 configurations[-1] = {
840 'strategies': [AdaptivityStrategy(useMPI=True)],
841 'num_procs': 1,
842 }
844 for num_procs in [4, 2]:
845 plotting_params = {'ls': ls[num_procs], 'label': f'adaptivity {num_procs} procs'}
846 configurations[num_procs] = {
847 'strategies': [AdaptivityStrategy(useMPI=True), AdaptivityRestartFirstStep(useMPI=True)],
848 'custom_description': desc,
849 'num_procs': num_procs,
850 'plotting_params': plotting_params,
851 }
853 elif mode == 'compare_strategies':
854 """
855 Compare the different SDC strategies.
856 """
857 from pySDC.projects.Resilience.strategies import (
858 AdaptivityStrategy,
859 kAdaptivityStrategy,
860 AdaptivityPolynomialError,
861 BaseStrategy,
862 )
864 newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
866 configurations[1] = {
867 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
868 }
869 configurations[2] = {
870 'strategies': [kAdaptivityStrategy(useMPI=True)],
871 }
872 configurations[0] = {
873 'custom_description': {
874 'step_params': {'maxiter': 5},
875 'sweeper_params': {'num_nodes': 3, 'quad_type': 'RADAU-RIGHT'},
876 },
877 'strategies': [
878 BaseStrategy(useMPI=True),
879 AdaptivityStrategy(useMPI=True),
880 ],
881 }
883 elif mode == 'RK_comp':
884 """
885 Compare parallel adaptive SDC to Runge-Kutta
886 """
887 from pySDC.projects.Resilience.strategies import (
888 AdaptivityStrategy,
889 ERKStrategy,
890 ESDIRKStrategy,
891 ARKStrategy,
892 AdaptivityPolynomialError,
893 ARK3_CFL_Strategy,
894 )
896 if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_RBC', 'run_GS']:
897 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
898 else:
899 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
900 generic_implicit_MPI as parallel_sweeper,
901 )
903 newton_inexactness = False if problem.__name__ in ['run_vdp', 'run_RBC', 'run_GS'] else True
905 desc = {}
906 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': "EE"}
907 desc['step_params'] = {'maxiter': 5}
908 num_procs_dt = {
909 'run_RBC': 1,
910 }.get(problem.__name__, 4)
912 desc_poly = {}
913 desc_poly['sweeper_class'] = parallel_sweeper
914 num_procs_dt_k = 3
916 ls = {
917 1: '--',
918 2: '--',
919 3: '-',
920 4: '-',
921 5: '-',
922 12: ':',
923 }
924 RK_strategies = []
925 if problem.__name__ in ['run_Lorenz']:
926 RK_strategies.append(ERKStrategy(useMPI=True))
927 desc_poly['sweeper_params'] = {'QI': 'MIN-SR-S', 'QE': 'PIC'}
928 desc['sweeper_params']['QI'] = 'MIN-SR-S'
929 desc['sweeper_params']['QE'] = 'PIC'
930 if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS']:
931 RK_strategies.append(ARKStrategy(useMPI=True))
932 elif problem.__name__ == 'run_RBC':
933 RK_strategies.append(ARK3_CFL_Strategy(useMPI=True))
934 desc['sweeper_params']['num_nodes'] = 2
935 desc['sweeper_params']['QI'] = 'LU'
936 desc['sweeper_params']['QE'] = 'PIC'
937 desc['step_params']['maxiter'] = 3
939 desc_poly['sweeper_params'] = {'num_nodes': 2, 'QI': 'MIN-SR-S'}
940 num_procs_dt_k = 2
941 else:
942 RK_strategies.append(ESDIRKStrategy(useMPI=True))
944 configurations[-1] = {
945 'strategies': RK_strategies,
946 'num_procs': 1,
947 }
948 if problem.__name__ == 'run_Lorenz':
949 configurations[3] = {
950 'custom_description': desc_poly,
951 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
952 'num_procs': 4,
953 'num_procs_sweeper': num_procs_dt_k,
954 'plotting_params': {
955 'label': rf'$\Delta t$-$k$-adaptivity $N$=4x{num_procs_dt_k}',
956 'ls': ls[num_procs_dt_k * 4],
957 },
958 }
959 else:
960 configurations[3] = {
961 'custom_description': desc_poly,
962 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
963 'num_procs': 1,
964 'num_procs_sweeper': num_procs_dt_k,
965 'plotting_params': {
966 'label': rf'$\Delta t$-$k$-adaptivity $N$=1x{num_procs_dt_k}',
967 'ls': ls[num_procs_dt_k],
968 },
969 }
970 if problem.__name__ in ['run_Lorenz']:
971 configurations[2] = {
972 'strategies': [AdaptivityStrategy(useMPI=True)],
973 'custom_description': {**desc, 'sweeper_class': parallel_sweeper},
974 'num_procs': num_procs_dt,
975 'num_procs_sweeper': num_procs_dt_k,
976 'plotting_params': {
977 'label': rf'$\Delta t$-adaptivity $N$={num_procs_dt}x3',
978 'ls': ls[num_procs_dt * num_procs_dt_k],
979 },
980 }
981 else:
982 configurations[2] = {
983 'strategies': [AdaptivityStrategy(useMPI=True)],
984 'custom_description': desc,
985 'num_procs': num_procs_dt,
986 'plotting_params': {'label': rf'$\Delta t$-adaptivity $N$={num_procs_dt}x1', 'ls': ls[num_procs_dt]},
987 }
989 elif mode == 'RK_comp_high_order_RBC':
990 """
991 Compare parallel adaptive SDC to Runge-Kutta at order five for RBC problem
992 """
993 from pySDC.projects.Resilience.strategies import (
994 AdaptivityStrategy,
995 ERKStrategy,
996 ESDIRKStrategy,
997 ARKStrategy,
998 AdaptivityPolynomialError,
999 ARK3_CFL_Strategy,
1000 )
1002 assert problem.__name__ == 'run_RBC'
1004 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1006 newton_inexactness = False
1008 desc = {}
1009 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': "EE"}
1010 desc['step_params'] = {'maxiter': 5}
1011 num_procs_dt = 1
1013 desc_poly = {}
1014 desc_poly['sweeper_class'] = parallel_sweeper
1015 num_procs_dt_k = 3
1017 ls = {
1018 1: '--',
1019 2: '--',
1020 3: '-',
1021 4: '-',
1022 5: '-',
1023 12: ':',
1024 }
1025 RK_strategies = [ARK3_CFL_Strategy(useMPI=True)]
1026 desc['sweeper_params']['num_nodes'] = 3
1027 desc['sweeper_params']['QI'] = 'LU'
1028 desc['sweeper_params']['QE'] = 'PIC'
1029 desc['step_params']['maxiter'] = 5
1031 desc_poly['sweeper_params'] = {'num_nodes': 3, 'QI': 'MIN-SR-S'}
1032 num_procs_dt_k = 3
1034 configurations[-1] = {
1035 'strategies': RK_strategies,
1036 'num_procs': 1,
1037 }
1038 configurations[3] = {
1039 'custom_description': desc_poly,
1040 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
1041 'num_procs': 1,
1042 'num_procs_sweeper': num_procs_dt_k,
1043 'plotting_params': {
1044 'label': rf'$\Delta t$-$k$-adaptivity $N$=1x{num_procs_dt_k}',
1045 'ls': ls[num_procs_dt_k],
1046 },
1047 }
1048 configurations[2] = {
1049 'strategies': [AdaptivityStrategy(useMPI=True)],
1050 'custom_description': desc,
1051 'num_procs': num_procs_dt,
1052 'plotting_params': {'label': rf'$\Delta t$-adaptivity $N$={num_procs_dt}x1', 'ls': ls[num_procs_dt]},
1053 }
1055 elif mode == 'parallel_efficiency':
1056 """
1057 Compare parallel runs of the step size adaptive SDC
1058 """
1059 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, AdaptivityPolynomialError
1061 if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS', 'run_RBC']:
1062 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1063 else:
1064 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1065 generic_implicit_MPI as parallel_sweeper,
1066 )
1068 desc = {}
1069 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': 'EE'}
1070 desc['step_params'] = {'maxiter': 5}
1072 if problem.__name__ in ['run_RBC']:
1073 desc['sweeper_params']['QE'] = 'PIC'
1074 desc['sweeper_params']['QI'] = 'LU'
1076 ls = {
1077 1: '-',
1078 2: '--',
1079 3: '-.',
1080 4: '--',
1081 5: ':',
1082 12: ':',
1083 }
1085 newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
1087 for num_procs in [4, 1]:
1088 plotting_params = (
1089 {'ls': ls[num_procs], 'label': fr'$\Delta t$-adaptivity $N$={num_procs}x1'} if num_procs > 1 else {}
1090 )
1091 configurations[num_procs] = {
1092 'strategies': [AdaptivityStrategy(useMPI=True)],
1093 'custom_description': desc.copy(),
1094 'num_procs': num_procs,
1095 'plotting_params': plotting_params.copy(),
1096 }
1097 configurations[num_procs * 100 + 79] = {
1098 'custom_description': {'sweeper_class': parallel_sweeper},
1099 'strategies': [
1100 AdaptivityPolynomialError(
1101 useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True
1102 )
1103 ],
1104 'num_procs_sweeper': 3,
1105 'num_procs': num_procs,
1106 'plotting_params': {
1107 'ls': ls.get(num_procs * 3, '-'),
1108 'label': rf'$\Delta t$-$k$-adaptivity $N$={num_procs}x3',
1109 },
1110 }
1112 configurations[200 + 79] = {
1113 'strategies': [
1114 AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True)
1115 ],
1116 'num_procs': 1,
1117 }
1118 elif mode == 'parallel_efficiency_dt':
1119 """
1120 Compare parallel runs of the step size adaptive SDC
1121 """
1122 from pySDC.projects.Resilience.strategies import AdaptivityStrategy
1124 if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS', 'run_RBC']:
1125 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1126 else:
1127 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1128 generic_implicit_MPI as parallel_sweeper,
1129 )
1131 desc = {}
1132 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': 'EE'}
1133 desc['step_params'] = {'maxiter': 5}
1135 if problem.__name__ in ['run_RBC']:
1136 desc['sweeper_params']['QE'] = 'PIC'
1137 desc['sweeper_params']['QI'] = 'LU'
1139 desc_serial = {
1140 'step_params': {'maxiter': 5},
1141 'sweeper_params': {'num_nodes': 3, 'quad_type': 'RADAU-RIGHT'},
1142 }
1144 ls = {
1145 1: '-',
1146 2: '--',
1147 3: '-.',
1148 4: '--',
1149 5: ':',
1150 12: ':',
1151 }
1153 newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
1155 for num_procs in [4, 1]:
1156 configurations[num_procs] = {
1157 'strategies': [AdaptivityStrategy(useMPI=True)],
1158 'custom_description': desc.copy() if num_procs > 1 else desc_serial,
1159 'num_procs': num_procs,
1160 'plotting_params': {
1161 'ls': ls.get(num_procs, '-'),
1162 'label': rf'$\Delta t$-adaptivity $N$={num_procs}x1',
1163 },
1164 }
1165 configurations[num_procs * 200 + 79] = {
1166 'custom_description': {
1167 'sweeper_class': parallel_sweeper,
1168 'sweeper_params': {'QI': 'MIN-SR-S', 'QE': 'PIC'},
1169 'step_params': {'maxiter': 5},
1170 },
1171 'strategies': [AdaptivityStrategy(useMPI=True)],
1172 'num_procs_sweeper': 3,
1173 'num_procs': num_procs,
1174 'plotting_params': {
1175 'ls': ls.get(num_procs * 3, '-'),
1176 'label': rf'$\Delta t$-adaptivity $N$={num_procs}x3',
1177 },
1178 }
1179 elif mode == 'parallel_efficiency_dt_k':
1180 """
1181 Compare parallel runs of the step size adaptive SDC
1182 """
1183 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
1185 if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS', 'run_RBC']:
1186 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1187 else:
1188 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1189 generic_implicit_MPI as parallel_sweeper,
1190 )
1192 ls = {
1193 1: '-',
1194 2: '--',
1195 3: '-.',
1196 4: '--',
1197 5: ':',
1198 12: ':',
1199 }
1201 QI = {
1202 (1, 3, 'run_Lorenz'): 'MIN-SR-NS',
1203 (1, 1, 'run_Lorenz'): 'MIN-SR-NS',
1204 (4, 1, 'run_Lorenz'): 'IE',
1205 }
1207 newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
1209 for num_procs in [4, 1]:
1210 configurations[num_procs * 100 + 79] = {
1211 'custom_description': {
1212 'sweeper_class': parallel_sweeper,
1213 'sweeper_params': {'QI': QI.get((num_procs, 3, problem.__name__), 'MIN-SR-S'), 'QE': 'PIC'},
1214 },
1215 'strategies': [
1216 AdaptivityPolynomialError(
1217 useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True
1218 )
1219 ],
1220 'num_procs_sweeper': 3,
1221 'num_procs': num_procs,
1222 'plotting_params': {
1223 'ls': ls.get(num_procs * 3, '-'),
1224 'label': rf'$\Delta t$-$k$-adaptivity $N$={num_procs}x3',
1225 },
1226 }
1227 configurations[num_procs * 200 + 79] = {
1228 'custom_description': {
1229 'sweeper_params': {'QI': QI.get((num_procs, 1, problem.__name__), 'MIN-SR-S'), 'QE': 'PIC'},
1230 },
1231 'strategies': [
1232 AdaptivityPolynomialError(
1233 useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True
1234 )
1235 ],
1236 'num_procs_sweeper': 1,
1237 'num_procs': num_procs,
1238 'plotting_params': {
1239 'ls': ls.get(num_procs, '-'),
1240 'label': rf'$\Delta t$-$k$-adaptivity $N$={num_procs}x1',
1241 },
1242 }
1243 elif mode == 'interpolate_between_restarts':
1244 """
1245 Compare adaptivity with interpolation between restarts and without
1246 """
1247 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
1249 i = 0
1250 for interpolate_between_restarts, handle, ls in zip(
1251 [True, False], ['Interpolation between restarts', 'regular'], ['--', '-']
1252 ):
1253 configurations[i] = {
1254 'strategies': [
1255 AdaptivityPolynomialError(interpolate_between_restarts=interpolate_between_restarts, useMPI=True)
1256 ],
1257 'plotting_params': {'ls': ls},
1258 'handle': handle,
1259 }
1260 i += 1
1261 elif mode == 'diagonal_SDC':
1262 """
1263 Run diagonal SDC with different number of nodes and ranks. You can use this to compute a speedup, but it's not strong scaling!
1264 """
1265 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
1267 if problem.__name__ in ['run_Schroedinger']:
1268 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1269 else:
1270 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1271 generic_implicit_MPI as parallel_sweeper,
1272 )
1274 for parallel in [False, True]:
1275 desc = {'sweeper_class': parallel_sweeper} if parallel else {}
1276 for num_nodes, ls in zip([3, 4, 2], ['-', '--', ':', '-.']):
1277 configurations[num_nodes + (99 if parallel else 0)] = {
1278 'custom_description': {**desc, 'sweeper_params': {'num_nodes': num_nodes}},
1279 'strategies': [
1280 AdaptivityPolynomialError(useMPI=True, newton_inexactness=True, linear_inexactness=True)
1281 ],
1282 'num_procs_sweeper': num_nodes if parallel else 1,
1283 'num_procs': 1,
1284 'handle': f'{num_nodes} nodes',
1285 'plotting_params': {
1286 'ls': ls,
1287 'label': f'{num_nodes} procs',
1288 # **{'color': 'grey' if parallel else None},
1289 },
1290 }
1292 elif mode[:13] == 'vdp_stiffness':
1293 """
1294 Run van der Pol with different parameter for the nonlinear term, which controls the stiffness.
1295 """
1296 from pySDC.projects.Resilience.strategies import (
1297 AdaptivityStrategy,
1298 ERKStrategy,
1299 ESDIRKStrategy,
1300 AdaptivityPolynomialError,
1301 )
1302 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1303 generic_implicit_MPI as parallel_sweeper,
1304 )
1306 Tends = {
1307 1000: 2000,
1308 100: 200,
1309 10: 20,
1310 0: 2,
1311 }
1312 mu = float(mode[14:])
1313 Tend = Tends[mu]
1315 problem_desc = {'problem_params': {'mu': mu}}
1317 desc = {}
1318 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'}
1319 desc['step_params'] = {'maxiter': 5}
1320 desc['problem_params'] = problem_desc['problem_params']
1322 ls = {
1323 1: '-',
1324 2: '--',
1325 3: '-.',
1326 4: ':',
1327 5: ':',
1328 'MIN-SR-S': '-',
1329 'MIN-SR-NS': '--',
1330 'MIN-SR-FLEX': '-.',
1331 }
1333 if mu < 100:
1334 configurations[2] = {
1335 'strategies': [ERKStrategy(useMPI=True)],
1336 'num_procs': 1,
1337 'handle': mode,
1338 'plotting_params': {'label': 'CP5(4)'},
1339 'custom_description': problem_desc,
1340 'Tend': Tend,
1341 }
1342 configurations[1] = {
1343 'strategies': [AdaptivityStrategy(useMPI=True)],
1344 'custom_description': desc,
1345 'num_procs': 4,
1346 'plotting_params': {'ls': ls[1], 'label': 'SDC $N$=4x1'},
1347 'handle': mode,
1348 'Tend': Tend,
1349 }
1350 configurations[4] = {
1351 'strategies': [ESDIRKStrategy(useMPI=True)],
1352 'num_procs': 1,
1353 'handle': mode,
1354 'plotting_params': {'label': 'ESDIRK5(3)'},
1355 'custom_description': problem_desc,
1356 'Tend': Tend,
1357 }
1358 for QI, i in zip(
1359 [
1360 'MIN-SR-S',
1361 # 'MIN-SR-FLEX',
1362 ],
1363 [9991, 12123127391, 1231723109247102731092],
1364 ):
1365 configurations[i] = {
1366 'custom_description': {
1367 'sweeper_params': {'num_nodes': 3, 'QI': QI},
1368 'problem_params': desc["problem_params"],
1369 'sweeper_class': parallel_sweeper,
1370 },
1371 'strategies': [
1372 AdaptivityPolynomialError(
1373 useMPI=True, newton_inexactness=False, linear_inexactness=False, max_slope=4
1374 )
1375 ],
1376 'num_procs_sweeper': 3,
1377 'num_procs': 1,
1378 'plotting_params': {
1379 'ls': ls.get(QI, '-'),
1380 'label': rf'$\Delta t$-$k$-adaptivity $N$={1}x3',
1381 },
1382 'handle': f'{mode}-{QI}',
1383 'Tend': Tend,
1384 }
1386 elif mode == 'inexactness':
1387 """
1388 Compare inexact SDC to exact SDC
1389 """
1390 from pySDC.projects.Resilience.strategies import (
1391 AdaptivityPolynomialError,
1392 )
1394 if problem.__name__ in ['run_Schroedinger']:
1395 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1396 else:
1397 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1398 generic_implicit_MPI as parallel_sweeper,
1399 )
1401 strategies = [
1402 AdaptivityPolynomialError,
1403 ]
1405 inexactness = {
1406 'newton_inexactness': True,
1407 'linear_inexactness': True,
1408 }
1409 no_inexactness = {
1410 'newton_inexactness': False,
1411 'linear_inexactness': False,
1412 'SDC_maxiter': 99,
1413 'use_restol_rel': False,
1414 }
1416 configurations[1] = {
1417 'custom_description': {'sweeper_class': parallel_sweeper},
1418 'strategies': [me(useMPI=True, **no_inexactness) for me in strategies],
1419 'num_procs_sweeper': 3,
1420 'handle': 'exact',
1421 'plotting_params': {'ls': '--'},
1422 }
1423 configurations[0] = {
1424 'custom_description': {'sweeper_class': parallel_sweeper},
1425 'strategies': [me(useMPI=True, **inexactness) for me in strategies],
1426 'handle': 'inexact',
1427 'num_procs_sweeper': 3,
1428 }
1429 elif mode == 'compare_adaptivity':
1430 """
1431 Compare various modes of adaptivity
1432 """
1433 # TODO: configurations not final!
1434 from pySDC.projects.Resilience.strategies import (
1435 # AdaptivityCollocationTypeStrategy,
1436 # AdaptivityCollocationRefinementStrategy,
1437 AdaptivityStrategy,
1438 # AdaptivityExtrapolationWithinQStrategy,
1439 ESDIRKStrategy,
1440 ARKStrategy,
1441 AdaptivityPolynomialError,
1442 )
1444 if problem.__name__ in ['run_Schroedinger']:
1445 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1446 else:
1447 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1448 generic_implicit_MPI as parallel_sweeper,
1449 )
1451 inexactness_params = {
1452 # 'double_adaptivity': True,
1453 'newton_inexactness': True,
1454 'linear_inexactness': True,
1455 }
1457 strategies = [
1458 AdaptivityPolynomialError,
1459 # AdaptivityCollocationTypeStrategy,
1460 # AdaptivityExtrapolationWithinQStrategy,
1461 ]
1463 # restol = None
1464 # for strategy in strategies:
1465 # strategy.restol = restol
1467 configurations[1] = {
1468 'custom_description': {'sweeper_class': parallel_sweeper},
1469 'strategies': [me(useMPI=True, **inexactness_params) for me in strategies],
1470 'handle': 'parallel',
1471 'num_procs_sweeper': 3,
1472 'plotting_params': {'ls': '-', 'label': '3 procs'},
1473 }
1474 configurations[2] = {
1475 'strategies': [me(useMPI=True, **inexactness_params) for me in strategies],
1476 'plotting_params': {'ls': '--'},
1477 }
1478 configurations[4] = {
1479 'custom_description': {'step_params': {'maxiter': 5}},
1480 'strategies': [AdaptivityStrategy(useMPI=True)],
1481 }
1483 desc_RK = {}
1484 configurations[-1] = {
1485 'strategies': [
1486 ARKStrategy(useMPI=True) if problem.__name__ == 'run_Schroedinger' else ESDIRKStrategy(useMPI=True),
1487 ],
1488 'num_procs': 1,
1489 'custom_description': desc_RK,
1490 }
1492 elif mode == 'preconditioners':
1493 """
1494 Compare different preconditioners
1495 """
1496 from pySDC.projects.Resilience.strategies import (
1497 AdaptivityStrategy,
1498 IterateStrategy,
1499 BaseStrategy,
1500 ESDIRKStrategy,
1501 ERKStrategy,
1502 AdaptivityPolynomialError,
1503 )
1505 inexacness = True
1506 strategies = [
1507 AdaptivityPolynomialError(
1508 useMPI=True, SDC_maxiter=29, newton_inexactness=inexacness, linear_inexactness=inexacness
1509 ),
1510 BaseStrategy(useMPI=True),
1511 ]
1513 desc = {}
1514 desc['sweeper_params'] = {
1515 'num_nodes': 3,
1516 }
1517 # desc['step_params'] = {'maxiter': 5}
1519 precons = ['IE', 'LU']
1520 ls = ['-.', '--', '-', ':']
1521 for i in range(len(precons) + 1):
1522 if i < len(precons):
1523 desc['sweeper_params']['QI'] = precons[i]
1524 handle = precons[i]
1525 else:
1526 handle = None
1527 configurations[i] = {
1528 'strategies': strategies,
1529 'custom_description': copy.deepcopy(desc),
1530 'handle': handle,
1531 'plotting_params': {'ls': ls[i]},
1532 }
1533 elif mode == 'RK_comp_high_order':
1534 """
1535 Compare higher order SDC than we can get with RKM to RKM
1536 """
1537 from pySDC.projects.Resilience.strategies import (
1538 AdaptivityStrategy,
1539 ERKStrategy,
1540 ESDIRKStrategy,
1541 ARKStrategy,
1542 AdaptivityPolynomialError,
1543 )
1545 if problem.__name__ in ['run_Schroedinger']:
1546 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1547 else:
1548 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1549 generic_implicit_MPI as parallel_sweeper,
1550 )
1552 desc = {}
1553 desc['sweeper_params'] = {'num_nodes': 4, 'QI': 'IE', 'QE': "EE"}
1554 desc['step_params'] = {'maxiter': 7}
1556 desc_poly = {}
1557 desc_poly['sweeper_class'] = parallel_sweeper
1558 desc_poly['sweeper_params'] = {'num_nodes': 4}
1560 ls = {
1561 1: '-',
1562 2: '--',
1563 3: '-.',
1564 4: ':',
1565 5: ':',
1566 }
1568 desc_RK = {}
1569 if problem.__name__ in ['run_Schroedinger']:
1570 desc_RK['problem_params'] = {'imex': True}
1572 configurations[3] = {
1573 'custom_description': desc_poly,
1574 'strategies': [AdaptivityPolynomialError(useMPI=True)],
1575 'num_procs': 1,
1576 'num_procs_sweeper': 4,
1577 }
1578 configurations[-1] = {
1579 'strategies': [
1580 ERKStrategy(useMPI=True),
1581 ARKStrategy(useMPI=True) if problem.__name__ in ['run_Schroedinger'] else ESDIRKStrategy(useMPI=True),
1582 ],
1583 'num_procs': 1,
1584 'custom_description': desc_RK,
1585 }
1587 configurations[2] = {
1588 'strategies': [AdaptivityStrategy(useMPI=True)],
1589 'custom_description': desc,
1590 'num_procs': 4,
1591 }
1592 elif mode == 'avoid_restarts':
1593 """
1594 Test how well avoiding restarts works.
1595 """
1596 from pySDC.projects.Resilience.strategies import (
1597 AdaptivityStrategy,
1598 AdaptivityAvoidRestartsStrategy,
1599 AdaptivityPolynomialStrategy,
1600 )
1602 desc = {'sweeper_params': {'QI': 'IE'}, 'step_params': {'maxiter': 3}}
1603 param_range = [1e-3, 1e-5]
1604 configurations[0] = {
1605 'strategies': [AdaptivityPolynomialStrategy(useMPI=True)],
1606 'plotting_params': {'ls': '--'},
1607 'custom_description': desc,
1608 'param_range': param_range,
1609 }
1610 configurations[1] = {
1611 'strategies': [AdaptivityAvoidRestartsStrategy(useMPI=True)],
1612 'plotting_params': {'ls': '-.'},
1613 'custom_description': desc,
1614 'param_range': param_range,
1615 }
1616 configurations[2] = {
1617 'strategies': [AdaptivityStrategy(useMPI=True)],
1618 'custom_description': desc,
1619 'param_range': param_range,
1620 }
1621 else:
1622 raise NotImplementedError(f'Don\'t know the mode "{mode}"!')
1624 return configurations
1627def get_fig(x=1, y=1, target='adaptivity', **kwargs): # pragma: no cover
1628 """
1629 Get a figure to plot in.
1631 Args:
1632 x (int): How many panels in horizontal direction you want
1633 y (int): How many panels in vertical direction you want
1634 target (str): Where the plot is supposed to end up
1636 Returns:
1637 matplotlib.pyplot.Figure
1638 """
1639 width = 1.0
1640 ratio = 1.0 if y == 2 else 0.5
1641 if target == 'adaptivity':
1642 journal = 'Springer_Numerical_Algorithms'
1643 elif target == 'thesis':
1644 journal = 'TUHH_thesis'
1645 elif target == 'talk':
1646 journal = 'JSC_beamer'
1647 else:
1648 raise NotImplementedError
1650 keyword_arguments = {
1651 'figsize': figsize_by_journal(journal, width, ratio),
1652 'layout': 'constrained',
1653 **kwargs,
1654 }
1655 return plt.subplots(y, x, **keyword_arguments)
1658def save_fig(
1659 fig, name, work_key, precision_key, legend=True, format='pdf', base_path='data', squares=True, ncols=None, **kwargs
1660): # pragma: no cover
1661 """
1662 Save a figure with a legend on the bottom.
1664 Args:
1665 fig (matplotlib.pyplot.Figure): Figure you want to save
1666 name (str): Name of the plot to put in the path
1667 work_key (str): The key in the recorded data you want on the x-axis
1668 precision_key (str): The key in the recorded data you want on the y-axis
1669 legend (bool): Put a legend or not
1670 format (str): Format to store the figure with
1671 base_path (str): Some path where all the files are stored
1672 squares (bool): Adjust aspect ratio to squares if true
1674 Returns:
1675 None
1676 """
1677 handles = []
1678 labels = []
1679 for ax in fig.get_axes():
1680 h, l = ax.get_legend_handles_labels()
1681 handles += [h[i] for i in range(len(h)) if l[i] not in labels]
1682 labels += [me for me in l if me not in labels]
1683 if squares:
1684 ax.set_box_aspect(1)
1685 # order = np.argsort([me[0] for me in labels])
1686 order = np.arange(len(labels))
1687 fig.legend(
1688 [handles[i] for i in order],
1689 [labels[i] for i in order],
1690 loc='outside lower center',
1691 ncols=ncols if ncols else 3 if len(handles) % 3 == 0 else 4,
1692 frameon=False,
1693 fancybox=True,
1694 handlelength=2.2,
1695 )
1697 path = f'{base_path}/wp-{name}-{work_key}-{precision_key}.{format}'
1698 fig.savefig(path, bbox_inches='tight', **kwargs)
1699 print(f'Stored figure \"{path}\"')
1702def all_problems(
1703 mode='compare_strategies', plotting=True, base_path='data', target='adaptivity', **kwargs
1704): # pragma: no cover
1705 """
1706 Make a plot comparing various strategies for all problems.
1708 Args:
1709 work_key (str): The key in the recorded data you want on the x-axis
1710 precision_key (str): The key in the recorded data you want on the y-axis
1712 Returns:
1713 None
1714 """
1716 if target == 'talk':
1717 fig, axs = get_fig(4, 1, target=target)
1718 else:
1719 fig, axs = get_fig(2, 2, target=target)
1721 shared_params = {
1722 'work_key': 'k_SDC',
1723 'precision_key': 'e_global',
1724 'num_procs': 1,
1725 'runs': 1,
1726 'comm_world': MPI.COMM_WORLD,
1727 'record': False,
1728 'plotting': plotting,
1729 **kwargs,
1730 }
1732 if target == 'adaptivity':
1733 problems = [run_vdp, run_quench, run_Schroedinger, run_AC]
1734 elif target in ['thesis', 'talk']:
1735 problems = [run_vdp, run_Lorenz, run_GS, run_RBC]
1736 else:
1737 raise NotImplementedError
1739 logger.log(26, f"Doing for all problems {mode}")
1740 for i in range(len(problems)):
1741 execute_configurations(
1742 **shared_params,
1743 problem=problems[i],
1744 ax=axs.flatten()[i],
1745 decorate=True,
1746 configurations=get_configs(mode, problems[i]),
1747 mode=mode,
1748 )
1750 if plotting and shared_params['comm_world'].rank == 0:
1751 ncols = {
1752 'parallel_efficiency': 2,
1753 'parallel_efficiency_dt': 2,
1754 'parallel_efficiency_dt_k': 2,
1755 'RK_comp': 2,
1756 }
1757 if target == 'talk':
1758 _ncols = 4
1759 else:
1760 _ncols = ncols.get(mode, None)
1762 if shared_params['work_key'] == 'param':
1763 for ax, prob in zip(fig.get_axes(), problems):
1764 add_param_order_lines(ax, prob)
1765 save_fig(
1766 fig=fig,
1767 name=mode,
1768 work_key=shared_params['work_key'],
1769 precision_key=shared_params['precision_key'],
1770 legend=True,
1771 base_path=base_path,
1772 ncols=_ncols,
1773 )
1776def add_param_order_lines(ax, problem):
1777 if problem.__name__ == 'run_vdp':
1778 yRfixed = 1e18
1779 yRdt = 1e-1
1780 yRdtk = 1e-4
1781 elif problem.__name__ == 'run_quench':
1782 yRfixed = 4e1
1783 yRdt = 1e4
1784 yRdtk = 1e4
1785 elif problem.__name__ == 'run_Schroedinger':
1786 yRfixed = 5
1787 yRdt = 1
1788 yRdtk = 1e-2
1789 elif problem.__name__ == 'run_AC':
1790 yRfixed = 1e8
1791 yRdt = 2e-2
1792 yRdtk = 1e-3
1793 elif problem.__name__ == 'run_Lorenz':
1794 yRfixed = 1e1
1795 yRdt = 2e-2
1796 yRdtk = 7e-4
1797 elif problem.__name__ == 'run_RBC':
1798 yRfixed = 1e-6
1799 yRdt = 4e-5
1800 yRdtk = 8e-6
1801 elif problem.__name__ == 'run_GS':
1802 yRfixed = 4e-3
1803 yRdt = 5e0
1804 yRdtk = 8e-1
1805 else:
1806 return None
1807 add_order_line(ax, 1, '--', yRdt, marker=None)
1808 add_order_line(ax, 5 / 4, ':', yRdtk, marker=None)
1809 add_order_line(ax, 5, '-.', yRfixed, marker=None)
1812def ODEs(mode='compare_strategies', plotting=True, base_path='data', **kwargs): # pragma: no cover
1813 """
1814 Make a plot comparing various strategies for the two ODEs.
1816 Args:
1817 work_key (str): The key in the recorded data you want on the x-axis
1818 precision_key (str): The key in the recorded data you want on the y-axis
1820 Returns:
1821 None
1822 """
1824 fig, axs = get_fig(x=2, y=1)
1826 shared_params = {
1827 'work_key': 'k_SDC',
1828 'precision_key': 'e_global',
1829 'num_procs': 1,
1830 'runs': 1,
1831 'comm_world': MPI.COMM_WORLD,
1832 'record': False,
1833 'plotting': plotting,
1834 **kwargs,
1835 }
1837 problems = [run_vdp, run_Lorenz]
1839 for i in range(len(problems)):
1840 execute_configurations(
1841 **shared_params,
1842 problem=problems[i],
1843 ax=axs.flatten()[i],
1844 decorate=i == 0,
1845 configurations=get_configs(mode, problems[i]),
1846 )
1848 if plotting and shared_params['comm_world'].rank == 0:
1849 save_fig(
1850 fig=fig,
1851 name=f'ODEs-{mode}',
1852 work_key=shared_params['work_key'],
1853 precision_key=shared_params['precision_key'],
1854 legend=True,
1855 base_path=base_path,
1856 )
1859def single_problem(mode, problem, plotting=True, base_path='data', target='thesis', **kwargs): # pragma: no cover
1860 """
1861 Make a plot for a single problem
1863 Args:
1864 mode (str): What you want to look at
1865 problem (function): A problem to run
1866 """
1867 if target == 'thesis':
1868 fig, ax = get_fig(1, 1, figsize=figsize_by_journal('TUHH_thesis', 0.7, 0.6))
1869 else:
1870 fig, ax = get_fig(1, 1, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1, 0.8))
1872 params = {
1873 'work_key': 'k_SDC',
1874 'precision_key': 'e_global',
1875 'num_procs': 1,
1876 'runs': 1,
1877 'comm_world': MPI.COMM_WORLD,
1878 'record': False,
1879 'plotting': plotting,
1880 **kwargs,
1881 }
1883 logger.log(26, f"Doing single problem {mode}")
1884 execute_configurations(
1885 **params, problem=problem, ax=ax, decorate=True, configurations=get_configs(mode, problem), mode=mode
1886 )
1888 if plotting:
1889 save_fig(
1890 fig=fig,
1891 name=f'{problem.__name__}-{mode}',
1892 work_key=params['work_key'],
1893 precision_key=params['precision_key'],
1894 legend=False,
1895 base_path=base_path,
1896 squares=target != 'thesis',
1897 )
1900def vdp_stiffness_plot(base_path='data', format='pdf', **kwargs): # pragma: no cover
1901 fig, axs = get_fig(3, 1, sharex=False, sharey=False)
1903 mus = [10, 100, 1000]
1905 for i in range(len(mus)):
1906 params = {
1907 'runs': 1,
1908 'problem': run_vdp,
1909 'record': False,
1910 'work_key': 't',
1911 'precision_key': 'e_global',
1912 'comm_world': MPI.COMM_WORLD,
1913 **kwargs,
1914 }
1915 params['num_procs'] = min(params['comm_world'].size, 5)
1916 params['plotting'] = params['comm_world'].rank == 0
1918 mode = f'vdp_stiffness-{mus[i]}'
1919 configurations = get_configs(mode=mode, problem=run_vdp)
1920 execute_configurations(**params, ax=axs.flatten()[i], decorate=True, configurations=configurations, mode=mode)
1921 axs.flatten()[i].set_title(rf'$\mu={ {mus[i]}} $')
1923 fig.suptitle('Van der Pol')
1924 if params['comm_world'].rank == 0:
1925 save_fig(
1926 fig=fig,
1927 name='vdp-stiffness',
1928 work_key=params['work_key'],
1929 precision_key=params['precision_key'],
1930 legend=False,
1931 base_path=base_path,
1932 format=format,
1933 )
1936def add_order_line(ax, order, ls, y_right=1.0, marker='.'):
1937 x_min = min([min(line.get_xdata()) for line in ax.get_lines()])
1938 x_max = max([max(line.get_xdata()) for line in ax.get_lines()])
1939 y_min = min([min(line.get_ydata()) for line in ax.get_lines()])
1940 y_max = max([max(line.get_ydata()) for line in ax.get_lines()])
1941 x = np.logspace(np.log10(x_min), np.log10(x_max), 100)
1942 y = y_right * (x / x_max) ** order
1943 mask = np.logical_and(y > y_min, y < y_max)
1944 ax.loglog(x[mask], y[mask], ls=ls, color='black', label=f'Order {order}', marker=marker, markevery=5)
1947def aggregate_parallel_efficiency_plot(): # pragma: no cover
1948 """
1949 Make a "weak scaling" plot for diagonal SDC
1950 """
1951 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
1953 fig, axs = plt.subplots(2, 2)
1955 _fig, _ax = plt.subplots(1, 1)
1956 num_procs = 1
1957 num_procs_sweeper = 2
1958 problem = run_quench
1960 num_procs_sweeper_list = [2, 3, 4]
1962 for problem, ax in zip([run_vdp, run_Lorenz, run_quench], axs.flatten()):
1963 speedup = []
1964 for num_procs_sweeper in num_procs_sweeper_list:
1965 s, e = plot_parallel_efficiency_diagonalSDC(
1966 ax=_ax,
1967 work_key='t',
1968 precision_key='e_global_rel',
1969 num_procs=num_procs,
1970 num_procs_sweeper=num_procs_sweeper,
1971 problem=problem,
1972 strategy=AdaptivityPolynomialError(),
1973 mode='diagonal_SDC',
1974 handle=f'{num_procs_sweeper} nodes',
1975 )
1976 speedup += [s]
1977 decorate_panel(ax, problem, work_key='nprocs', precision_key='')
1979 ax.plot(num_procs_sweeper_list, speedup, label='speedup')
1980 ax.plot(
1981 num_procs_sweeper_list,
1982 [speedup[i] / num_procs_sweeper_list[i] for i in range(len(speedup))],
1983 label='parallel efficiency',
1984 )
1986 fig.tight_layout()
1987 save_fig(fig, 'parallel_efficiency', 'nprocs', 'speedup')
1990if __name__ == "__main__":
1991 comm_world = MPI.COMM_WORLD
1993 import argparse
1995 parser = argparse.ArgumentParser()
1996 parser.add_argument('--mode', type=str, default='compare_strategies')
1997 parser.add_argument('--record', type=str, choices=['True', 'False'], default='True')
1998 parser.add_argument('--plotting', type=str, choices=['True', 'False'], default='True')
1999 parser.add_argument('--runs', type=int, default=5)
2000 parser.add_argument(
2001 '--problem', type=str, choices=['vdp', 'RBC', 'AC', 'quench', 'Lorenz', 'Schroedinger', 'GS'], default='vdp'
2002 )
2003 parser.add_argument('--work_key', type=str, default='t')
2004 parser.add_argument('--precision_key', type=str, default='e_global_rel')
2005 parser.add_argument('--logger_level', type=int, default='25')
2007 args = parser.parse_args()
2009 problems = {
2010 'Lorenz': run_Lorenz,
2011 'vdp': run_vdp,
2012 'Schroedinger': run_Schroedinger,
2013 'quench': run_quench,
2014 'AC': run_AC,
2015 'RBC': run_RBC,
2016 'GS': run_GS,
2017 }
2019 params = {
2020 **vars(args),
2021 'record': args.record == 'True',
2022 'plotting': args.plotting == 'True' and comm_world.rank == 0,
2023 'problem': problems[args.problem],
2024 }
2026 LOGGER_LEVEL = params.pop('logger_level')
2028 single_problem(**params)
2030 if comm_world.rank == 0:
2031 plt.show()