Coverage for pySDC/tutorial/step_3/A_getting_statistics.py: 100%
42 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 A: Getting statistics
12#
13# We run the heat equation from [Step 2](../step_2/C_using_pySDCs_frontend) again, now for eight time steps, and
14# look at what the controller returns besides the solution: the **statistics**. pySDC records residuals, iteration
15# counts, timings and more while it runs, and this is how we find out what happened.
17# %%
18import matplotlib.pyplot as plt
19import numpy as np
21from pySDC.helpers.stats_helper import get_list_of_types, get_sorted
22from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
23from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_forced
24from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
26# initialize level parameters
27level_params = {'restol': 1e-10, 'dt': 0.1}
29# initialize sweeper parameters
30sweeper_params = {'quad_type': 'RADAU-RIGHT', 'num_nodes': 3}
32# initialize problem parameters
33problem_params = {
34 'nu': 0.1, # diffusion coefficient
35 'freq': 4, # frequency for the test value
36 'nvars': 1023, # number of degrees of freedom
37 'bc': 'dirichlet-zero', # boundary conditions
38}
40# initialize step parameters
41step_params = {'maxiter': 20}
43# initialize controller parameters (<-- this is new!)
44controller_params = {'logger_level': 30} # reduce verbosity of each run
46# Fill description dictionary for easy hierarchy creation
47description = {
48 'problem_class': heatNd_forced,
49 'problem_params': problem_params,
50 'sweeper_class': imex_1st_order,
51 'sweeper_params': sweeper_params,
52 'level_params': level_params,
53 'step_params': step_params,
54}
56# instantiate the controller
57controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
59# set time parameters
60t0 = 0.1
61Tend = 0.9
63# get initial values on finest level
64P = controller.MS[0].levels[0].prob
65uinit = P.u_exact(t0)
67# call main function to get things done...
68uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
70# %% [markdown]
71# The controller parameter `logger_level` controls how much the run prints. It takes the
72# [levels of Python's logging module](https://docs.python.org/3/library/logging.html#logging-levels): 20 (`INFO`,
73# the default) shows every iteration, 30 (`WARNING`) only what goes wrong. That is why this run was silent.
74#
75# ## What is in `stats`
76#
77# `stats` is a dictionary. Each value is one number, and its key says where and when that number was recorded:
79# %%
80key = next(iter(stats))
81print(f'{len(stats)} entries, for example')
82print(f' key: {key}')
83print(f' value: {stats[key]}')
85# %% [markdown]
86# The keys are tuples that record the process, the time, the level, the iteration, the sweep, the type of the
87# value and more, so every value carries a kind of time stamp. This makes `stats` a little complex, but it also
88# means that anyone can add their own entries without changing its definition, as we do in
89# [Part B](B_adding_statistics). Which types a run recorded, including those added by users, `get_list_of_types`
90# tells:
92# %%
93print('List of registered statistic types:', get_list_of_types(stats))
95# %% [markdown]
96# The entries are made by **hooks**, see the `Hooks` class in `pySDC/core/hooks.py` and the default ones in
97# `pySDC/core/default_hook.py`.
98#
99# ## Filtering and sorting
100#
101# `get_sorted` picks the entries that match the keywords we give it, and returns them sorted by one field of the key,
102# as a list of tuples (value of that field, value). Here: all residuals recorded in the first time step, sorted by
103# iteration.
105# %%
106# filter statistics by first time interval and type (residual)
107residuals = get_sorted(stats, time=0.1, type='residual_post_iteration', sortby='iter')
109for item in residuals:
110 print('Residual in iteration %2i: %8.4e' % item)
112# %% [markdown]
113# And the number of iterations for every time step, sorted by time. `get_sorted` is a shortcut for
114# `sort_stats(filter_stats(...))`, with the same arguments.
116# %%
117# get and convert filtered statistics to list of iterations count, sorted by time
118iter_counts = get_sorted(stats, type='niter', sortby='time')
120for item in iter_counts:
121 print('Number of iterations at time %4.2f: %2i' % item)
123# %% [markdown]
124# :::{warning}
125# Filtering compares exactly. The first step starts at `t0 = 0.1`, which is stored as given, but later times are
126# sums of time steps: the third step starts at `0.30000000000000004`, so `time=0.3` finds nothing. To pick a later
127# step, sort by time and compare with a tolerance instead:
128# :::
130# %%
131print('filtered with time=0.3:', get_sorted(stats, time=0.3, type='niter'))
132print('compared with np.isclose:', [item for item in iter_counts if np.isclose(item[0], 0.3)])
134# %% [markdown]
135# With a little more filtering we can see all eight steps converge the same way:
137# %% tags=["hide-input"]
138fig, ax = plt.subplots(figsize=(6, 3.5))
139for time, _ in iter_counts:
140 # the exact time of each step, straight from the keys, so no filtering by computed times
141 residual = sorted(
142 (key.iter, value) for key, value in stats.items() if key.type == 'residual_post_iteration' and key.time == time
143 )
144 ax.semilogy(*zip(*residual, strict=True), label=f'$t = {time:.1f}$')
145ax.axhline(level_params['restol'], color='k', ls='--')
146ax.set_xlabel('iteration')
147ax.set_ylabel('residual')
148ax.legend(frameon=False, ncol=2, fontsize=8)
149ax.grid(alpha=0.3)
150fig.tight_layout()
152# %% [markdown]
153# ## Summary
154#
155# - The controller returns `stats`, a dictionary whose keys record where and when each value was taken.
156# - `get_list_of_types` shows what is in there; `get_sorted` filters and sorts it.
157# - Filter by time only with times you did not compute; otherwise, compare with a tolerance.
158#
159# The check the tests run:
161# %%
162assert all(item[1] == 12 for item in iter_counts), f'ERROR: number of iterations are not as expected, got {iter_counts}'