# Step-3: Statistics and a new sweeper¶

In this step, we will show how to work with the statistics pySDC generates. We will also introduce a new problem as well as a new sweeper.

## Part A: Getting statistics¶

In the first part, we run our standard heat equation problem again, but focus on the `stats` dictionary the controller returns in addition to the final solution. This `stats` dictionary is a little bit complex, as its keys are tuples which contain the process, time, level, iteration and type of entry for the values. This way, each value has some sort of time stamp attached to it, so that it is clear when this value was added. To help dealing with the statistics, we make use of the `get_sorted` routine. This filters the dictionary, following the keys we provide. Here, we would like to have all residuals logged during time 0.1 (i.e. for all iterations in the first time step). Then, it is converted to a list of tuples, where the first part is the item defined by the parameter `sortby` and the second part is the value. Here, we would like to have a list of iterations and residuals to see how SDC converged over the iterations.

Important things to note:

• Here, we use the `controller_params` to control the logging verbosity, see here for more details.

• Admittedly, the `stats` dictionary is a complex thing, but it allows users to add other details to it without changing its definition (see next part).

• For more details on how the entries of `stats` are created, please check the `Hooks` class.

• We also use the `get_list_of_type` function to show what kind of values are registered in the statistics. This can be helpful, if users register their own types, see below.

```from pathlib import Path

from pySDC.helpers.stats_helper import get_list_of_types, get_sorted

from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_forced
from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order

def main():
"""
A simple test program to describe how to get statistics of a run
"""

# run simulation
stats = run_simulation()

Path("data").mkdir(parents=True, exist_ok=True)
f = open('data/step_3_A_out.txt', 'w')
out = 'List of registered statistic types: %s' % get_list_of_types(stats)
f.write(out + '\n')
print(out)

# filter statistics by first time interval and type (residual)
residuals = get_sorted(stats, time=0.1, type='residual_post_iteration', sortby='iter')

for item in residuals:
out = 'Residual in iteration %2i: %8.4e' % item
f.write(out + '\n')
print(out)

# get and convert filtered statistics to list of iterations count, sorted by time
# the get_sorted function is just a shortcut for sort_stats(filter_stats()) with all the same arguments
iter_counts = get_sorted(stats, type='niter', sortby='time')

for item in iter_counts:
out = 'Number of iterations at time %4.2f: %2i' % item
f.write(out + '\n')
print(out)

f.close()

assert all(item[1] == 12 for item in iter_counts), (
'ERROR: number of iterations are not as expected, got %s' % iter_counts
)

def run_simulation():
"""
A simple test program to run IMEX SDC for a single time step
"""
# initialize level parameters
level_params = {}
level_params['restol'] = 1e-10
level_params['dt'] = 0.1

# initialize sweeper parameters
sweeper_params = {}
sweeper_params['num_nodes'] = 3

# initialize problem parameters
problem_params = {}
problem_params['nu'] = 0.1  # diffusion coefficient
problem_params['freq'] = 4  # frequency for the test value
problem_params['nvars'] = 1023  # number of degrees of freedom
problem_params['bc'] = 'dirichlet-zero'  # boundary conditions

# initialize step parameters
step_params = {}
step_params['maxiter'] = 20

# initialize controller parameters (<-- this is new!)
controller_params = {}
controller_params['logger_level'] = 30  # reduce verbosity of each run

# Fill description dictionary for easy hierarchy creation
description = {}
description['problem_class'] = heatNd_forced
description['problem_params'] = problem_params
description['sweeper_class'] = imex_1st_order
description['sweeper_params'] = sweeper_params
description['level_params'] = level_params
description['step_params'] = step_params

# instantiate the controller (no controller parameters used here)
controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)

# set time parameters
t0 = 0.1
Tend = 0.9

# get initial values on finest level
P = controller.MS[0].levels[0].prob
uinit = P.u_exact(t0)

# call main function to get things done...
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)

return stats

if __name__ == "__main__":
main()
```

Results:

```List of registered statistic types: ['timing_setup', 'timing_comm', 'residual_post_sweep', 'timing_sweep', 'residual_post_iteration', 'timing_iteration', 'timing_step', 'niter', 'residual_post_step', '_recomputed', 'timing_run', 'restart']
Residual in iteration  1: 4.1119e-03
Residual in iteration  2: 6.6844e-04
Residual in iteration  3: 8.8038e-05
Residual in iteration  4: 1.2171e-05
Residual in iteration  5: 1.3827e-06
Residual in iteration  6: 6.3645e-07
Residual in iteration  7: 1.6895e-07
Residual in iteration  8: 3.5261e-08
Residual in iteration  9: 6.0742e-09
Residual in iteration 10: 8.2757e-10
Residual in iteration 11: 1.1908e-10
Residual in iteration 12: 1.4679e-11
Number of iterations at time 0.10: 12
Number of iterations at time 0.20: 12
Number of iterations at time 0.30: 12
Number of iterations at time 0.40: 12
Number of iterations at time 0.50: 12
Number of iterations at time 0.60: 12
Number of iterations at time 0.70: 12
Number of iterations at time 0.80: 12
```

We now extend the statistics of pySDC by user-defined entries. To make things much more interesting (and complicated), we introduce a new problem class (`PenningTrap_3D`) as well a new sweeper (`boris_2nd_order`). This is accompanied by new data types, namely `particles` and `fields`. Details on the idea of the Boris solver for particles moving in electromagnetic fields can be found in this paper. Yet, one important measure for this kind of problem is the total energy of the system. We would like to compute this after each time step and therefore, we define a new hook class called `particle_hooks`. Here, all necessary computations are done and the value is added to the statistics for later processing.

Important things to note:

• In order to extend (and not replace) pySDC’s statistics, make sure that your custom hook class calls `super` each time a function is called.

• User-defined statistics can also be used from within the problem class: simply define a problem attribute (e.g. the number of GMRES iterations in the spatial solver) and access it during the hook call using the level.

```import numpy as np
from pathlib import Path

from pySDC.helpers.stats_helper import get_sorted

from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
from pySDC.implementations.problem_classes.PenningTrap_3D import penningtrap
from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order
from pySDC.tutorial.step_3.HookClass_Particles import particle_hook

def main():
"""
A simple test program to retrieve user-defined statistics from a run
"""
Path("data").mkdir(parents=True, exist_ok=True)

err, stats = run_penning_trap_simulation()

# filter statistics type (etot)
energy = get_sorted(stats, type='etot', sortby='iter')

# get base energy and show difference
base_energy = energy[0][1]
f = open('data/step_3_B_out.txt', 'a')
for item in energy:
out = 'Total energy and deviation in iteration %2i: %12.10f -- %12.8e' % (
item[0],
item[1],
abs(base_energy - item[1]),
)
f.write(out + '\n')
print(out)
f.close()

assert abs(base_energy - energy[-1][1]) < 15, 'ERROR: energy deviated too much, got %s' % (
base_energy - energy[-1][1]
)
assert err < 5e-04, "ERROR: solution is not as exact as expected, got %s" % err

def run_penning_trap_simulation():
"""
A simple test program to run IMEX SDC for a single time step of the penning trap example
"""
# initialize level parameters
level_params = dict()
level_params['restol'] = 1e-08
level_params['dt'] = 1.0 / 16

# initialize sweeper parameters
sweeper_params = dict()
sweeper_params['num_nodes'] = 3

# initialize problem parameters for the Penning trap
problem_params = dict()
problem_params['omega_E'] = 4.9  # E-field frequency
problem_params['omega_B'] = 25.0  # B-field frequency
problem_params['u0'] = np.array([[10, 0, 0], [100, 0, 100], [1], [1]], dtype=object)  # initial center of positions
problem_params['nparts'] = 1  # number of particles in the trap
problem_params['sig'] = 0.1  # smoothing parameter for the forces

# initialize step parameters
step_params = dict()
step_params['maxiter'] = 20

# initialize controller parameters
controller_params = dict()
controller_params['hook_class'] = particle_hook  # specialized hook class for more statistics and output
controller_params['log_to_file'] = True
controller_params['fname'] = 'data/step_3_B_out.txt'

# Fill description dictionary for easy hierarchy creation
description = dict()
description['problem_class'] = penningtrap
description['problem_params'] = problem_params
description['sweeper_class'] = boris_2nd_order
description['sweeper_params'] = sweeper_params
description['level_params'] = level_params
description['step_params'] = step_params

# instantiate the controller (no controller parameters used here)
controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)

# set time parameters
t0 = 0.0
Tend = level_params['dt']

# get initial values on finest level
P = controller.MS[0].levels[0].prob
uinit = P.u_init()

# call main function to get things done...
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)

# compute error compared to know exact solution for one particle
uex = P.u_exact(Tend)
err = np.linalg.norm(uex.pos - uend.pos, np.inf) / np.linalg.norm(uex.pos, np.inf)

return err, stats

if __name__ == "__main__":
main()
```

