Coverage for pySDC / projects / Resilience / work_precision.py: 0%
486 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-02-12 11:13 +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'], ['--', '-'], strict=True
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], ['-', '--', ':', '-.'], strict=True):
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 [
1366 9991,
1367 # 12123127391, 1231723109247102731092
1368 ],
1369 strict=True,
1370 ):
1371 configurations[i] = {
1372 'custom_description': {
1373 'sweeper_params': {'num_nodes': 3, 'QI': QI},
1374 'problem_params': desc["problem_params"],
1375 'sweeper_class': parallel_sweeper,
1376 },
1377 'strategies': [
1378 AdaptivityPolynomialError(
1379 useMPI=True, newton_inexactness=False, linear_inexactness=False, max_slope=4
1380 )
1381 ],
1382 'num_procs_sweeper': 3,
1383 'num_procs': 1,
1384 'plotting_params': {
1385 'ls': ls.get(QI, '-'),
1386 'label': rf'$\Delta t$-$k$-adaptivity $N$={1}x3',
1387 },
1388 'handle': f'{mode}-{QI}',
1389 'Tend': Tend,
1390 }
1392 elif mode == 'inexactness':
1393 """
1394 Compare inexact SDC to exact SDC
1395 """
1396 from pySDC.projects.Resilience.strategies import (
1397 AdaptivityPolynomialError,
1398 )
1400 if problem.__name__ in ['run_Schroedinger']:
1401 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1402 else:
1403 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1404 generic_implicit_MPI as parallel_sweeper,
1405 )
1407 strategies = [
1408 AdaptivityPolynomialError,
1409 ]
1411 inexactness = {
1412 'newton_inexactness': True,
1413 'linear_inexactness': True,
1414 }
1415 no_inexactness = {
1416 'newton_inexactness': False,
1417 'linear_inexactness': False,
1418 'SDC_maxiter': 99,
1419 'use_restol_rel': False,
1420 }
1422 configurations[1] = {
1423 'custom_description': {'sweeper_class': parallel_sweeper},
1424 'strategies': [me(useMPI=True, **no_inexactness) for me in strategies],
1425 'num_procs_sweeper': 3,
1426 'handle': 'exact',
1427 'plotting_params': {'ls': '--'},
1428 }
1429 configurations[0] = {
1430 'custom_description': {'sweeper_class': parallel_sweeper},
1431 'strategies': [me(useMPI=True, **inexactness) for me in strategies],
1432 'handle': 'inexact',
1433 'num_procs_sweeper': 3,
1434 }
1435 elif mode == 'compare_adaptivity':
1436 """
1437 Compare various modes of adaptivity
1438 """
1439 # TODO: configurations not final!
1440 from pySDC.projects.Resilience.strategies import (
1441 # AdaptivityCollocationTypeStrategy,
1442 # AdaptivityCollocationRefinementStrategy,
1443 AdaptivityStrategy,
1444 # AdaptivityExtrapolationWithinQStrategy,
1445 ESDIRKStrategy,
1446 ARKStrategy,
1447 AdaptivityPolynomialError,
1448 )
1450 if problem.__name__ in ['run_Schroedinger']:
1451 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1452 else:
1453 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1454 generic_implicit_MPI as parallel_sweeper,
1455 )
1457 inexactness_params = {
1458 # 'double_adaptivity': True,
1459 'newton_inexactness': True,
1460 'linear_inexactness': True,
1461 }
1463 strategies = [
1464 AdaptivityPolynomialError,
1465 # AdaptivityCollocationTypeStrategy,
1466 # AdaptivityExtrapolationWithinQStrategy,
1467 ]
1469 # restol = None
1470 # for strategy in strategies:
1471 # strategy.restol = restol
1473 configurations[1] = {
1474 'custom_description': {'sweeper_class': parallel_sweeper},
1475 'strategies': [me(useMPI=True, **inexactness_params) for me in strategies],
1476 'handle': 'parallel',
1477 'num_procs_sweeper': 3,
1478 'plotting_params': {'ls': '-', 'label': '3 procs'},
1479 }
1480 configurations[2] = {
1481 'strategies': [me(useMPI=True, **inexactness_params) for me in strategies],
1482 'plotting_params': {'ls': '--'},
1483 }
1484 configurations[4] = {
1485 'custom_description': {'step_params': {'maxiter': 5}},
1486 'strategies': [AdaptivityStrategy(useMPI=True)],
1487 }
1489 desc_RK = {}
1490 configurations[-1] = {
1491 'strategies': [
1492 ARKStrategy(useMPI=True) if problem.__name__ == 'run_Schroedinger' else ESDIRKStrategy(useMPI=True),
1493 ],
1494 'num_procs': 1,
1495 'custom_description': desc_RK,
1496 }
1498 elif mode == 'preconditioners':
1499 """
1500 Compare different preconditioners
1501 """
1502 from pySDC.projects.Resilience.strategies import (
1503 AdaptivityStrategy,
1504 IterateStrategy,
1505 BaseStrategy,
1506 ESDIRKStrategy,
1507 ERKStrategy,
1508 AdaptivityPolynomialError,
1509 )
1511 inexacness = True
1512 strategies = [
1513 AdaptivityPolynomialError(
1514 useMPI=True, SDC_maxiter=29, newton_inexactness=inexacness, linear_inexactness=inexacness
1515 ),
1516 BaseStrategy(useMPI=True),
1517 ]
1519 desc = {}
1520 desc['sweeper_params'] = {
1521 'num_nodes': 3,
1522 }
1523 # desc['step_params'] = {'maxiter': 5}
1525 precons = ['IE', 'LU']
1526 ls = ['-.', '--', '-', ':']
1527 for i in range(len(precons) + 1):
1528 if i < len(precons):
1529 desc['sweeper_params']['QI'] = precons[i]
1530 handle = precons[i]
1531 else:
1532 handle = None
1533 configurations[i] = {
1534 'strategies': strategies,
1535 'custom_description': copy.deepcopy(desc),
1536 'handle': handle,
1537 'plotting_params': {'ls': ls[i]},
1538 }
1539 elif mode == 'RK_comp_high_order':
1540 """
1541 Compare higher order SDC than we can get with RKM to RKM
1542 """
1543 from pySDC.projects.Resilience.strategies import (
1544 AdaptivityStrategy,
1545 ERKStrategy,
1546 ESDIRKStrategy,
1547 ARKStrategy,
1548 AdaptivityPolynomialError,
1549 )
1551 if problem.__name__ in ['run_Schroedinger']:
1552 from pySDC.implementations.sweeper_classes.imex_1st_order_MPI import imex_1st_order_MPI as parallel_sweeper
1553 else:
1554 from pySDC.implementations.sweeper_classes.generic_implicit_MPI import (
1555 generic_implicit_MPI as parallel_sweeper,
1556 )
1558 desc = {}
1559 desc['sweeper_params'] = {'num_nodes': 4, 'QI': 'IE', 'QE': "EE"}
1560 desc['step_params'] = {'maxiter': 7}
1562 desc_poly = {}
1563 desc_poly['sweeper_class'] = parallel_sweeper
1564 desc_poly['sweeper_params'] = {'num_nodes': 4}
1566 ls = {
1567 1: '-',
1568 2: '--',
1569 3: '-.',
1570 4: ':',
1571 5: ':',
1572 }
1574 desc_RK = {}
1575 if problem.__name__ in ['run_Schroedinger']:
1576 desc_RK['problem_params'] = {'imex': True}
1578 configurations[3] = {
1579 'custom_description': desc_poly,
1580 'strategies': [AdaptivityPolynomialError(useMPI=True)],
1581 'num_procs': 1,
1582 'num_procs_sweeper': 4,
1583 }
1584 configurations[-1] = {
1585 'strategies': [
1586 ERKStrategy(useMPI=True),
1587 ARKStrategy(useMPI=True) if problem.__name__ in ['run_Schroedinger'] else ESDIRKStrategy(useMPI=True),
1588 ],
1589 'num_procs': 1,
1590 'custom_description': desc_RK,
1591 }
1593 configurations[2] = {
1594 'strategies': [AdaptivityStrategy(useMPI=True)],
1595 'custom_description': desc,
1596 'num_procs': 4,
1597 }
1598 elif mode == 'avoid_restarts':
1599 """
1600 Test how well avoiding restarts works.
1601 """
1602 from pySDC.projects.Resilience.strategies import (
1603 AdaptivityStrategy,
1604 AdaptivityAvoidRestartsStrategy,
1605 AdaptivityPolynomialStrategy,
1606 )
1608 desc = {'sweeper_params': {'QI': 'IE'}, 'step_params': {'maxiter': 3}}
1609 param_range = [1e-3, 1e-5]
1610 configurations[0] = {
1611 'strategies': [AdaptivityPolynomialStrategy(useMPI=True)],
1612 'plotting_params': {'ls': '--'},
1613 'custom_description': desc,
1614 'param_range': param_range,
1615 }
1616 configurations[1] = {
1617 'strategies': [AdaptivityAvoidRestartsStrategy(useMPI=True)],
1618 'plotting_params': {'ls': '-.'},
1619 'custom_description': desc,
1620 'param_range': param_range,
1621 }
1622 configurations[2] = {
1623 'strategies': [AdaptivityStrategy(useMPI=True)],
1624 'custom_description': desc,
1625 'param_range': param_range,
1626 }
1627 else:
1628 raise NotImplementedError(f'Don\'t know the mode "{mode}"!')
1630 return configurations
1633def get_fig(x=1, y=1, target='adaptivity', **kwargs): # pragma: no cover
1634 """
1635 Get a figure to plot in.
1637 Args:
1638 x (int): How many panels in horizontal direction you want
1639 y (int): How many panels in vertical direction you want
1640 target (str): Where the plot is supposed to end up
1642 Returns:
1643 matplotlib.pyplot.Figure
1644 """
1645 width = 1.0
1646 ratio = 1.0 if y == 2 else 0.5
1647 if target == 'adaptivity':
1648 journal = 'Springer_Numerical_Algorithms'
1649 elif target == 'thesis':
1650 journal = 'TUHH_thesis'
1651 elif target == 'talk':
1652 journal = 'JSC_beamer'
1653 else:
1654 raise NotImplementedError
1656 keyword_arguments = {
1657 'figsize': figsize_by_journal(journal, width, ratio),
1658 'layout': 'constrained',
1659 **kwargs,
1660 }
1661 return plt.subplots(y, x, **keyword_arguments)
1664def save_fig(
1665 fig, name, work_key, precision_key, legend=True, format='pdf', base_path='data', squares=True, ncols=None, **kwargs
1666): # pragma: no cover
1667 """
1668 Save a figure with a legend on the bottom.
1670 Args:
1671 fig (matplotlib.pyplot.Figure): Figure you want to save
1672 name (str): Name of the plot to put in the path
1673 work_key (str): The key in the recorded data you want on the x-axis
1674 precision_key (str): The key in the recorded data you want on the y-axis
1675 legend (bool): Put a legend or not
1676 format (str): Format to store the figure with
1677 base_path (str): Some path where all the files are stored
1678 squares (bool): Adjust aspect ratio to squares if true
1680 Returns:
1681 None
1682 """
1683 handles = []
1684 labels = []
1685 for ax in fig.get_axes():
1686 h, l = ax.get_legend_handles_labels()
1687 handles += [h[i] for i in range(len(h)) if l[i] not in labels]
1688 labels += [me for me in l if me not in labels]
1689 if squares:
1690 ax.set_box_aspect(1)
1691 # order = np.argsort([me[0] for me in labels])
1692 order = np.arange(len(labels))
1693 fig.legend(
1694 [handles[i] for i in order],
1695 [labels[i] for i in order],
1696 loc='outside lower center',
1697 ncols=ncols if ncols else 3 if len(handles) % 3 == 0 else 4,
1698 frameon=False,
1699 fancybox=True,
1700 handlelength=2.2,
1701 )
1703 path = f'{base_path}/wp-{name}-{work_key}-{precision_key}.{format}'
1704 fig.savefig(path, bbox_inches='tight', **kwargs)
1705 print(f'Stored figure \"{path}\"')
1708def all_problems(
1709 mode='compare_strategies', plotting=True, base_path='data', target='adaptivity', **kwargs
1710): # pragma: no cover
1711 """
1712 Make a plot comparing various strategies for all problems.
1714 Args:
1715 work_key (str): The key in the recorded data you want on the x-axis
1716 precision_key (str): The key in the recorded data you want on the y-axis
1718 Returns:
1719 None
1720 """
1722 if target == 'talk':
1723 fig, axs = get_fig(4, 1, target=target)
1724 else:
1725 fig, axs = get_fig(2, 2, target=target)
1727 shared_params = {
1728 'work_key': 'k_SDC',
1729 'precision_key': 'e_global',
1730 'num_procs': 1,
1731 'runs': 1,
1732 'comm_world': MPI.COMM_WORLD,
1733 'record': False,
1734 'plotting': plotting,
1735 **kwargs,
1736 }
1738 if target == 'adaptivity':
1739 problems = [run_vdp, run_quench, run_Schroedinger, run_AC]
1740 elif target in ['thesis', 'talk']:
1741 problems = [run_vdp, run_Lorenz, run_GS, run_RBC]
1742 else:
1743 raise NotImplementedError
1745 logger.log(26, f"Doing for all problems {mode}")
1746 for i in range(len(problems)):
1747 execute_configurations(
1748 **shared_params,
1749 problem=problems[i],
1750 ax=axs.flatten()[i],
1751 decorate=True,
1752 configurations=get_configs(mode, problems[i]),
1753 mode=mode,
1754 )
1756 if plotting and shared_params['comm_world'].rank == 0:
1757 ncols = {
1758 'parallel_efficiency': 2,
1759 'parallel_efficiency_dt': 2,
1760 'parallel_efficiency_dt_k': 2,
1761 'RK_comp': 2,
1762 }
1763 if target == 'talk':
1764 _ncols = 4
1765 else:
1766 _ncols = ncols.get(mode, None)
1768 if shared_params['work_key'] == 'param':
1769 for ax, prob in zip(fig.get_axes(), problems, strict=True):
1770 add_param_order_lines(ax, prob)
1771 save_fig(
1772 fig=fig,
1773 name=mode,
1774 work_key=shared_params['work_key'],
1775 precision_key=shared_params['precision_key'],
1776 legend=True,
1777 base_path=base_path,
1778 ncols=_ncols,
1779 )
1782def add_param_order_lines(ax, problem):
1783 if problem.__name__ == 'run_vdp':
1784 yRfixed = 1e18
1785 yRdt = 1e-1
1786 yRdtk = 1e-4
1787 elif problem.__name__ == 'run_quench':
1788 yRfixed = 4e1
1789 yRdt = 1e4
1790 yRdtk = 1e4
1791 elif problem.__name__ == 'run_Schroedinger':
1792 yRfixed = 5
1793 yRdt = 1
1794 yRdtk = 1e-2
1795 elif problem.__name__ == 'run_AC':
1796 yRfixed = 1e8
1797 yRdt = 2e-2
1798 yRdtk = 1e-3
1799 elif problem.__name__ == 'run_Lorenz':
1800 yRfixed = 1e1
1801 yRdt = 2e-2
1802 yRdtk = 7e-4
1803 elif problem.__name__ == 'run_RBC':
1804 yRfixed = 1e-6
1805 yRdt = 4e-5
1806 yRdtk = 8e-6
1807 elif problem.__name__ == 'run_GS':
1808 yRfixed = 4e-3
1809 yRdt = 5e0
1810 yRdtk = 8e-1
1811 else:
1812 return None
1813 add_order_line(ax, 1, '--', yRdt, marker=None)
1814 add_order_line(ax, 5 / 4, ':', yRdtk, marker=None)
1815 add_order_line(ax, 5, '-.', yRfixed, marker=None)
1818def ODEs(mode='compare_strategies', plotting=True, base_path='data', **kwargs): # pragma: no cover
1819 """
1820 Make a plot comparing various strategies for the two ODEs.
1822 Args:
1823 work_key (str): The key in the recorded data you want on the x-axis
1824 precision_key (str): The key in the recorded data you want on the y-axis
1826 Returns:
1827 None
1828 """
1830 fig, axs = get_fig(x=2, y=1)
1832 shared_params = {
1833 'work_key': 'k_SDC',
1834 'precision_key': 'e_global',
1835 'num_procs': 1,
1836 'runs': 1,
1837 'comm_world': MPI.COMM_WORLD,
1838 'record': False,
1839 'plotting': plotting,
1840 **kwargs,
1841 }
1843 problems = [run_vdp, run_Lorenz]
1845 for i in range(len(problems)):
1846 execute_configurations(
1847 **shared_params,
1848 problem=problems[i],
1849 ax=axs.flatten()[i],
1850 decorate=i == 0,
1851 configurations=get_configs(mode, problems[i]),
1852 )
1854 if plotting and shared_params['comm_world'].rank == 0:
1855 save_fig(
1856 fig=fig,
1857 name=f'ODEs-{mode}',
1858 work_key=shared_params['work_key'],
1859 precision_key=shared_params['precision_key'],
1860 legend=True,
1861 base_path=base_path,
1862 )
1865def single_problem(mode, problem, plotting=True, base_path='data', target='thesis', **kwargs): # pragma: no cover
1866 """
1867 Make a plot for a single problem
1869 Args:
1870 mode (str): What you want to look at
1871 problem (function): A problem to run
1872 """
1873 if target == 'thesis':
1874 fig, ax = get_fig(1, 1, figsize=figsize_by_journal('TUHH_thesis', 0.7, 0.6))
1875 else:
1876 fig, ax = get_fig(1, 1, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1, 0.8))
1878 params = {
1879 'work_key': 'k_SDC',
1880 'precision_key': 'e_global',
1881 'num_procs': 1,
1882 'runs': 1,
1883 'comm_world': MPI.COMM_WORLD,
1884 'record': False,
1885 'plotting': plotting,
1886 **kwargs,
1887 }
1889 logger.log(26, f"Doing single problem {mode}")
1890 execute_configurations(
1891 **params, problem=problem, ax=ax, decorate=True, configurations=get_configs(mode, problem), mode=mode
1892 )
1894 if plotting:
1895 save_fig(
1896 fig=fig,
1897 name=f'{problem.__name__}-{mode}',
1898 work_key=params['work_key'],
1899 precision_key=params['precision_key'],
1900 legend=False,
1901 base_path=base_path,
1902 squares=target != 'thesis',
1903 )
1906def vdp_stiffness_plot(base_path='data', format='pdf', **kwargs): # pragma: no cover
1907 fig, axs = get_fig(3, 1, sharex=False, sharey=False)
1909 mus = [10, 100, 1000]
1911 for i in range(len(mus)):
1912 params = {
1913 'runs': 1,
1914 'problem': run_vdp,
1915 'record': False,
1916 'work_key': 't',
1917 'precision_key': 'e_global',
1918 'comm_world': MPI.COMM_WORLD,
1919 **kwargs,
1920 }
1921 params['num_procs'] = min(params['comm_world'].size, 5)
1922 params['plotting'] = params['comm_world'].rank == 0
1924 mode = f'vdp_stiffness-{mus[i]}'
1925 configurations = get_configs(mode=mode, problem=run_vdp)
1926 execute_configurations(**params, ax=axs.flatten()[i], decorate=True, configurations=configurations, mode=mode)
1927 axs.flatten()[i].set_title(rf'$\mu={{{mus[i]}}}$')
1929 fig.suptitle('Van der Pol')
1930 if params['comm_world'].rank == 0:
1931 save_fig(
1932 fig=fig,
1933 name='vdp-stiffness',
1934 work_key=params['work_key'],
1935 precision_key=params['precision_key'],
1936 legend=False,
1937 base_path=base_path,
1938 format=format,
1939 )
1942def add_order_line(ax, order, ls, y_right=1.0, marker='.'):
1943 x_min = min([min(line.get_xdata()) for line in ax.get_lines()])
1944 x_max = max([max(line.get_xdata()) for line in ax.get_lines()])
1945 y_min = min([min(line.get_ydata()) for line in ax.get_lines()])
1946 y_max = max([max(line.get_ydata()) for line in ax.get_lines()])
1947 x = np.logspace(np.log10(x_min), np.log10(x_max), 100)
1948 y = y_right * (x / x_max) ** order
1949 mask = np.logical_and(y > y_min, y < y_max)
1950 ax.loglog(x[mask], y[mask], ls=ls, color='black', label=f'Order {order}', marker=marker, markevery=5)
1953def aggregate_parallel_efficiency_plot(): # pragma: no cover
1954 """
1955 Make a "weak scaling" plot for diagonal SDC
1956 """
1957 from pySDC.projects.Resilience.strategies import AdaptivityPolynomialError
1959 fig, axs = plt.subplots(2, 2)
1961 _fig, _ax = plt.subplots(1, 1)
1962 num_procs = 1
1963 num_procs_sweeper = 2
1964 problem = run_quench
1966 num_procs_sweeper_list = [2, 3, 4]
1968 for problem, ax in zip([run_vdp, run_Lorenz, run_quench], axs.flatten(), strict=True):
1969 speedup = []
1970 for num_procs_sweeper in num_procs_sweeper_list:
1971 s, e = plot_parallel_efficiency_diagonalSDC(
1972 ax=_ax,
1973 work_key='t',
1974 precision_key='e_global_rel',
1975 num_procs=num_procs,
1976 num_procs_sweeper=num_procs_sweeper,
1977 problem=problem,
1978 strategy=AdaptivityPolynomialError(),
1979 mode='diagonal_SDC',
1980 handle=f'{num_procs_sweeper} nodes',
1981 )
1982 speedup += [s]
1983 decorate_panel(ax, problem, work_key='nprocs', precision_key='')
1985 ax.plot(num_procs_sweeper_list, speedup, label='speedup')
1986 ax.plot(
1987 num_procs_sweeper_list,
1988 [speedup[i] / num_procs_sweeper_list[i] for i in range(len(speedup))],
1989 label='parallel efficiency',
1990 )
1992 fig.tight_layout()
1993 save_fig(fig, 'parallel_efficiency', 'nprocs', 'speedup')
1996if __name__ == "__main__":
1997 comm_world = MPI.COMM_WORLD
1999 import argparse
2001 parser = argparse.ArgumentParser()
2002 parser.add_argument('--mode', type=str, default='compare_strategies')
2003 parser.add_argument('--record', type=str, choices=['True', 'False'], default='True')
2004 parser.add_argument('--plotting', type=str, choices=['True', 'False'], default='True')
2005 parser.add_argument('--runs', type=int, default=5)
2006 parser.add_argument(
2007 '--problem', type=str, choices=['vdp', 'RBC', 'AC', 'quench', 'Lorenz', 'Schroedinger', 'GS'], default='vdp'
2008 )
2009 parser.add_argument('--work_key', type=str, default='t')
2010 parser.add_argument('--precision_key', type=str, default='e_global_rel')
2011 parser.add_argument('--logger_level', type=int, default='25')
2013 args = parser.parse_args()
2015 problems = {
2016 'Lorenz': run_Lorenz,
2017 'vdp': run_vdp,
2018 'Schroedinger': run_Schroedinger,
2019 'quench': run_quench,
2020 'AC': run_AC,
2021 'RBC': run_RBC,
2022 'GS': run_GS,
2023 }
2025 params = {
2026 **vars(args),
2027 'record': args.record == 'True',
2028 'plotting': args.plotting == 'True' and comm_world.rank == 0,
2029 'problem': problems[args.problem],
2030 }
2032 LOGGER_LEVEL = params.pop('logger_level')
2034 single_problem(**params)
2036 if comm_world.rank == 0:
2037 plt.show()