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

```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)])
uexact_fine = Pfine.dtype_u(Pfine.init)
uexact_fine[:] = np.sin(np.pi * Pfine.freq * xvalues_fine)

# set exact coarse solution as source
xvalues_coarse = np.array([(i + 1) * Pcoarse.dx for i in range(Pcoarse.nvars)])
uexact_coarse = Pfine.dtype_u(Pcoarse.init)
uexact_coarse[:] = np.sin(np.pi * Pcoarse.freq * 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.

```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['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, L.sweep.coll.num_nodes)
f.write(out + '\n')
print(out)
assert L.prob.nvars == 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['num_nodes'] = 5
sweeper_params_sdc['QI'] = 'LU'

sweeper_params_mlsdc = dict()
sweeper_params_mlsdc['node_type'] = 'LEGENDRE'
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.levels.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')
niter_mlsdc = get_sorted(stats_mlsdc, type='niter', sortby='time')

# 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.96227149e-08 -- 3.95404452e-08
Difference SDC vs. MLSDC: 8.22714952e-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.

```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
for item in energy_sdc:
out = 'Total energy and relative deviation in iteration %2i: %12.10f -- %12.8e' % (
item,
item,
abs(base_energy - item) / 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,
item,
abs(base_energy - item) / 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,
item,
abs(base_energy - item) / base_energy,
)
f.write(out + '\n')
print(out)
f.close()

assert (
abs(energy_sdc[-1] - energy_mlsdc[-1]) / base_energy < 6e-10
), 'ERROR: energy deviated too much between SDC and MLSDC, got %s' % (
abs(energy_sdc[-1] - energy_mlsdc[-1]) / base_energy
)
assert (
abs(energy_mlsdc[-1] - energy_mlsdc_finter[-1]) / base_energy < 8e-10
), 'ERROR: energy deviated too much after using finter, got %s' % (
abs(energy_mlsdc[-1] - energy_mlsdc_finter[-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['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], , ], 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.levels.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:   5.67701618 --   6.41298294 --   7.35475068
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.9425660004 -- 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
```