Coverage for pySDC/tutorial/step_5/B_my_first_PFASST_run.py: 100%
50 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1# ---
2# jupyter:
3# jupytext:
4# formats: py:percent
5# kernelspec:
6# display_name: Python 3
7# name: python3
8# ---
10# %% [markdown]
11# # Part B: My first PFASST run
12#
13# With the multistep multilevel hierarchy of [Part A](A_multistep_multilevel_hierarchy), we are ready for PFASST,
14# starting with the heat equation. One of the most important characteristics of a parallel-in-time method is how it
15# behaves when more time steps are computed in parallel for the same problem, i.e. with `dt` and `Tend` fixed. So we
16# run the 16 time steps of $[0, 4]$ with 1, 2, 4, 8 and 16 parallel steps.
18# %%
19import matplotlib.pyplot as plt
20import numpy as np
22from pySDC.helpers.stats_helper import get_sorted
23from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
24from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_forced
25from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
26from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
28# initialize level parameters
29level_params = {'restol': 1e-10, 'dt': 0.25}
31# initialize sweeper parameters
32sweeper_params = {
33 'quad_type': 'RADAU-RIGHT',
34 'num_nodes': [3],
35 'QI': 'LU', # For the IMEX sweeper, the LU-trick can be activated for the implicit part
36}
38# initialize problem parameters
39problem_params = {
40 'nu': 0.1, # diffusion coefficient
41 'freq': 8, # frequency for the test value
42 'nvars': [511, 255], # number of degrees of freedom for each level
43 'bc': 'dirichlet-zero', # boundary conditions
44}
46# initialize step parameters
47step_params = {'maxiter': 50}
49# initialize space transfer parameters
50space_transfer_params = {'rorder': 2, 'iorder': 6}
52# initialize controller parameters
53controller_params = {
54 'logger_level': 30,
55 'predict_type': 'pfasst_burnin', # PFASST's coarse-level predictor, see step 4, part C
56}
58# fill description dictionary for easy step instantiation
59description = {
60 'problem_class': heatNd_forced, # pass problem class
61 'problem_params': problem_params, # pass problem parameters
62 'sweeper_class': imex_1st_order, # pass sweeper
63 'sweeper_params': sweeper_params, # pass sweeper parameters
64 'level_params': level_params, # pass level parameters
65 'step_params': step_params, # pass step parameters
66 'space_transfer_class': mesh_to_mesh, # pass spatial transfer class
67 'space_transfer_params': space_transfer_params, # pass parameters for spatial transfer
68}
70# set time parameters
71t0 = 0.0
72Tend = 4.0
74# set up list of parallel time-steps to run PFASST with
75nsteps = int(Tend / level_params['dt'])
76num_proc_list = [2**i for i in range(int(np.log2(nsteps) + 1))]
77print('numbers of parallel time steps:', num_proc_list)
79# %% [markdown]
80# ## The runs
81#
82# One controller for each number of parallel steps; the rest is what we did in [Step 3](../step_3), for each run:
83# the error, and statistics of the iteration counts.
85# %% tags=["scroll-output"]
86results = {}
87# loop over different number of processes and check results
88for num_proc in num_proc_list:
89 print('Working with %2i processes...' % num_proc)
90 # instantiate controller
91 controller = controller_nonMPI(num_procs=num_proc, controller_params=controller_params, description=description)
93 # get initial values on finest level
94 P = controller.MS[0].levels[0].prob
95 uinit = P.u_exact(t0)
97 # call main function to get things done...
98 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
100 # compute exact solution and compare
101 uex = P.u_exact(Tend)
102 err = abs(uex - uend)
104 # filter statistics by type (number of iterations)
105 iter_counts = get_sorted(stats, type='niter', sortby='time')
107 # compute and print statistics
108 for item in iter_counts:
109 print('Number of iterations for time %4.2f: %2i' % item)
110 print()
111 niters = np.array([item[1] for item in iter_counts])
112 print(' Mean number of iterations: %4.2f' % np.mean(niters))
113 print(' Range of values for number of iterations: %2i ' % np.ptp(niters))
114 print(' Position of max/min number of iterations: %2i -- %2i' % (int(np.argmax(niters)), int(np.argmin(niters))))
115 print(' Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters))))
116 print()
117 print()
119 results[num_proc] = (err, niters)
121# %% [markdown]
122# The iteration counts of all runs at a glance, one row per number of parallel steps:
124# %% tags=["hide-input"]
125fig, ax = plt.subplots(figsize=(8, 2.6), constrained_layout=True)
126image = ax.imshow([results[n][1] for n in num_proc_list], cmap='viridis', aspect='auto', vmin=0)
127ax.set_yticks(range(len(num_proc_list)))
128ax.set_yticklabels(num_proc_list)
129ax.set_ylabel('parallel steps')
130ax.set_xlabel('time step')
131fig.colorbar(image, label='iterations')
133# %% [markdown]
134# PFASST performs very well here: whether the 16 steps are computed one after another or all at once, they need
135# about 5 iterations each.
136#
137# :::{admonition} Important things to note
138# - In the IMEX sweeper, `QI = 'LU'` activates the LU-trick for the implicit part. For stiff parabolic problems
139# with Gauss-Radau nodes, this is usually a very good idea.
140# - As usual for MLSDC and PFASST, success depends heavily on the parameters. A more complicated problem, a more
141# or less stiff one, a different order of the spatial interpolation, etc. can give completely different results.
142# [Part C](C_advection_and_PFASST) has an example.
143# :::
144#
145# The checks the tests run, for every number of parallel steps:
147# %%
148for err, niters in results.values():
149 assert err < 1.3505e-04, f"ERROR: error is too high, got {err}"
150 assert np.ptp(niters) <= 1, f"ERROR: range of number of iterations is too high, got {np.ptp(niters)}"
151 assert np.mean(niters) <= 5.0, f"ERROR: mean number of iterations is too high, got {np.mean(niters)}"