Step-4: Multilevel SDC¶

In this step, we will show how pySDC creates a multilevel hierarchy and how MLSDC can be run and tested.

Part A: Spatial transfer operators¶

For a multilevel hierarchy, we need transfer operators. The user, having knowledge of the data types, will have to provide a space_transfer_class which deals with restriction and interpolation in the spatial dimension. In this part, we simply set up two problems with two different resolutions and check the order of interpolation (4 in this case).

Important things to note:

As for the sweeper and many other things, the user does not have to deal with the instantiation of the transfer class, when one of the controllers is used. This is for demonstrational purpose only.
MLSDC (and PFASST) rely on high-order interpolation in space. When using Lagrange-based interpolation, orders above 4-6 are recommended. Restriction, however, can be of order 2 and is thus not tested here.

Full code: pySDC/tutorial/step_4/A_spatial_transfer_operators.py

from collections import namedtuple
from pathlib import Path

import numpy as np

from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
from pySDC.tutorial.step_1.B_spatial_accuracy_check import get_accuracy_order

# setup id for gathering the results (will sort by nvars)
ID = namedtuple('ID', 'nvars_fine')


def main():
    """
    A simple test program to test interpolation order in space
    """

    # initialize problem parameters
    problem_params = dict()
    problem_params['nu'] = 0.1  # diffusion coefficient
    problem_params['freq'] = 3  # frequency for the test value
    problem_params['bc'] = 'dirichlet-zero'  # boundary conditions

    # initialize transfer parameters
    space_transfer_params = dict()
    space_transfer_params['rorder'] = 2
    space_transfer_params['iorder'] = 4

    nvars_fine_list = [2**p - 1 for p in range(5, 10)]

    # set up dictionary to store results (plus lists)
    results = dict()
    results['nvars_list'] = nvars_fine_list

    for nvars_fine in nvars_fine_list:
        print('Working on nvars_fine = %4i...' % nvars_fine)

        # instantiate fine problem
        problem_params['nvars'] = nvars_fine  # number of degrees of freedom
        Pfine = heatNd_unforced(**problem_params)

        # instantiate coarse problem using half of the DOFs
        problem_params['nvars'] = int((nvars_fine + 1) / 2.0 - 1)
        Pcoarse = heatNd_unforced(**problem_params)

        # instantiate spatial interpolation
        T = mesh_to_mesh(fine_prob=Pfine, coarse_prob=Pcoarse, params=space_transfer_params)

        # set exact fine solution to compare with
        xvalues_fine = np.array([(i + 1) * Pfine.dx for i in range(Pfine.nvars[0])])
        uexact_fine = Pfine.dtype_u(Pfine.init)
        uexact_fine[:] = np.sin(np.pi * Pfine.freq[0] * xvalues_fine)

        # set exact coarse solution as source
        xvalues_coarse = np.array([(i + 1) * Pcoarse.dx for i in range(Pcoarse.nvars[0])])
        uexact_coarse = Pfine.dtype_u(Pcoarse.init)
        uexact_coarse[:] = np.sin(np.pi * Pcoarse.freq[0] * xvalues_coarse)

        # do the interpolation/prolongation
        uinter = T.prolong(uexact_coarse)

        # compute error and store
        id = ID(nvars_fine=nvars_fine)
        results[id] = abs(uinter - uexact_fine)

    # print out and check
    print('Running order checks...')
    orders = get_accuracy_order(results)
    Path("data").mkdir(parents=True, exist_ok=True)
    f = open('data/step_4_A_out.txt', 'w')
    for p in range(len(orders)):
        out = 'Expected order %2i, got order %5.2f, deviation of %5.2f%%' % (
            space_transfer_params['iorder'],
            orders[p],
            100 * abs(space_transfer_params['iorder'] - orders[p]) / space_transfer_params['iorder'],
        )
        f.write(out + '\n')
        print(out)
        assert (
            abs(space_transfer_params['iorder'] - orders[p]) / space_transfer_params['iorder'] < 0.05
        ), 'ERROR: did not get expected orders for interpolation, got %s' % str(orders[p])
    f.close()
    print('...got what we expected!')


if __name__ == "__main__":
    main()