Results:

```2023-09-24 09:59:39,812 - controller - Controller - welcome_message - 146 - INFO: Welcome to the one and only, really very astonishing and 87.3% bug free
_____ _____   _____
/ ____|  __ \ / ____|
_ __  _   _| (___ | |  | | |
| '_ \| | | |\___ \| |  | | |
| |_) | |_| |____) | |__| | |____
| .__/ \__, |_____/|_____/ \_____|
| |     __/ |
|_|    |___/

2023-09-24 09:59:39,812 - controller - Controller - dump_setup - 160 - INFO: Setup overview (--> user-defined, -> dependency) -- BEGIN
2023-09-24 09:59:39,813 - controller - Controller - dump_setup - 226 - INFO: ----------------------------------------------------------------------------------------------------

Controller: <class 'pySDC.implementations.controller_classes.controller_nonMPI.controller_nonMPI'>
mssdc_jac = True
predict_type = None
all_to_done = False
logger_level = 20
--> log_to_file = True
dump_setup = True
--> fname = data/step_3_B_out.txt
use_iteration_estimator = False
--> hook_class = <class 'pySDC.tutorial.step_3.HookClass_Particles.particle_hook'>

Step: <class 'pySDC.core.Step.step'>
--> maxiter = 20
Level: <class 'pySDC.core.Level.level'>
Level  0
-->         dt = 0.0625
dt_initial = 0.0625
-->         restol = 1e-08
nsweeps = 1
residual_type = full_abs
-->         Problem: <class 'pySDC.implementations.problem_classes.PenningTrap_3D.penningtrap'>
-->             omega_E = 4.9
-->             nparts = 1
-->             omega_B = 25.0
-->             u0 = [list([10, 0, 0]) list([100, 0, 100]) list([1]) list([1])]
-->             sig = 0.1
-->             Data type u: <class 'pySDC.implementations.datatype_classes.particles.particles'>
-->             Data type f: <class 'pySDC.implementations.datatype_classes.particles.fields'>
-->             Sweeper: <class 'pySDC.implementations.sweeper_classes.boris_2nd_order.boris_2nd_order'>
do_coll_update = False
skip_residual_computation = ()
-->                 num_nodes = 3
QI = IE
QE = EE
-->                 Collocation: <class 'pySDC.core.Collocation.CollBase'>

Active convergence controllers:
|  # | order | convergence controller
----+----+-------+---------------------------------------------------------------------------------------
|  0 |    95 | BasicRestartingNonMPI
-> |  1 |   100 | SpreadStepSizesBlockwiseNonMPI
|  2 |   200 | CheckConvergence

2023-09-24 09:59:39,813 - controller - Controller - dump_setup - 229 - INFO: ----------------------------------------------------------------------------------------------------
2023-09-24 09:59:39,813 - controller - Controller - dump_setup - 231 - INFO: Setup overview (--> user-defined, -> dependency) -- END

2023-09-24 09:59:39,822 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration:  1 -- Sweep:  1 -- residual: 3.53203678e+00
2023-09-24 09:59:39,830 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration:  2 -- Sweep:  1 -- residual: 2.09852117e-01
2023-09-24 09:59:39,838 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration:  3 -- Sweep:  1 -- residual: 3.50301513e-02
2023-09-24 09:59:39,846 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration:  4 -- Sweep:  1 -- residual: 4.67724741e-03
2023-09-24 09:59:39,854 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration:  5 -- Sweep:  1 -- residual: 7.95583202e-04
2023-09-24 09:59:39,863 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration:  6 -- Sweep:  1 -- residual: 1.11405073e-04
2023-09-24 09:59:39,871 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration:  7 -- Sweep:  1 -- residual: 1.26902403e-05
2023-09-24 09:59:39,879 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration:  8 -- Sweep:  1 -- residual: 1.16534547e-06
2023-09-24 09:59:39,888 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration:  9 -- Sweep:  1 -- residual: 1.66968007e-07
2023-09-24 09:59:39,896 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration: 10 -- Sweep:  1 -- residual: 2.09407887e-08
2023-09-24 09:59:39,904 - hooks - default_hook - post_sweep - 170 - INFO: Process  0 on time 0.000000 at stage         IT_FINE: Level: 0 -- Iteration: 11 -- Sweep:  1 -- residual: 2.17123386e-09
2023-09-24 09:59:39,908 - hooks - default_hook - post_run - 340 - INFO: Finished run after 0.09s
Total energy and deviation in iteration  0: 8799.5000000000 -- 0.00000000e+00
Total energy and deviation in iteration 11: 8785.0038936088 -- 1.44961064e+01
```

