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.

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 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.

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-11-16 14:26:15,548 - controller - controller - welcome_message - 147 - INFO: Welcome to the one and only, really very astonishing and 87.3% bug free
                                 _____ _____   _____ 
                                / ____|  __ \ / ____|
                    _ __  _   _| (___ | |  | | |     
                   | '_ \| | | |\___ \| |  | | |     
                   | |_) | |_| |____) | |__| | |____ 
                   | .__/ \__, |_____/|_____/ \_____|
                   | |     __/ |                     
                   |_|    |___/                      
                                                     
2024-11-16 14:26:15,548 - controller - controller - dump_setup - 161 - INFO: Setup overview (--> user-defined, -> dependency) -- BEGIN
2024-11-16 14:26:15,548 - 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-11-16 14:26:15,548 - controller - controller - dump_setup - 231 - INFO: ----------------------------------------------------------------------------------------------------
2024-11-16 14:26:15,548 - controller - controller - dump_setup - 233 - INFO: Setup overview (--> user-defined, -> dependency) -- END

2024-11-16 14:26:16,907 - 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-11-16 14:26:16,913 - 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-11-16 14:26:16,920 - 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-11-16 14:26:16,926 - 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-11-16 14:26:16,933 - 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-11-16 14:26:16,940 - 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-11-16 14:26:16,946 - 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-11-16 14:26:16,953 - 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-11-16 14:26:16,959 - 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-11-16 14:26:16,966 - 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-11-16 14:26:16,973 - 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-11-16 14:26:16,975 - hooks - log_timings - post_run - 290 - INFO: Finished run after 1.43e+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 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.

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