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_paramsto control the logging verbosity, see here for more details.Admittedly, the
statsdictionary 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
statsare created, please check theHooksclass.We also use the
get_list_of_typefunction to show what kind of values are registered in the statistics. This can be helpful, if users register their own types, see below.
Full code: pySDC/tutorial/step_3/A_getting_statistics.py
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['quad_type'] = 'RADAU-RIGHT'
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: ['residual_post_sweep', 'residual_post_iteration', 'niter', 'residual_post_step', '_recomputed', 'timing_setup', 'timing_comm', 'timing_sweep', 'timing_iteration', 'timing_step', '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.5260e-08
Residual in iteration 9: 6.0725e-09
Residual in iteration 10: 8.2734e-10
Residual in iteration 11: 1.1893e-10
Residual in iteration 12: 1.4850e-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
Part B: Adding statistics¶
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
supereach 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.
Full code: pySDC/tutorial/step_3/B_adding_statistics.py
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['quad_type'] = 'RADAU-RIGHT'
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:
2024-12-20 14:25:22,294 - controller - controller - welcome_message - 147 - INFO: Welcome to the one and only, really very astonishing and 87.3% bug free
_____ _____ _____
/ ____| __ \ / ____|
_ __ _ _| (___ | | | | |
| '_ \| | | |\___ \| | | | |
| |_) | |_| |____) | |__| | |____
| .__/ \__, |_____/|_____/ \_____|
| | __/ |
|_| |___/
2024-12-20 14:25:22,294 - controller - controller - dump_setup - 161 - INFO: Setup overview (--> user-defined, -> dependency) -- BEGIN
2024-12-20 14:25:22,294 - controller - controller - dump_setup - 228 - INFO: ----------------------------------------------------------------------------------------------------
Controller: <class 'pySDC.implementations.controller_classes.controller_nonMPI.controller_nonMPI'>
all_to_done = False
dump_setup = True
--> fname = data/step_3_B_out.txt
--> hook_class = [<class 'pySDC.implementations.hooks.default_hook.DefaultHooks'>, <class 'pySDC.implementations.hooks.log_timings.CPUTimings'>, <class 'pySDC.tutorial.step_3.HookClass_Particles.particle_hook'>]
--> log_to_file = True
logger_level = 20
mssdc_jac = True
predict_type = None
use_iteration_estimator = False
Step: <class 'pySDC.core.step.Step'>
--> maxiter = 20
Number of steps: None
Level: <class 'pySDC.core.level.Level'>
Level 0
--> dt = 0.0625
dt_initial = 0.0625
nsweeps = 1
residual_type = full_abs
--> restol = 1e-08
--> Problem: <class 'pySDC.implementations.problem_classes.PenningTrap_3D.penningtrap'>
--> nparts = 1
--> omega_B = 25.0
--> omega_E = 4.9
--> sig = 0.1
--> u0 = [list([10, 0, 0]) list([100, 0, 100]) list([1]) list([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'>
QE = EE
QI = IE
do_coll_update = False
initial_guess = spread
--> num_nodes = 3
--> quad_type = RADAU-RIGHT
skip_residual_computation = ()
--> Collocation: <class 'pySDC.core.collocation.CollBase'>
Active convergence controllers:
| # | order | convergence controller
----+----+-------+---------------------------------------------------------------------------------------
| 0 | 95 | BasicRestartingNonMPI
-> | 1 | 100 | SpreadStepSizesBlockwiseNonMPI
| 2 | 200 | CheckConvergence
2024-12-20 14:25:22,294 - controller - controller - dump_setup - 231 - INFO: ----------------------------------------------------------------------------------------------------
2024-12-20 14:25:22,294 - controller - controller - dump_setup - 233 - INFO: Setup overview (--> user-defined, -> dependency) -- END
2024-12-20 14:25:23,810 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 1 -- Sweep: 1 -- residual: 3.53203678e+00
2024-12-20 14:25:23,817 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 2 -- Sweep: 1 -- residual: 2.09852117e-01
2024-12-20 14:25:23,824 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 3 -- Sweep: 1 -- residual: 3.50301513e-02
2024-12-20 14:25:23,831 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 4 -- Sweep: 1 -- residual: 4.67724741e-03
2024-12-20 14:25:23,838 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 5 -- Sweep: 1 -- residual: 7.95583202e-04
2024-12-20 14:25:23,845 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 6 -- Sweep: 1 -- residual: 1.11405073e-04
2024-12-20 14:25:23,852 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 7 -- Sweep: 1 -- residual: 1.26902403e-05
2024-12-20 14:25:23,859 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 8 -- Sweep: 1 -- residual: 1.16534547e-06
2024-12-20 14:25:23,866 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 9 -- Sweep: 1 -- residual: 1.66968022e-07
2024-12-20 14:25:23,873 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 10 -- Sweep: 1 -- residual: 2.09408171e-08
2024-12-20 14:25:23,880 - hooks - default_hook - post_sweep - 22 - INFO: Process 0 on time 0.000000 at stage IT_FINE: Level: 0 -- Iteration: 11 -- Sweep: 1 -- residual: 2.17123386e-09
2024-12-20 14:25:23,883 - hooks - log_timings - post_run - 290 - INFO: Finished run after 1.59e+00s
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
statsdictionaries is not straightforward. Yet, putting them into a meta-dictionary is useful (as done here). Alternatively, eachstatscan be processed on the fly and only the relevant information can be stored.
Full code: pySDC/tutorial/step_3/C_study_collocations.py
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['RADAU-RIGHT'] = 15
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
quad_types = ['RADAU-RIGHT', 'GAUSS', 'LOBATTO']
stats_dict = dict()
for qtype in quad_types:
sweeper_params['quad_type'] = qtype
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.33862938e-05
Energy deviation for LOBATTO: 9.32710645e-06