Results:

Expected order  4, got order  4.08, deviation of  1.91%
Expected order  4, got order  3.95, deviation of  1.35%
Expected order  4, got order  3.98, deviation of  0.62%
Expected order  4, got order  3.99, deviation of  0.30%

Part B: Multilevel hierarchy¶

In this example, we demonstrate how the step class creates the space-time hierarchy dynamically, depending on the description and its parameters. pySDC supports two different generic coarsening strategies: coarsening in space and coarsening in the collocation order. To enable collocation-based coarsening, we simply replace the num_nodes parameter by a list, where the first entry corresponds to the finest level. For spatial coarsening, the problem parameter nvars is replaced by a list, too. During the step setup, these dictionaries with list entries are transformed into lists of dictionaries corresponding to the levels (3 in this case). A third generic way of creating multiple levels is to replace an entry in the description by a list, e.g. a list of problem classes. The first entry of each list will always belong to the finest level.

Important things to note:

Not all lists must have the same length: The longest list defines the number of levels and if other lists are shorter, the levels get the last entry in these lists (3 nodes on level 1 and 2 in this example).
As for most other parameters, space_transfer_class and space_transfer_params are part of the description of the problem.
For advanced users: it is also possible to pass parameters to the base_transfer by specifying base_transfer_params or even replace this by defining base_transfer_class in the description.

Full code: pySDC/tutorial/step_4/B_multilevel_hierarchy.py

from pathlib import Path

from pySDC.core.Step import step

from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh


def main():
    """
    A simple test program to set up a full step hierarchy
    """

    # initialize level parameters
    level_params = dict()
    level_params['restol'] = 1e-10
    level_params['dt'] = 0.1

    # initialize sweeper parameters
    sweeper_params = dict()
    sweeper_params['quad_type'] = 'RADAU-RIGHT'
    sweeper_params['num_nodes'] = [5, 3]
    sweeper_params['QI'] = 'LU'

    # initialize problem parameters
    problem_params = dict()
    problem_params['nu'] = 0.1  # diffusion coefficient
    problem_params['freq'] = 4  # frequency for the test value
    problem_params['nvars'] = [31, 15, 7]  # number of degrees of freedom for each level
    problem_params['bc'] = 'dirichlet-zero'  # boundary conditions

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

    # initialize space transfer parameters
    space_transfer_params = dict()
    space_transfer_params['rorder'] = 2
    space_transfer_params['iorder'] = 2

    # fill description dictionary for easy step instantiation
    description = dict()
    description['problem_class'] = heatNd_unforced
    description['problem_params'] = problem_params
    description['sweeper_class'] = generic_implicit
    description['sweeper_params'] = sweeper_params
    description['level_params'] = level_params
    description['step_params'] = step_params
    description['space_transfer_class'] = mesh_to_mesh
    description['space_transfer_params'] = space_transfer_params

    # now the description contains more or less everything we need to create a step with multiple levels
    S = step(description=description)

    # print out and check
    Path("data").mkdir(parents=True, exist_ok=True)
    f = open('data/step_4_B_out.txt', 'w')
    for l in range(len(S.levels)):
        L = S.levels[l]
        out = 'Level %2i: nvars = %4i -- nnodes = %2i' % (l, L.prob.nvars[0], L.sweep.coll.num_nodes)
        f.write(out + '\n')
        print(out)
        assert L.prob.nvars[0] == problem_params['nvars'][min(l, len(problem_params['nvars']) - 1)], (
            "ERROR: number of DOFs is not correct on this level, got %s" % L.prob.nvars
        )
        assert L.sweep.coll.num_nodes == sweeper_params['num_nodes'][min(l, len(sweeper_params['num_nodes']) - 1)], (
            "ERROR: number of nodes is not correct on this level, got %s" % L.sweep.coll.num_nodes
        )
    f.close()


if __name__ == "__main__":
    main()

Results:

Level  0: nvars =   31 -- nnodes =  5
Level  1: nvars =   15 -- nnodes =  3
Level  2: nvars =    7 -- nnodes =  3

Part C: SDC vs. MLSDC¶