## Part C: Studying collocation node types¶

Having first experience with the Boris-SDC solver, we see that the energy deviates quite a bit already for this simple setup. We now test this for three different collocation classes to demonstrate how a parameter study can be done with pySDC. To this end, we describe the whole setup up to the parameter we would like to vary. Using a simple list of parameters, we construct the controller each time using a different collocation class. While slightly inefficient, this makes sure that the variables and statistics are reset each time.

Important things to note:

• Interestingly, the energy does not deviate a lot for Gauss-Lobatto and Gauss-Legendre nodes. This is related to their symmetric nature (in contrast to Gauss-Radau).

• Using a lower residual tolerance actually improves the energy conservation even further, at least for symmetric collocation nodes. Rule of thumb: conservation up to residual tolerance is achieved.

• Working with multiple `stats` dictionaries is not straightforward. Yet, putting them into a meta-dictionary is useful (as done here). Alternatively, each `stats` can be processed on the fly and only the relevant information can be stored.

```from pathlib import Path
import numpy as np

from pySDC.helpers.stats_helper import get_sorted
from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
from pySDC.implementations.problem_classes.PenningTrap_3D import penningtrap
from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order
from pySDC.tutorial.step_3.HookClass_Particles import particle_hook

def main():
"""
A simple test program to show the energy deviation for different quadrature nodes
"""
stats_dict = run_simulation()

ediff = dict()
Path("data").mkdir(parents=True, exist_ok=True)
f = open('data/step_3_C_out.txt', 'w')
for cclass, stats in stats_dict.items():
# filter and convert/sort statistics by etot and iterations
energy = get_sorted(stats, type='etot', sortby='iter')
# compare base and final energy
base_energy = energy[0][1]
final_energy = energy[-1][1]
ediff[cclass] = abs(base_energy - final_energy)
out = "Energy deviation for %s: %12.8e" % (cclass, ediff[cclass])
f.write(out + '\n')
print(out)
f.close()

# set expected differences and check
ediff_expect = dict()
ediff_expect['LOBATTO'] = 1e-05
ediff_expect['GAUSS'] = 3e-05
for k, v in ediff.items():
assert v < ediff_expect[k], "ERROR: energy deviated too much, got %s" % ediff[k]

def run_simulation():
"""
A simple test program to run IMEX SDC for a single time step
"""
# initialize level parameters
level_params = dict()
level_params['restol'] = 1e-06
level_params['dt'] = 1.0 / 16

# initialize sweeper parameters
sweeper_params = dict()
sweeper_params['num_nodes'] = 3

# initialize problem parameters
problem_params = dict()
problem_params['omega_E'] = 4.9
problem_params['omega_B'] = 25.0
problem_params['u0'] = np.array([[10, 0, 0], [100, 0, 100], [1], [1]], dtype=object)
problem_params['nparts'] = 1
problem_params['sig'] = 0.1

# initialize step parameters
step_params = dict()
step_params['maxiter'] = 20

# initialize controller parameters
controller_params = dict()
controller_params['hook_class'] = particle_hook  # specialized hook class for more statistics and output
controller_params['logger_level'] = 30  # reduce verbosity of each run

# Fill description dictionary for easy hierarchy creation
description = dict()
description['problem_class'] = penningtrap
description['problem_params'] = problem_params
description['sweeper_class'] = boris_2nd_order
description['level_params'] = level_params
description['step_params'] = step_params

# assemble and loop over list of collocation classes
stats_dict = dict()
description['sweeper_params'] = sweeper_params

# instantiate the controller (no controller parameters used here)
controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)

# set time parameters
t0 = 0.0
Tend = level_params['dt']

# get initial values on finest level
P = controller.MS[0].levels[0].prob
uinit = P.u_init()

# call main function to get things done...
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)

# gather stats in dictionary, collocation classes being the keys
stats_dict[qtype] = stats

return stats_dict

if __name__ == "__main__":
main()
```

Results:

```Energy deviation for RADAU-RIGHT: 1.44960920e+01
Energy deviation for GAUSS: 2.33862975e-05
Energy deviation for LOBATTO: 9.32710282e-06
```