Coverage for pySDC/projects/Resilience/work_precision.py: 0%
486 statements
« prev ^ index » next coverage.py v7.8.0, created at 2025-04-18 08:18 +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", "ESDIRKStrategy", "ERKStrategy"]:
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 )
582 elif problem.__name__ == 'run_Lorenz':
583 if mode == 'parallel_efficiency_dt_k':
584 ax.set_xticks(
585 ticks=[
586 0.1,
587 0.2,
588 0.3,
589 0.4,
590 5e-1,
591 6e-1,
592 ],
593 labels=['', r'$2\times 10^{-1}$', '', r'$4\times 10^{-1}$', '', ''],
594 minor=True,
595 )
598def plot_parallel_efficiency_diagonalSDC(
599 ax, work_key, precision_key, num_procs_sweeper, num_procs=1, **kwargs
600): # pragma: no cover
601 serial_data = load(
602 num_procs=num_procs,
603 num_procs_sweeper=1,
604 **kwargs,
605 )
606 parallel_data = load(
607 num_procs=num_procs,
608 num_procs_sweeper=num_procs_sweeper,
609 **kwargs,
610 )
611 serial_work, serial_precision = extract_data(serial_data, work_key, precision_key)
612 parallel_work, parallel_precision = extract_data(parallel_data, work_key, precision_key)
613 # assert np.allclose(serial_precision, parallel_precision)
615 serial_work = np.asarray(serial_work)
616 parallel_work = np.asarray(parallel_work)
618 # ax.loglog(serial_work, serial_precision)
619 # ax.loglog(parallel_work, parallel_precision)
621 speedup = serial_work / parallel_work
622 parallel_efficiency = np.median(speedup) / num_procs_sweeper
623 ax.plot(serial_precision, speedup)
624 ax.set_xscale('log')
625 ax.set_ylabel('speedup')
627 if 't' in [work_key, precision_key]:
628 meta = parallel_data.get('meta', {})
630 if meta.get('hostname', None) in ['thomas-work']:
631 ax.text(0.1, 0.1, "Laptop timings!", transform=ax.transAxes)
632 if meta.get('runs', None) == 1:
633 ax.text(0.1, 0.2, "No sampling!", transform=ax.transAxes)
635 return np.median(speedup), parallel_efficiency
638def decorate_panel(ax, problem, work_key, precision_key, num_procs=1, title_only=False): # pragma: no cover
639 """
640 Decorate a plot
642 Args:
643 ax (matplotlib.pyplot.axes): Somewhere to plot
644 problem (function): A problem to run
645 work_key (str): The key in the recorded data you want on the x-axis
646 precision_key (str): The key in the recorded data you want on the y-axis
647 num_procs (int): Number of processes for the time communicator
648 title_only (bool): Put only the title on top, or do the whole shebang
650 Returns:
651 None
652 """
653 labels = {
654 'k_SDC': 'SDC iterations',
655 'k_SDC_no_restart': 'SDC iterations (restarts excluded)',
656 'k_Newton': 'Newton iterations',
657 'k_Newton_no_restart': 'Newton iterations (restarts excluded)',
658 'k_rhs': 'right hand side evaluations',
659 'k_factorizations': 'matrix factorizations',
660 't': 'wall clock time / s',
661 'e_global': 'global error',
662 'e_global_rel': 'relative global error',
663 'e_local_max': 'max. local error',
664 'restart': 'restarts',
665 'dt_max': r'$\Delta t_\mathrm{max}$',
666 'dt_min': r'$\Delta t_\mathrm{min}$',
667 'dt_sigma': r'$\sigma(\Delta t)$',
668 'dt_mean': r'$\bar{\Delta t}$',
669 'param': 'accuracy parameter',
670 'u0_increment_max': r'$\| \Delta u_0 \|_{\infty} $',
671 'u0_increment_mean': r'$\bar{\Delta u_0}$',
672 'u0_increment_max_no_restart': r'$\| \Delta u_0 \|_{\infty} $ (restarts excluded)',
673 'u0_increment_mean_no_restart': r'$\bar{\Delta u_0}$ (restarts excluded)',
674 'k_linear': 'Linear solver iterations',
675 'speedup': 'Speedup',
676 'nprocs': r'$N_\mathrm{procs}$',
677 '': '',
678 }
680 if not title_only:
681 ax.set_xlabel(labels.get(work_key, 'work'))
682 ax.set_ylabel(labels.get(precision_key, 'precision'))
683 # ax.legend(frameon=False)
685 titles = {
686 'run_vdp': 'Van der Pol',
687 'run_Lorenz': 'Lorenz attractor',
688 'run_Schroedinger': r'Schr\"odinger',
689 'run_quench': 'Quench',
690 'run_AC': 'Allen-Cahn',
691 'run_RBC': 'Rayleigh-Benard',
692 'run_GS': 'Gray-Scott',
693 }
694 ax.set_title(titles.get(problem.__name__, ''))
697def execute_configurations(
698 problem,
699 configurations,
700 work_key,
701 precision_key,
702 num_procs,
703 ax,
704 decorate,
705 record,
706 runs,
707 comm_world,
708 plotting,
709 Tend=None,
710 num_procs_sweeper=1,
711 mode='',
712):
713 """
714 Run for multiple configurations.
716 Args:
717 problem (function): A problem to run
718 configurations (dict): The configurations you want to run with
719 work_key (str): The key in the recorded data you want on the x-axis
720 precision_key (str): The key in the recorded data you want on the y-axis
721 num_procs (int): Number of processes for the time communicator
722 ax (matplotlib.pyplot.axes): Somewhere to plot
723 decorate (bool): Whether to decorate fully or only put the title
724 record (bool): Whether to only plot or also record the data first
725 runs (int): Number of runs you want to do
726 comm_world (mpi4py.MPI.Intracomm): Communicator that is available for the entire script
727 plotting (bool): Whether to plot something
728 num_procs_sweeper (int): Number of processes for the sweeper
729 mode (str): What you want to look at
731 Returns:
732 None
733 """
734 for _, config in configurations.items():
735 for strategy in config['strategies']:
736 shared_args = {
737 'problem': problem,
738 'strategy': strategy,
739 'handle': config.get('handle', ''),
740 'num_procs': config.get('num_procs', num_procs),
741 'num_procs_sweeper': config.get('num_procs_sweeper', num_procs_sweeper),
742 }
743 if record:
744 logger.debug('Recording work precision')
745 record_work_precision(
746 **shared_args,
747 custom_description=config.get('custom_description', {}),
748 runs=runs,
749 comm_world=comm_world,
750 problem_args=config.get('problem_args', {}),
751 param_range=config.get('param_range', None),
752 hooks=config.get('hooks', None),
753 Tend=config.get('Tend') if Tend is None else Tend,
754 mode=mode,
755 )
756 if plotting and comm_world.rank == 0:
757 logger.debug('Plotting')
758 plot_work_precision(
759 **shared_args,
760 work_key=work_key,
761 precision_key=precision_key,
762 ax=ax,
763 plotting_params=config.get('plotting_params', {}),
764 comm_world=comm_world,
765 mode=mode,
766 )
768 if comm_world.rank == 0:
769 decorate_panel(
770 ax=ax,
771 problem=problem,
772 work_key=work_key,
773 precision_key=precision_key,
774 num_procs=num_procs,
775 title_only=not decorate,
776 )
779def get_configs(mode, problem):
780 """
781 Get configurations for work-precision plots. These are dictionaries containing strategies and handles and so on.
783 Args:
784 mode (str): The of the configurations you want to retrieve
785 problem (function): A problem to run
787 Returns:
788 dict: Configurations
789 """
790 configurations = {}
791 if mode == 'regular':
792 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, BaseStrategy, IterateStrategy
794 handle = 'regular'
795 configurations[0] = {
796 'handle': handle,
797 'strategies': [AdaptivityStrategy(useMPI=True), BaseStrategy(useMPI=True), IterateStrategy(useMPI=True)],
798 }
799 elif mode == 'step_size_limiting':
800 from pySDC.implementations.convergence_controller_classes.step_size_limiter import StepSizeLimiter
801 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, ESDIRKStrategy
803 limits = [
804 25.0,
805 50.0,
806 ]
807 colors = ['teal', 'magenta']
808 markers = ['v', 'x']
809 markersize = 9
810 for i in range(len(limits)):
811 configurations[i] = {
812 'custom_description': {'convergence_controllers': {StepSizeLimiter: {'dt_max': limits[i]}}},
813 'handle': f'steplimiter{limits[i]:.0f}',
814 'strategies': [AdaptivityStrategy(useMPI=True)],
815 'plotting_params': {
816 'color': colors[i],
817 'marker': markers[i],
818 'label': rf'$\Delta t \leq { {limits[i]:.0f}} $',
819 # 'ls': '',
820 'markersize': markersize,
821 },
822 'num_procs': 1,
823 }
824 configurations[99] = {
825 'custom_description': {},
826 'handle': 'no limits',
827 'plotting_params': {
828 'label': 'no limiter',
829 # 'ls': '',
830 'markersize': markersize,
831 },
832 'strategies': [AdaptivityStrategy(useMPI=True)],
833 'num_procs': 1,
834 }
835 elif mode == 'dynamic_restarts':
836 """
837 Compare Block Gauss-Seidel SDC with restarting the first step in the block or the first step that exceeded the error threshold.
838 """
839 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, AdaptivityRestartFirstStep
841 desc = {}
842 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'}
843 desc['step_params'] = {'maxiter': 5}
845 ls = {
846 1: '-',
847 2: '--',
848 3: '-.',
849 4: ':',
850 5: ':',
851 }
853 configurations[-1] = {
854 'strategies': [AdaptivityStrategy(useMPI=True)],
855 'num_procs': 1,
856 }
858 for num_procs in [4, 2]:
859 plotting_params = {'ls': ls[num_procs], 'label': f'adaptivity {num_procs} procs'}
860 configurations[num_procs] = {
861 'strategies': [AdaptivityStrategy(useMPI=True), AdaptivityRestartFirstStep(useMPI=True)],
862 'custom_description': desc,
863 'num_procs': num_procs,
864 'plotting_params': plotting_params,
865 }
867 elif mode == 'compare_strategies':
868 """
869 Compare the different SDC strategies.
870 """
871 from pySDC.projects.Resilience.strategies import (
872 AdaptivityStrategy,
873 kAdaptivityStrategy,
874 AdaptivityPolynomialError,
875 BaseStrategy,
876 )
878 newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
880 configurations[1] = {
881 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
882 }
883 configurations[2] = {
884 'strategies': [kAdaptivityStrategy(useMPI=True)],
885 }
886 configurations[0] = {
887 'custom_description': {
888 'step_params': {'maxiter': 5},
889 'sweeper_params': {'num_nodes': 3, 'quad_type': 'RADAU-RIGHT'},
890 },
891 'strategies': [
892 BaseStrategy(useMPI=True),
893 AdaptivityStrategy(useMPI=True),
894 ],
895 }
897 elif mode == 'RK_comp':
898 """
899 Compare parallel adaptive SDC to Runge-Kutta
900 """
901 from pySDC.projects.Resilience.strategies import (
902 AdaptivityStrategy,
903 ERKStrategy,
904 ESDIRKStrategy,
905 ARKStrategy,
906 AdaptivityPolynomialError,
907 ARK3_CFL_Strategy,
908 )
910 if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_RBC', 'run_GS']:
911 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
912 else:
913 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
914 generic_implicit_MPI as parallel_sweeper,
915 )
917 newton_inexactness = False if problem.__name__ in ['run_vdp', 'run_RBC', 'run_GS'] else True
919 desc = {}
920 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': "EE"}
921 desc['step_params'] = {'maxiter': 5}
922 num_procs_dt = {
923 'run_RBC': 1,
924 }.get(problem.__name__, 4)
926 desc_poly = {}
927 desc_poly['sweeper_class'] = parallel_sweeper
928 num_procs_dt_k = 3
930 ls = {
931 1: '--',
932 2: '--',
933 3: '-',
934 4: '-',
935 5: '-',
936 12: ':',
937 }
938 RK_strategies = []
939 if problem.__name__ in ['run_Lorenz']:
940 RK_strategies.append(ERKStrategy(useMPI=True))
941 desc_poly['sweeper_params'] = {'QI': 'MIN-SR-S', 'QE': 'PIC'}
942 desc['sweeper_params']['QI'] = 'MIN-SR-S'
943 desc['sweeper_params']['QE'] = 'PIC'
944 if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS']:
945 RK_strategies.append(ARKStrategy(useMPI=True))
946 elif problem.__name__ == 'run_RBC':
947 RK_strategies.append(ARK3_CFL_Strategy(useMPI=True))
948 desc['sweeper_params']['num_nodes'] = 2
949 desc['sweeper_params']['QI'] = 'LU'
950 desc['sweeper_params']['QE'] = 'PIC'
951 desc['step_params']['maxiter'] = 3
953 desc_poly['sweeper_params'] = {'num_nodes': 2, 'QI': 'MIN-SR-S'}
954 num_procs_dt_k = 2
955 else:
956 RK_strategies.append(ESDIRKStrategy(useMPI=True))
958 configurations[-1] = {
959 'strategies': RK_strategies,
960 'num_procs': 1,
961 }
962 configurations[3] = {
963 'custom_description': desc_poly,
964 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
965 'num_procs': 1,
966 'num_procs_sweeper': num_procs_dt_k,
967 'plotting_params': {
968 'label': rf'$\Delta t$-$k$-adaptivity $N$=1x{num_procs_dt_k}',
969 'ls': ls[num_procs_dt_k],
970 },
971 }
972 if problem.__name__ in ['run_Lorenz']:
973 configurations[2] = {
974 'strategies': [AdaptivityStrategy(useMPI=True)],
975 'custom_description': {**desc, 'sweeper_class': parallel_sweeper},
976 'num_procs': num_procs_dt,
977 'num_procs_sweeper': num_procs_dt_k,
978 'plotting_params': {
979 'label': rf'$\Delta t$-adaptivity $N$={num_procs_dt}x3',
980 'ls': ls[num_procs_dt * num_procs_dt_k],
981 },
982 }
983 else:
984 configurations[2] = {
985 'strategies': [AdaptivityStrategy(useMPI=True)],
986 'custom_description': desc,
987 'num_procs': num_procs_dt,
988 'plotting_params': {'label': rf'$\Delta t$-adaptivity $N$={num_procs_dt}x1', 'ls': ls[num_procs_dt]},
989 }
991 elif mode == 'RK_comp_high_order_RBC':
992 """
993 Compare parallel adaptive SDC to Runge-Kutta at order five for RBC problem
994 """
995 from pySDC.projects.Resilience.strategies import (
996 AdaptivityStrategy,
997 ERKStrategy,
998 ESDIRKStrategy,
999 ARKStrategy,
1000 AdaptivityPolynomialError,
1001 ARK3_CFL_Strategy,
1002 )
1004 assert problem.__name__ == 'run_RBC'
1006 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1008 newton_inexactness = False
1010 desc = {}
1011 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': "EE"}
1012 desc['step_params'] = {'maxiter': 5}
1013 num_procs_dt = 1
1015 desc_poly = {}
1016 desc_poly['sweeper_class'] = parallel_sweeper
1017 num_procs_dt_k = 3
1019 ls = {
1020 1: '--',
1021 2: '--',
1022 3: '-',
1023 4: '-',
1024 5: '-',
1025 12: ':',
1026 }
1027 RK_strategies = [ARK3_CFL_Strategy(useMPI=True)]
1028 desc['sweeper_params']['num_nodes'] = 3
1029 desc['sweeper_params']['QI'] = 'LU'
1030 desc['sweeper_params']['QE'] = 'PIC'
1031 desc['step_params']['maxiter'] = 5
1033 desc_poly['sweeper_params'] = {'num_nodes': 3, 'QI': 'MIN-SR-S'}
1034 num_procs_dt_k = 3
1036 configurations[-1] = {
1037 'strategies': RK_strategies,
1038 'num_procs': 1,
1039 }
1040 configurations[3] = {
1041 'custom_description': desc_poly,
1042 'strategies': [AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness)],
1043 'num_procs': 1,
1044 'num_procs_sweeper': num_procs_dt_k,
1045 'plotting_params': {
1046 'label': rf'$\Delta t$-$k$-adaptivity $N$=1x{num_procs_dt_k}',
1047 'ls': ls[num_procs_dt_k],
1048 },
1049 }
1050 configurations[2] = {
1051 'strategies': [AdaptivityStrategy(useMPI=True)],
1052 'custom_description': desc,
1053 'num_procs': num_procs_dt,
1054 'plotting_params': {'label': rf'$\Delta t$-adaptivity $N$={num_procs_dt}x1', 'ls': ls[num_procs_dt]},
1055 }
1057 elif mode == 'parallel_efficiency':
1058 """
1059 Compare parallel runs of the step size adaptive SDC
1060 """
1061 from pySDC.projects.Resilience.strategies import AdaptivityStrategy, AdaptivityPolynomialError
1063 if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS', 'run_RBC']:
1064 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1065 else:
1066 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1067 generic_implicit_MPI as parallel_sweeper,
1068 )
1070 desc = {}
1071 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': 'EE'}
1072 desc['step_params'] = {'maxiter': 5}
1074 if problem.__name__ in ['run_RBC']:
1075 desc['sweeper_params']['QE'] = 'PIC'
1076 desc['sweeper_params']['QI'] = 'LU'
1078 ls = {
1079 1: '-',
1080 2: '--',
1081 3: '-.',
1082 4: '--',
1083 5: ':',
1084 12: ':',
1085 }
1087 newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
1089 for num_procs in [4, 1]:
1090 plotting_params = (
1091 {'ls': ls[num_procs], 'label': fr'$\Delta t$-adaptivity $N$={num_procs}x1'} if num_procs > 1 else {}
1092 )
1093 configurations[num_procs] = {
1094 'strategies': [AdaptivityStrategy(useMPI=True)],
1095 'custom_description': desc.copy(),
1096 'num_procs': num_procs,
1097 'plotting_params': plotting_params.copy(),
1098 }
1099 configurations[num_procs * 100 + 79] = {
1100 'custom_description': {'sweeper_class': parallel_sweeper},
1101 'strategies': [
1102 AdaptivityPolynomialError(
1103 useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True
1104 )
1105 ],
1106 'num_procs_sweeper': 3,
1107 'num_procs': num_procs,
1108 'plotting_params': {
1109 'ls': ls.get(num_procs * 3, '-'),
1110 'label': rf'$\Delta t$-$k$-adaptivity $N$={num_procs}x3',
1111 },
1112 }
1114 configurations[200 + 79] = {
1115 'strategies': [
1116 AdaptivityPolynomialError(useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True)
1117 ],
1118 'num_procs': 1,
1119 }
1120 elif mode == 'parallel_efficiency_dt':
1121 """
1122 Compare parallel runs of the step size adaptive SDC
1123 """
1124 from pySDC.projects.Resilience.strategies import AdaptivityStrategy
1126 if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS', 'run_RBC']:
1127 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1128 else:
1129 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1130 generic_implicit_MPI as parallel_sweeper,
1131 )
1133 desc = {}
1134 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE', 'QE': 'EE'}
1135 desc['step_params'] = {'maxiter': 5}
1137 if problem.__name__ in ['run_RBC']:
1138 desc['sweeper_params']['QE'] = 'PIC'
1139 desc['sweeper_params']['QI'] = 'LU'
1141 desc_serial = {
1142 'step_params': {'maxiter': 5},
1143 'sweeper_params': {'num_nodes': 3, 'quad_type': 'RADAU-RIGHT'},
1144 }
1146 ls = {
1147 1: '-',
1148 2: '--',
1149 3: '-.',
1150 4: '--',
1151 5: ':',
1152 12: ':',
1153 }
1155 newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
1157 for num_procs in [4, 1]:
1158 configurations[num_procs] = {
1159 'strategies': [AdaptivityStrategy(useMPI=True)],
1160 'custom_description': desc.copy() if num_procs > 1 else desc_serial,
1161 'num_procs': num_procs,
1162 'plotting_params': {
1163 'ls': ls.get(num_procs, '-'),
1164 'label': rf'$\Delta t$-adaptivity $N$={num_procs}x1',
1165 },
1166 }
1167 configurations[num_procs * 200 + 79] = {
1168 'custom_description': {
1169 'sweeper_class': parallel_sweeper,
1170 'sweeper_params': {'QI': 'MIN-SR-S', 'QE': 'PIC'},
1171 'step_params': {'maxiter': 5},
1172 },
1173 'strategies': [AdaptivityStrategy(useMPI=True)],
1174 'num_procs_sweeper': 3,
1175 'num_procs': num_procs,
1176 'plotting_params': {
1177 'ls': ls.get(num_procs * 3, '-'),
1178 'label': rf'$\Delta t$-adaptivity $N$={num_procs}x3',
1179 },
1180 }
1181 elif mode == 'parallel_efficiency_dt_k':
1182 """
1183 Compare parallel runs of the step size adaptive SDC
1184 """
1185 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
1187 if problem.__name__ in ['run_Schroedinger', 'run_AC', 'run_GS', 'run_RBC']:
1188 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1189 else:
1190 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1191 generic_implicit_MPI as parallel_sweeper,
1192 )
1194 ls = {
1195 1: '-',
1196 2: '--',
1197 3: '-.',
1198 4: '--',
1199 5: ':',
1200 12: ':',
1201 }
1203 QI = {
1204 (1, 3, 'run_Lorenz'): 'MIN-SR-NS',
1205 (1, 1, 'run_Lorenz'): 'MIN-SR-NS',
1206 (4, 1, 'run_Lorenz'): 'MIN-SR-NS',
1207 }
1209 newton_inexactness = False if problem.__name__ in ['run_vdp'] else True
1211 for num_procs in [4, 1]:
1212 configurations[num_procs * 100 + 79] = {
1213 'custom_description': {
1214 'sweeper_class': parallel_sweeper,
1215 'sweeper_params': {'QI': QI.get((num_procs, 3, problem.__name__), 'MIN-SR-S'), 'QE': 'PIC'},
1216 },
1217 'strategies': [
1218 AdaptivityPolynomialError(
1219 useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True
1220 )
1221 ],
1222 'num_procs_sweeper': 3,
1223 'num_procs': num_procs,
1224 'plotting_params': {
1225 'ls': ls.get(num_procs * 3, '-'),
1226 'label': rf'$\Delta t$-$k$-adaptivity $N$={num_procs}x3',
1227 },
1228 }
1229 configurations[num_procs * 200 + 79] = {
1230 'custom_description': {
1231 'sweeper_params': {'QI': QI.get((num_procs, 1, problem.__name__), 'MIN-SR-S'), 'QE': 'PIC'},
1232 },
1233 'strategies': [
1234 AdaptivityPolynomialError(
1235 useMPI=True, newton_inexactness=newton_inexactness, linear_inexactness=True
1236 )
1237 ],
1238 'num_procs_sweeper': 1,
1239 'num_procs': num_procs,
1240 'plotting_params': {
1241 'ls': ls.get(num_procs, '-'),
1242 'label': rf'$\Delta t$-$k$-adaptivity $N$={num_procs}x1',
1243 },
1244 }
1245 elif mode == 'interpolate_between_restarts':
1246 """
1247 Compare adaptivity with interpolation between restarts and without
1248 """
1249 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
1251 i = 0
1252 for interpolate_between_restarts, handle, ls in zip(
1253 [True, False], ['Interpolation between restarts', 'regular'], ['--', '-']
1254 ):
1255 configurations[i] = {
1256 'strategies': [
1257 AdaptivityPolynomialError(interpolate_between_restarts=interpolate_between_restarts, useMPI=True)
1258 ],
1259 'plotting_params': {'ls': ls},
1260 'handle': handle,
1261 }
1262 i += 1
1263 elif mode == 'diagonal_SDC':
1264 """
1265 Run diagonal SDC with different number of nodes and ranks. You can use this to compute a speedup, but it's not strong scaling!
1266 """
1267 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
1269 if problem.__name__ in ['run_Schroedinger']:
1270 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1271 else:
1272 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1273 generic_implicit_MPI as parallel_sweeper,
1274 )
1276 for parallel in [False, True]:
1277 desc = {'sweeper_class': parallel_sweeper} if parallel else {}
1278 for num_nodes, ls in zip([3, 4, 2], ['-', '--', ':', '-.']):
1279 configurations[num_nodes + (99 if parallel else 0)] = {
1280 'custom_description': {**desc, 'sweeper_params': {'num_nodes': num_nodes}},
1281 'strategies': [
1282 AdaptivityPolynomialError(useMPI=True, newton_inexactness=True, linear_inexactness=True)
1283 ],
1284 'num_procs_sweeper': num_nodes if parallel else 1,
1285 'num_procs': 1,
1286 'handle': f'{num_nodes} nodes',
1287 'plotting_params': {
1288 'ls': ls,
1289 'label': f'{num_nodes} procs',
1290 # **{'color': 'grey' if parallel else None},
1291 },
1292 }
1294 elif mode[:13] == 'vdp_stiffness':
1295 """
1296 Run van der Pol with different parameter for the nonlinear term, which controls the stiffness.
1297 """
1298 from pySDC.projects.Resilience.strategies import (
1299 AdaptivityStrategy,
1300 ERKStrategy,
1301 ESDIRKStrategy,
1302 AdaptivityPolynomialError,
1303 )
1304 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1305 generic_implicit_MPI as parallel_sweeper,
1306 )
1308 Tends = {
1309 1000: 2000,
1310 100: 200,
1311 10: 20,
1312 0: 2,
1313 }
1314 mu = float(mode[14:])
1315 Tend = Tends[mu]
1317 problem_desc = {'problem_params': {'mu': mu}}
1319 desc = {}
1320 desc['sweeper_params'] = {'num_nodes': 3, 'QI': 'IE'}
1321 desc['step_params'] = {'maxiter': 5}
1322 desc['problem_params'] = problem_desc['problem_params']
1324 ls = {
1325 1: '-',
1326 2: '--',
1327 3: '-.',
1328 4: ':',
1329 5: ':',
1330 'MIN-SR-S': '-',
1331 'MIN-SR-NS': '--',
1332 'MIN-SR-FLEX': '-.',
1333 }
1335 if mu < 100:
1336 configurations[2] = {
1337 'strategies': [ERKStrategy(useMPI=True)],
1338 'num_procs': 1,
1339 'handle': mode,
1340 'plotting_params': {'label': 'CP5(4)'},
1341 'custom_description': problem_desc,
1342 'Tend': Tend,
1343 }
1344 configurations[1] = {
1345 'strategies': [AdaptivityStrategy(useMPI=True)],
1346 'custom_description': desc,
1347 'num_procs': 4,
1348 'plotting_params': {'ls': ls[1], 'label': 'SDC $N$=4x1'},
1349 'handle': mode,
1350 'Tend': Tend,
1351 }
1352 configurations[4] = {
1353 'strategies': [ESDIRKStrategy(useMPI=True)],
1354 'num_procs': 1,
1355 'handle': mode,
1356 'plotting_params': {'label': 'ESDIRK5(3)'},
1357 'custom_description': problem_desc,
1358 'Tend': Tend,
1359 }
1360 for QI, i in zip(
1361 [
1362 'MIN-SR-S',
1363 # 'MIN-SR-FLEX',
1364 ],
1365 [9991, 12123127391, 1231723109247102731092],
1366 ):
1367 configurations[i] = {
1368 'custom_description': {
1369 'sweeper_params': {'num_nodes': 3, 'QI': QI},
1370 'problem_params': desc["problem_params"],
1371 'sweeper_class': parallel_sweeper,
1372 },
1373 'strategies': [
1374 AdaptivityPolynomialError(
1375 useMPI=True, newton_inexactness=False, linear_inexactness=False, max_slope=4
1376 )
1377 ],
1378 'num_procs_sweeper': 3,
1379 'num_procs': 1,
1380 'plotting_params': {
1381 'ls': ls.get(QI, '-'),
1382 'label': rf'$\Delta t$-$k$-adaptivity $N$={1}x3',
1383 },
1384 'handle': f'{mode}-{QI}',
1385 'Tend': Tend,
1386 }
1388 elif mode == 'inexactness':
1389 """
1390 Compare inexact SDC to exact SDC
1391 """
1392 from pySDC.projects.Resilience.strategies import (
1393 AdaptivityPolynomialError,
1394 )
1396 if problem.__name__ in ['run_Schroedinger']:
1397 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1398 else:
1399 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1400 generic_implicit_MPI as parallel_sweeper,
1401 )
1403 strategies = [
1404 AdaptivityPolynomialError,
1405 ]
1407 inexactness = {
1408 'newton_inexactness': True,
1409 'linear_inexactness': True,
1410 }
1411 no_inexactness = {
1412 'newton_inexactness': False,
1413 'linear_inexactness': False,
1414 'SDC_maxiter': 99,
1415 'use_restol_rel': False,
1416 }
1418 configurations[1] = {
1419 'custom_description': {'sweeper_class': parallel_sweeper},
1420 'strategies': [me(useMPI=True, **no_inexactness) for me in strategies],
1421 'num_procs_sweeper': 3,
1422 'handle': 'exact',
1423 'plotting_params': {'ls': '--'},
1424 }
1425 configurations[0] = {
1426 'custom_description': {'sweeper_class': parallel_sweeper},
1427 'strategies': [me(useMPI=True, **inexactness) for me in strategies],
1428 'handle': 'inexact',
1429 'num_procs_sweeper': 3,
1430 }
1431 elif mode == 'compare_adaptivity':
1432 """
1433 Compare various modes of adaptivity
1434 """
1435 # TODO: configurations not final!
1436 from pySDC.projects.Resilience.strategies import (
1437 # AdaptivityCollocationTypeStrategy,
1438 # AdaptivityCollocationRefinementStrategy,
1439 AdaptivityStrategy,
1440 # AdaptivityExtrapolationWithinQStrategy,
1441 ESDIRKStrategy,
1442 ARKStrategy,
1443 AdaptivityPolynomialError,
1444 )
1446 if problem.__name__ in ['run_Schroedinger']:
1447 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1448 else:
1449 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1450 generic_implicit_MPI as parallel_sweeper,
1451 )
1453 inexactness_params = {
1454 # 'double_adaptivity': True,
1455 'newton_inexactness': True,
1456 'linear_inexactness': True,
1457 }
1459 strategies = [
1460 AdaptivityPolynomialError,
1461 # AdaptivityCollocationTypeStrategy,
1462 # AdaptivityExtrapolationWithinQStrategy,
1463 ]
1465 # restol = None
1466 # for strategy in strategies:
1467 # strategy.restol = restol
1469 configurations[1] = {
1470 'custom_description': {'sweeper_class': parallel_sweeper},
1471 'strategies': [me(useMPI=True, **inexactness_params) for me in strategies],
1472 'handle': 'parallel',
1473 'num_procs_sweeper': 3,
1474 'plotting_params': {'ls': '-', 'label': '3 procs'},
1475 }
1476 configurations[2] = {
1477 'strategies': [me(useMPI=True, **inexactness_params) for me in strategies],
1478 'plotting_params': {'ls': '--'},
1479 }
1480 configurations[4] = {
1481 'custom_description': {'step_params': {'maxiter': 5}},
1482 'strategies': [AdaptivityStrategy(useMPI=True)],
1483 }
1485 desc_RK = {}
1486 configurations[-1] = {
1487 'strategies': [
1488 ARKStrategy(useMPI=True) if problem.__name__ == 'run_Schroedinger' else ESDIRKStrategy(useMPI=True),
1489 ],
1490 'num_procs': 1,
1491 'custom_description': desc_RK,
1492 }
1494 elif mode == 'preconditioners':
1495 """
1496 Compare different preconditioners
1497 """
1498 from pySDC.projects.Resilience.strategies import (
1499 AdaptivityStrategy,
1500 IterateStrategy,
1501 BaseStrategy,
1502 ESDIRKStrategy,
1503 ERKStrategy,
1504 AdaptivityPolynomialError,
1505 )
1507 inexacness = True
1508 strategies = [
1509 AdaptivityPolynomialError(
1510 useMPI=True, SDC_maxiter=29, newton_inexactness=inexacness, linear_inexactness=inexacness
1511 ),
1512 BaseStrategy(useMPI=True),
1513 ]
1515 desc = {}
1516 desc['sweeper_params'] = {
1517 'num_nodes': 3,
1518 }
1519 # desc['step_params'] = {'maxiter': 5}
1521 precons = ['IE', 'LU']
1522 ls = ['-.', '--', '-', ':']
1523 for i in range(len(precons) + 1):
1524 if i < len(precons):
1525 desc['sweeper_params']['QI'] = precons[i]
1526 handle = precons[i]
1527 else:
1528 handle = None
1529 configurations[i] = {
1530 'strategies': strategies,
1531 'custom_description': copy.deepcopy(desc),
1532 'handle': handle,
1533 'plotting_params': {'ls': ls[i]},
1534 }
1535 elif mode == 'RK_comp_high_order':
1536 """
1537 Compare higher order SDC than we can get with RKM to RKM
1538 """
1539 from pySDC.projects.Resilience.strategies import (
1540 AdaptivityStrategy,
1541 ERKStrategy,
1542 ESDIRKStrategy,
1543 ARKStrategy,
1544 AdaptivityPolynomialError,
1545 )
1547 if problem.__name__ in ['run_Schroedinger']:
1548 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1549 else:
1550 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1551 generic_implicit_MPI as parallel_sweeper,
1552 )
1554 desc = {}
1555 desc['sweeper_params'] = {'num_nodes': 4, 'QI': 'IE', 'QE': "EE"}
1556 desc['step_params'] = {'maxiter': 7}
1558 desc_poly = {}
1559 desc_poly['sweeper_class'] = parallel_sweeper
1560 desc_poly['sweeper_params'] = {'num_nodes': 4}
1562 ls = {
1563 1: '-',
1564 2: '--',
1565 3: '-.',
1566 4: ':',
1567 5: ':',
1568 }
1570 desc_RK = {}
1571 if problem.__name__ in ['run_Schroedinger']:
1572 desc_RK['problem_params'] = {'imex': True}
1574 configurations[3] = {
1575 'custom_description': desc_poly,
1576 'strategies': [AdaptivityPolynomialError(useMPI=True)],
1577 'num_procs': 1,
1578 'num_procs_sweeper': 4,
1579 }
1580 configurations[-1] = {
1581 'strategies': [
1582 ERKStrategy(useMPI=True),
1583 ARKStrategy(useMPI=True) if problem.__name__ in ['run_Schroedinger'] else ESDIRKStrategy(useMPI=True),
1584 ],
1585 'num_procs': 1,
1586 'custom_description': desc_RK,
1587 }
1589 configurations[2] = {
1590 'strategies': [AdaptivityStrategy(useMPI=True)],
1591 'custom_description': desc,
1592 'num_procs': 4,
1593 }
1594 elif mode == 'avoid_restarts':
1595 """
1596 Test how well avoiding restarts works.
1597 """
1598 from pySDC.projects.Resilience.strategies import (
1599 AdaptivityStrategy,
1600 AdaptivityAvoidRestartsStrategy,
1601 AdaptivityPolynomialStrategy,
1602 )
1604 desc = {'sweeper_params': {'QI': 'IE'}, 'step_params': {'maxiter': 3}}
1605 param_range = [1e-3, 1e-5]
1606 configurations[0] = {
1607 'strategies': [AdaptivityPolynomialStrategy(useMPI=True)],
1608 'plotting_params': {'ls': '--'},
1609 'custom_description': desc,
1610 'param_range': param_range,
1611 }
1612 configurations[1] = {
1613 'strategies': [AdaptivityAvoidRestartsStrategy(useMPI=True)],
1614 'plotting_params': {'ls': '-.'},
1615 'custom_description': desc,
1616 'param_range': param_range,
1617 }
1618 configurations[2] = {
1619 'strategies': [AdaptivityStrategy(useMPI=True)],
1620 'custom_description': desc,
1621 'param_range': param_range,
1622 }
1623 else:
1624 raise NotImplementedError(f'Don\'t know the mode "{mode}"!')
1626 return configurations
1629def get_fig(x=1, y=1, target='adaptivity', **kwargs): # pragma: no cover
1630 """
1631 Get a figure to plot in.
1633 Args:
1634 x (int): How many panels in horizontal direction you want
1635 y (int): How many panels in vertical direction you want
1636 target (str): Where the plot is supposed to end up
1638 Returns:
1639 matplotlib.pyplot.Figure
1640 """
1641 width = 1.0
1642 ratio = 1.0 if y == 2 else 0.5
1643 if target == 'adaptivity':
1644 journal = 'Springer_Numerical_Algorithms'
1645 elif target == 'thesis':
1646 journal = 'TUHH_thesis'
1647 elif target == 'talk':
1648 journal = 'JSC_beamer'
1649 else:
1650 raise NotImplementedError
1652 keyword_arguments = {
1653 'figsize': figsize_by_journal(journal, width, ratio),
1654 'layout': 'constrained',
1655 **kwargs,
1656 }
1657 return plt.subplots(y, x, **keyword_arguments)
1660def save_fig(
1661 fig, name, work_key, precision_key, legend=True, format='pdf', base_path='data', squares=True, ncols=None, **kwargs
1662): # pragma: no cover
1663 """
1664 Save a figure with a legend on the bottom.
1666 Args:
1667 fig (matplotlib.pyplot.Figure): Figure you want to save
1668 name (str): Name of the plot to put in the path
1669 work_key (str): The key in the recorded data you want on the x-axis
1670 precision_key (str): The key in the recorded data you want on the y-axis
1671 legend (bool): Put a legend or not
1672 format (str): Format to store the figure with
1673 base_path (str): Some path where all the files are stored
1674 squares (bool): Adjust aspect ratio to squares if true
1676 Returns:
1677 None
1678 """
1679 handles = []
1680 labels = []
1681 for ax in fig.get_axes():
1682 h, l = ax.get_legend_handles_labels()
1683 handles += [h[i] for i in range(len(h)) if l[i] not in labels]
1684 labels += [me for me in l if me not in labels]
1685 if squares:
1686 ax.set_box_aspect(1)
1687 # order = np.argsort([me[0] for me in labels])
1688 order = np.arange(len(labels))
1689 fig.legend(
1690 [handles[i] for i in order],
1691 [labels[i] for i in order],
1692 loc='outside lower center',
1693 ncols=ncols if ncols else 3 if len(handles) % 3 == 0 else 4,
1694 frameon=False,
1695 fancybox=True,
1696 handlelength=2.2,
1697 )
1699 path = f'{base_path}/wp-{name}-{work_key}-{precision_key}.{format}'
1700 fig.savefig(path, bbox_inches='tight', **kwargs)
1701 print(f'Stored figure \"{path}\"')
1704def all_problems(
1705 mode='compare_strategies', plotting=True, base_path='data', target='adaptivity', **kwargs
1706): # pragma: no cover
1707 """
1708 Make a plot comparing various strategies for all problems.
1710 Args:
1711 work_key (str): The key in the recorded data you want on the x-axis
1712 precision_key (str): The key in the recorded data you want on the y-axis
1714 Returns:
1715 None
1716 """
1718 if target == 'talk':
1719 fig, axs = get_fig(4, 1, target=target)
1720 else:
1721 fig, axs = get_fig(2, 2, target=target)
1723 shared_params = {
1724 'work_key': 'k_SDC',
1725 'precision_key': 'e_global',
1726 'num_procs': 1,
1727 'runs': 1,
1728 'comm_world': MPI.COMM_WORLD,
1729 'record': False,
1730 'plotting': plotting,
1731 **kwargs,
1732 }
1734 if target == 'adaptivity':
1735 problems = [run_vdp, run_quench, run_Schroedinger, run_AC]
1736 elif target in ['thesis', 'talk']:
1737 problems = [run_vdp, run_Lorenz, run_GS, run_RBC]
1738 else:
1739 raise NotImplementedError
1741 logger.log(26, f"Doing for all problems {mode}")
1742 for i in range(len(problems)):
1743 execute_configurations(
1744 **shared_params,
1745 problem=problems[i],
1746 ax=axs.flatten()[i],
1747 decorate=True,
1748 configurations=get_configs(mode, problems[i]),
1749 mode=mode,
1750 )
1752 if plotting and shared_params['comm_world'].rank == 0:
1753 ncols = {
1754 'parallel_efficiency': 2,
1755 'parallel_efficiency_dt': 2,
1756 'parallel_efficiency_dt_k': 2,
1757 'RK_comp': 2,
1758 }
1759 if target == 'talk':
1760 _ncols = 4
1761 else:
1762 _ncols = ncols.get(mode, None)
1764 if shared_params['work_key'] == 'param':
1765 for ax, prob in zip(fig.get_axes(), problems):
1766 add_param_order_lines(ax, prob)
1767 save_fig(
1768 fig=fig,
1769 name=mode,
1770 work_key=shared_params['work_key'],
1771 precision_key=shared_params['precision_key'],
1772 legend=True,
1773 base_path=base_path,
1774 ncols=_ncols,
1775 )
1778def add_param_order_lines(ax, problem):
1779 if problem.__name__ == 'run_vdp':
1780 yRfixed = 1e18
1781 yRdt = 1e-1
1782 yRdtk = 1e-4
1783 elif problem.__name__ == 'run_quench':
1784 yRfixed = 4e1
1785 yRdt = 1e4
1786 yRdtk = 1e4
1787 elif problem.__name__ == 'run_Schroedinger':
1788 yRfixed = 5
1789 yRdt = 1
1790 yRdtk = 1e-2
1791 elif problem.__name__ == 'run_AC':
1792 yRfixed = 1e8
1793 yRdt = 2e-2
1794 yRdtk = 1e-3
1795 elif problem.__name__ == 'run_Lorenz':
1796 yRfixed = 1e1
1797 yRdt = 2e-2
1798 yRdtk = 7e-4
1799 elif problem.__name__ == 'run_RBC':
1800 yRfixed = 1e-6
1801 yRdt = 4e-5
1802 yRdtk = 8e-6
1803 elif problem.__name__ == 'run_GS':
1804 yRfixed = 4e-3
1805 yRdt = 5e0
1806 yRdtk = 8e-1
1807 else:
1808 return None
1809 add_order_line(ax, 1, '--', yRdt, marker=None)
1810 add_order_line(ax, 5 / 4, ':', yRdtk, marker=None)
1811 add_order_line(ax, 5, '-.', yRfixed, marker=None)
1814def ODEs(mode='compare_strategies', plotting=True, base_path='data', **kwargs): # pragma: no cover
1815 """
1816 Make a plot comparing various strategies for the two ODEs.
1818 Args:
1819 work_key (str): The key in the recorded data you want on the x-axis
1820 precision_key (str): The key in the recorded data you want on the y-axis
1822 Returns:
1823 None
1824 """
1826 fig, axs = get_fig(x=2, y=1)
1828 shared_params = {
1829 'work_key': 'k_SDC',
1830 'precision_key': 'e_global',
1831 'num_procs': 1,
1832 'runs': 1,
1833 'comm_world': MPI.COMM_WORLD,
1834 'record': False,
1835 'plotting': plotting,
1836 **kwargs,
1837 }
1839 problems = [run_vdp, run_Lorenz]
1841 for i in range(len(problems)):
1842 execute_configurations(
1843 **shared_params,
1844 problem=problems[i],
1845 ax=axs.flatten()[i],
1846 decorate=i == 0,
1847 configurations=get_configs(mode, problems[i]),
1848 )
1850 if plotting and shared_params['comm_world'].rank == 0:
1851 save_fig(
1852 fig=fig,
1853 name=f'ODEs-{mode}',
1854 work_key=shared_params['work_key'],
1855 precision_key=shared_params['precision_key'],
1856 legend=True,
1857 base_path=base_path,
1858 )
1861def single_problem(mode, problem, plotting=True, base_path='data', target='thesis', **kwargs): # pragma: no cover
1862 """
1863 Make a plot for a single problem
1865 Args:
1866 mode (str): What you want to look at
1867 problem (function): A problem to run
1868 """
1869 if target == 'thesis':
1870 fig, ax = get_fig(1, 1, figsize=figsize_by_journal('TUHH_thesis', 0.7, 0.6))
1871 else:
1872 fig, ax = get_fig(1, 1, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1, 0.8))
1874 params = {
1875 'work_key': 'k_SDC',
1876 'precision_key': 'e_global',
1877 'num_procs': 1,
1878 'runs': 1,
1879 'comm_world': MPI.COMM_WORLD,
1880 'record': False,
1881 'plotting': plotting,
1882 **kwargs,
1883 }
1885 logger.log(26, f"Doing single problem {mode}")
1886 execute_configurations(
1887 **params, problem=problem, ax=ax, decorate=True, configurations=get_configs(mode, problem), mode=mode
1888 )
1890 if plotting:
1891 save_fig(
1892 fig=fig,
1893 name=f'{problem.__name__}-{mode}',
1894 work_key=params['work_key'],
1895 precision_key=params['precision_key'],
1896 legend=False,
1897 base_path=base_path,
1898 squares=target != 'thesis',
1899 )
1902def vdp_stiffness_plot(base_path='data', format='pdf', **kwargs): # pragma: no cover
1903 fig, axs = get_fig(3, 1, sharex=False, sharey=False)
1905 mus = [10, 100, 1000]
1907 for i in range(len(mus)):
1908 params = {
1909 'runs': 1,
1910 'problem': run_vdp,
1911 'record': False,
1912 'work_key': 't',
1913 'precision_key': 'e_global',
1914 'comm_world': MPI.COMM_WORLD,
1915 **kwargs,
1916 }
1917 params['num_procs'] = min(params['comm_world'].size, 5)
1918 params['plotting'] = params['comm_world'].rank == 0
1920 mode = f'vdp_stiffness-{mus[i]}'
1921 configurations = get_configs(mode=mode, problem=run_vdp)
1922 execute_configurations(**params, ax=axs.flatten()[i], decorate=True, configurations=configurations, mode=mode)
1923 axs.flatten()[i].set_title(rf'$\mu={ {mus[i]}} $')
1925 fig.suptitle('Van der Pol')
1926 if params['comm_world'].rank == 0:
1927 save_fig(
1928 fig=fig,
1929 name='vdp-stiffness',
1930 work_key=params['work_key'],
1931 precision_key=params['precision_key'],
1932 legend=False,
1933 base_path=base_path,
1934 format=format,
1935 )
1938def add_order_line(ax, order, ls, y_right=1.0, marker='.'):
1939 x_min = min([min(line.get_xdata()) for line in ax.get_lines()])
1940 x_max = max([max(line.get_xdata()) for line in ax.get_lines()])
1941 y_min = min([min(line.get_ydata()) for line in ax.get_lines()])
1942 y_max = max([max(line.get_ydata()) for line in ax.get_lines()])
1943 x = np.logspace(np.log10(x_min), np.log10(x_max), 100)
1944 y = y_right * (x / x_max) ** order
1945 mask = np.logical_and(y > y_min, y < y_max)
1946 ax.loglog(x[mask], y[mask], ls=ls, color='black', label=f'Order {order}', marker=marker, markevery=5)
1949def aggregate_parallel_efficiency_plot(): # pragma: no cover
1950 """
1951 Make a "weak scaling" plot for diagonal SDC
1952 """
1953 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
1955 fig, axs = plt.subplots(2, 2)
1957 _fig, _ax = plt.subplots(1, 1)
1958 num_procs = 1
1959 num_procs_sweeper = 2
1960 problem = run_quench
1962 num_procs_sweeper_list = [2, 3, 4]
1964 for problem, ax in zip([run_vdp, run_Lorenz, run_quench], axs.flatten()):
1965 speedup = []
1966 for num_procs_sweeper in num_procs_sweeper_list:
1967 s, e = plot_parallel_efficiency_diagonalSDC(
1968 ax=_ax,
1969 work_key='t',
1970 precision_key='e_global_rel',
1971 num_procs=num_procs,
1972 num_procs_sweeper=num_procs_sweeper,
1973 problem=problem,
1974 strategy=AdaptivityPolynomialError(),
1975 mode='diagonal_SDC',
1976 handle=f'{num_procs_sweeper} nodes',
1977 )
1978 speedup += [s]
1979 decorate_panel(ax, problem, work_key='nprocs', precision_key='')
1981 ax.plot(num_procs_sweeper_list, speedup, label='speedup')
1982 ax.plot(
1983 num_procs_sweeper_list,
1984 [speedup[i] / num_procs_sweeper_list[i] for i in range(len(speedup))],
1985 label='parallel efficiency',
1986 )
1988 fig.tight_layout()
1989 save_fig(fig, 'parallel_efficiency', 'nprocs', 'speedup')
1992if __name__ == "__main__":
1993 comm_world = MPI.COMM_WORLD
1995 import argparse
1997 parser = argparse.ArgumentParser()
1998 parser.add_argument('--mode', type=str, default='compare_strategies')
1999 parser.add_argument('--record', type=str, choices=['True', 'False'], default='True')
2000 parser.add_argument('--plotting', type=str, choices=['True', 'False'], default='True')
2001 parser.add_argument('--runs', type=int, default=5)
2002 parser.add_argument(
2003 '--problem', type=str, choices=['vdp', 'RBC', 'AC', 'quench', 'Lorenz', 'Schroedinger', 'GS'], default='vdp'
2004 )
2005 parser.add_argument('--work_key', type=str, default='t')
2006 parser.add_argument('--precision_key', type=str, default='e_global_rel')
2007 parser.add_argument('--logger_level', type=int, default='25')
2009 args = parser.parse_args()
2011 problems = {
2012 'Lorenz': run_Lorenz,
2013 'vdp': run_vdp,
2014 'Schroedinger': run_Schroedinger,
2015 'quench': run_quench,
2016 'AC': run_AC,
2017 'RBC': run_RBC,
2018 'GS': run_GS,
2019 }
2021 params = {
2022 **vars(args),
2023 'record': args.record == 'True',
2024 'plotting': args.plotting == 'True' and comm_world.rank == 0,
2025 'problem': problems[args.problem],
2026 }
2028 LOGGER_LEVEL = params.pop('logger_level')
2030 single_problem(**params)
2032 if comm_world.rank == 0:
2033 plt.show()