After we have seen how a hierarchy is created, we now run MLSDC and compare the results with SDC. We create two different descriptions and controllers and run first SDC and then MLSDC for the simple unforced heat equation. We see that the results are pretty much the same, while MLSDC only takes about half as many iterations.

Important things to note:

In this case, the number of iterations is halved when using MLSDC. This is the best case and in many situations, this cannot be achieved. In particular, the interpolation order is crucial.
Using the controller parameter predict, we can turn the coarse level predictor on (default) or off in the case of MLSDC or PFASST.
While MLSDC looks less expensive, the number of evaluations of the right-hand side of the ODE is basically the same: This is due to the fact that after each coarse grid correction (i.e. after the interpolation), the right-hand side needs to be re-evaluated on the finer level to be prepared for the next sweep. One way of improving this is to do the interpolation also in the right-hand side itself. This can be achieved by specifying base_transfer_params['finter] = True and by passing it to the description. See also Part D.

Full code: pySDC/tutorial/step_4/C_SDC_vs_MLSDC.py

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.HeatEquation_ND_FD import heatNd_unforced
from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh


def main():
    """
    A simple test program to compare SDC and MLSDC
    """

    # initialize level parameters
    level_params = dict()
    level_params['restol'] = 1e-09
    level_params['dt'] = 0.1

    # initialize sweeper parameters
    sweeper_params_sdc = dict()
    sweeper_params_sdc['node_type'] = 'LEGENDRE'
    sweeper_params_sdc['quad_type'] = 'RADAU-RIGHT'
    sweeper_params_sdc['num_nodes'] = 5
    sweeper_params_sdc['QI'] = 'LU'

    sweeper_params_mlsdc = dict()
    sweeper_params_mlsdc['node_type'] = 'LEGENDRE'
    sweeper_params_mlsdc['quad_type'] = 'RADAU-RIGHT'
    sweeper_params_mlsdc['num_nodes'] = [5, 3, 2]
    sweeper_params_mlsdc['QI'] = 'LU'

    # initialize problem parameters
    problem_params_sdc = dict()
    problem_params_sdc['nu'] = 0.1  # diffusion coefficient
    problem_params_sdc['freq'] = 4  # frequency for the test value
    problem_params_sdc['nvars'] = 1023  # number of degrees of freedom for each level
    problem_params_sdc['bc'] = 'dirichlet-zero'  # boundary conditions

    problem_params_mlsdc = dict()
    problem_params_mlsdc['nu'] = 0.1  # diffusion coefficient
    problem_params_mlsdc['freq'] = 4  # frequency for the test value
    problem_params_mlsdc['nvars'] = [1023, 511, 255]  # number of degrees of freedom for each level
    problem_params_mlsdc['bc'] = 'dirichlet-zero'  # boundary conditions

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

    # initialize space transfer parameters
    space_transfer_params = dict()
    space_transfer_params['rorder'] = 2
    space_transfer_params['iorder'] = 6

    # initialize controller parameters
    controller_params = dict()
    controller_params['logger_level'] = 30

    # fill description dictionary for SDC
    description_sdc = dict()
    description_sdc['problem_class'] = heatNd_unforced
    description_sdc['problem_params'] = problem_params_sdc
    description_sdc['sweeper_class'] = generic_implicit
    description_sdc['sweeper_params'] = sweeper_params_sdc
    description_sdc['level_params'] = level_params
    description_sdc['step_params'] = step_params

    # fill description dictionary for MLSDC
    description_mlsdc = dict()
    description_mlsdc['problem_class'] = heatNd_unforced
    description_mlsdc['problem_params'] = problem_params_mlsdc
    description_mlsdc['sweeper_class'] = generic_implicit
    description_mlsdc['sweeper_params'] = sweeper_params_mlsdc
    description_mlsdc['level_params'] = level_params
    description_mlsdc['step_params'] = step_params
    description_mlsdc['space_transfer_class'] = mesh_to_mesh
    description_mlsdc['space_transfer_params'] = space_transfer_params

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

    # set time parameters
    t0 = 0.0
    Tend = 0.1

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

    # call main functions to get things done...
    uend_sdc, stats_sdc = controller_sdc.run(u0=uinit, t0=t0, Tend=Tend)
    uend_mlsdc, stats_mlsdc = controller_mlsdc.run(u0=uinit, t0=t0, Tend=Tend)

    # get number of iterations for both
    niter_sdc = get_sorted(stats_sdc, type='niter', sortby='time')[0][1]
    niter_mlsdc = get_sorted(stats_mlsdc, type='niter', sortby='time')[0][1]

    # compute exact solution and compare both
    uex = P.u_exact(Tend)
    err_sdc = abs(uex - uend_sdc)
    err_mlsdc = abs(uex - uend_mlsdc)
    diff = abs(uend_mlsdc - uend_sdc)

    # print out and check
    Path("data").mkdir(parents=True, exist_ok=True)
    f = open('data/step_4_C_out.txt', 'a')
    out = 'Error SDC and MLSDC: %12.8e -- %12.8e' % (err_sdc, err_mlsdc)
    f.write(out + '\n')
    print(out)
    out = 'Difference SDC vs. MLSDC: %12.8e' % diff
    f.write(out + '\n')
    print(out)
    out = 'Number of iterations SDC and MLSDC: %2i -- %2i' % (niter_sdc, niter_mlsdc)
    f.write(out + '\n')
    print(out)

    assert diff < 6e-10, "ERROR: difference between MLSDC and SDC is higher than expected, got %s" % diff
    assert niter_sdc - niter_mlsdc <= 6, "ERROR: MLSDC required more iterations than expected, got %s" % niter_mlsdc


if __name__ == "__main__":
    main()

Results:

Error SDC and MLSDC: 3.96232037e-08 -- 3.95409337e-08
Difference SDC vs. MLSDC: 8.22700796e-11
Number of iterations SDC and MLSDC: 12 --  6

Part D: MLSDC with particles¶

For this example, we return to the Boris solver and show how coarsening can be done beyond degrees of freedom in space or collocation nodes. Here, we use the setup from step_3, Parts B and C. We run three different versions of SDC or MLSDC: the standard SDC, the standard MLSDC and MSDC with interpolation of the right-hand side. For coarsening, we replace the problem class by a simpler version: the coarse evaluation of the forces omits the particle-particle interaction and only takes external forces into account. This is done simply by replacing the problem class by a list of two problem classes in the description. In the results, we can see that all versions produce more or less the same energies, where MLSDC without f-interpolation takes about half as many iterations and with f-interpolation slightly more. We also check the timings of the three runs: although MLSDC requires much fewer iterations, it takes longer to run. This is due to the fact that the right-hand side of the ODE (i.e. the costly force evaluation) is required after each interpolation! To this end, we also use f-interpolation, which increases the iteration counts a bit but leads to a slightly reduced runtime.

Important things to note:

Again, the number of MLSDC iterations is highly sensitive to the interplay of all the different parameters (number of particles, smoothing parameter, number of nodes, residual tolerance etc.). It is by far not trivial to get a speedup at all (although a reasonable setup has been chosen here).
Of course, this type of coarsening can be combined with the generic coarsening strategies. Yet, using a reduced number of collocation nodes leads to increased iteration counts here.

Full code: pySDC/tutorial/step_4/D_MLSDC_with_particles.py

import time
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.implementations.transfer_classes.TransferParticles_NoCoarse import particles_to_particles
from pySDC.tutorial.step_3.HookClass_Particles import particle_hook
from pySDC.tutorial.step_4.PenningTrap_3D_coarse import penningtrap_coarse


def main():
    """
    A simple test program to compare SDC with two flavors of MLSDC for particle dynamics
    """

    # run SDC, MLSDC and MLSDC plus f-interpolation and compare
    stats_sdc, time_sdc = run_penning_trap_simulation(mlsdc=False)
    stats_mlsdc, time_mlsdc = run_penning_trap_simulation(mlsdc=True)
    stats_mlsdc_finter, time_mlsdc_finter = run_penning_trap_simulation(mlsdc=True, finter=True)

    Path("data").mkdir(parents=True, exist_ok=True)
    f = open('data/step_4_D_out.txt', 'w')
    out = 'Timings for SDC, MLSDC and MLSDC+finter: %12.8f -- %12.8f -- %12.8f' % (
        time_sdc,
        time_mlsdc,
        time_mlsdc_finter,
    )
    f.write(out + '\n')
    print(out)

    # sort and convert stats to list, sorted by iteration numbers (only pre- and after-step are present here)
    energy_sdc = get_sorted(stats_sdc, type='etot', sortby='iter')
    energy_mlsdc = get_sorted(stats_mlsdc, type='etot', sortby='iter')
    energy_mlsdc_finter = get_sorted(stats_mlsdc_finter, type='etot', sortby='iter')

    # get base energy and show differences
    base_energy = energy_sdc[0][1]
    for item in energy_sdc:
        out = 'Total energy and relative deviation in iteration %2i: %12.10f -- %12.8e' % (
            item[0],
            item[1],
            abs(base_energy - item[1]) / base_energy,
        )
        f.write(out + '\n')
        print(out)
    for item in energy_mlsdc:
        out = 'Total energy and relative deviation in iteration %2i: %12.10f -- %12.8e' % (
            item[0],
            item[1],
            abs(base_energy - item[1]) / base_energy,
        )
        f.write(out + '\n')
        print(out)
    for item in energy_mlsdc_finter:
        out = 'Total energy and relative deviation in iteration %2i: %12.10f -- %12.8e' % (
            item[0],
            item[1],
            abs(base_energy - item[1]) / base_energy,
        )
        f.write(out + '\n')
        print(out)
    f.close()

    assert (
        abs(energy_sdc[-1][1] - energy_mlsdc[-1][1]) / base_energy < 6e-10
    ), 'ERROR: energy deviated too much between SDC and MLSDC, got %s' % (
        abs(energy_sdc[-1][1] - energy_mlsdc[-1][1]) / base_energy
    )
    assert (
        abs(energy_mlsdc[-1][1] - energy_mlsdc_finter[-1][1]) / base_energy < 8e-10
    ), 'ERROR: energy deviated too much after using finter, got %s' % (
        abs(energy_mlsdc[-1][1] - energy_mlsdc_finter[-1][1]) / base_energy
    )


def run_penning_trap_simulation(mlsdc, finter=False):
    """
    A simple test program to run IMEX SDC for a single time step
    """
    # initialize level parameters
    level_params = dict()
    level_params['restol'] = 1e-07
    level_params['dt'] = 1.0 / 8

    # initialize sweeper parameters
    sweeper_params = dict()
    sweeper_params['quad_type'] = 'RADAU-RIGHT'
    sweeper_params['num_nodes'] = 5

    # 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'] = 50  # 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['logger_level'] = 30

    transfer_params = dict()
    transfer_params['finter'] = finter

    # Fill description dictionary for easy hierarchy creation
    description = dict()
    if mlsdc:
        # MLSDC: provide list of two problem classes: one for the fine, one for the coarse level
        description['problem_class'] = [penningtrap, penningtrap_coarse]
    else:
        # SDC: provide only one problem class
        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
    description['space_transfer_class'] = particles_to_particles
    description['base_transfer_params'] = transfer_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 and time main function to get things done...
    start_time = time.perf_counter()
    uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
    end_time = time.perf_counter() - start_time

    return stats, end_time


if __name__ == "__main__":
    main()

Results:

Timings for SDC, MLSDC and MLSDC+finter:   3.58500926 --   4.00883495 --   4.62404321
Total energy and relative deviation in iteration  0: 407936.7556966486 -- 0.00000000e+00
Total energy and relative deviation in iteration 12: 406977.9425667246 -- 2.35039652e-03
Total energy and relative deviation in iteration  0: 407936.7556966486 -- 0.00000000e+00
Total energy and relative deviation in iteration  6: 406977.9425660003 -- 2.35039652e-03
Total energy and relative deviation in iteration  0: 407936.7556966486 -- 0.00000000e+00
Total energy and relative deviation in iteration  7: 406977.9428639794 -- 2.35039579e-03

Step-4: Multilevel SDC¶

Part A: Spatial transfer operators¶

Part B: Multilevel hierarchy¶

Part C: SDC vs. MLSDC¶

Part D: MLSDC with particles¶

Table of Contents

Previous topic

Next topic

This Page