Coverage for pySDC/tutorial/step_4/D_MLSDC_with_particles.py: 100%
55 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1# ---
2# jupyter:
3# jupytext:
4# formats: py:percent
5# kernelspec:
6# display_name: Python 3
7# name: python3
8# ---
10# %% [markdown]
11# # Part D: MLSDC with particles
12#
13# Back to the Boris solver and the Penning trap of [Step 3](../step_3/B_adding_statistics), now with 50 particles,
14# to show that coarsening can go beyond the number of unknowns or collocation nodes. The costly part of this
15# problem is the force: every particle interacts with every other one, $O(N^2)$ work. On the coarse level, we use
16# a simpler problem, `penningtrap_coarse` in `PenningTrap_3D_coarse.py` next to this tutorial, which leaves out the
17# particle-particle interaction and only computes the external fields, $O(N)$ work. For this, the problem class in
18# the description simply becomes a list of two classes.
19#
20# The transfer class `particles_to_particles` is a dummy: both levels have the same particles, so there is nothing
21# to restrict or interpolate but the values themselves.
23# %%
24import time
26import matplotlib.pyplot as plt
27import numpy as np
29from pySDC.helpers.stats_helper import get_sorted
30from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
31from pySDC.implementations.problem_classes.PenningTrap_3D import penningtrap
32from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order
33from pySDC.implementations.transfer_classes.TransferParticles_NoCoarse import particles_to_particles
34from pySDC.tutorial.step_3.HookClass_Particles import particle_hook
35from pySDC.tutorial.step_4.PenningTrap_3D_coarse import penningtrap_coarse
38def run_penning_trap_simulation(mlsdc, finter=False):
39 """
40 Run one step of the Penning trap with SDC or MLSDC, and return the statistics, the run time and the number of
41 right-hand side evaluations on each level
42 """
43 # initialize level parameters
44 level_params = {'restol': 1e-07, 'dt': 1.0 / 8}
46 # initialize sweeper parameters
47 sweeper_params = {'quad_type': 'RADAU-RIGHT', 'num_nodes': 5}
49 # initialize problem parameters for the Penning trap
50 problem_params = {
51 'omega_E': 4.9, # E-field frequency
52 'omega_B': 25.0, # B-field frequency
53 'u0': np.array([[10, 0, 0], [100, 0, 100], [1], [1]], dtype=object), # initial center of positions
54 'nparts': 50, # number of particles in the trap
55 'sig': 0.1, # smoothing parameter for the forces
56 }
58 # initialize step parameters
59 step_params = {'maxiter': 20}
61 # initialize controller parameters
62 controller_params = {
63 'hook_class': particle_hook, # specialized hook class for more statistics and output
64 'logger_level': 30,
65 }
67 transfer_params = {'finter': finter}
69 # Fill description dictionary for easy hierarchy creation
70 description = {
71 # MLSDC: a list of two problem classes, one for the fine and one for the coarse level
72 'problem_class': [penningtrap, penningtrap_coarse] if mlsdc else penningtrap,
73 'problem_params': problem_params,
74 'sweeper_class': boris_2nd_order,
75 'sweeper_params': sweeper_params,
76 'level_params': level_params,
77 'step_params': step_params,
78 'space_transfer_class': particles_to_particles,
79 'base_transfer_params': transfer_params,
80 }
82 # instantiate the controller
83 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
85 # set time parameters
86 t0 = 0.0
87 Tend = level_params['dt']
89 # get initial values on finest level
90 P = controller.MS[0].levels[0].prob
91 uinit = P.u_init()
93 # call and time main function to get things done...
94 start_time = time.perf_counter()
95 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
96 end_time = time.perf_counter() - start_time
98 rhs = [L.prob.work_counters['rhs'].niter for L in controller.MS[0].levels]
99 return stats, end_time, rhs
102# %% [markdown]
103# ## Three runs
104#
105# SDC, MLSDC, and MLSDC with interpolation of the right-hand side (`finter`, see
106# [Part C](C_SDC_vs_MLSDC)).
108# %%
109# run SDC, MLSDC and MLSDC plus f-interpolation and compare
110stats_sdc, time_sdc, rhs_sdc = run_penning_trap_simulation(mlsdc=False)
111stats_mlsdc, time_mlsdc, rhs_mlsdc = run_penning_trap_simulation(mlsdc=True)
112stats_mlsdc_finter, time_mlsdc_finter, rhs_mlsdc_finter = run_penning_trap_simulation(mlsdc=True, finter=True)
114runs = {
115 'SDC': (stats_sdc, time_sdc, rhs_sdc),
116 'MLSDC': (stats_mlsdc, time_mlsdc, rhs_mlsdc),
117 'MLSDC+finter': (stats_mlsdc_finter, time_mlsdc_finter, rhs_mlsdc_finter),
118}
119for name, (stats, runtime, rhs) in runs.items():
120 niter = get_sorted(stats, type='niter')[0][1]
121 print(f'{name:13s}: {niter:2d} iterations, {runtime:6.3f} s, right-hand side evaluations per level {rhs}')
123# %% [markdown]
124# ## Energy
125#
126# All three produce more or less the same energy:
128# %%
129# sort and convert stats to list, sorted by iteration numbers (only pre- and after-step are present here)
130energy_sdc = get_sorted(stats_sdc, type='etot', sortby='iter')
131energy_mlsdc = get_sorted(stats_mlsdc, type='etot', sortby='iter')
132energy_mlsdc_finter = get_sorted(stats_mlsdc_finter, type='etot', sortby='iter')
134# get base energy and show differences
135base_energy = energy_sdc[0][1]
136for energy in [energy_sdc, energy_mlsdc, energy_mlsdc_finter]:
137 for item in energy:
138 print(
139 'Total energy and relative deviation in iteration %2i: %12.10f -- %12.8e'
140 % (item[0], item[1], abs(base_energy - item[1]) / base_energy)
141 )
143# %% tags=["hide-input"]
144fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3))
145names = list(runs)
146ax1.bar(names, [get_sorted(runs[n][0], type='niter')[0][1] for n in names], color='#4f6bed')
147ax1.set_title('iterations', fontsize=10)
148ax2.bar(names, [runs[n][2][0] for n in names], color='#e8743b', label='fine: all forces, $O(N^2)$')
149ax2.bar(
150 names,
151 [sum(runs[n][2][1:]) for n in names],
152 bottom=[runs[n][2][0] for n in names],
153 color='#19a979',
154 label='coarse: external fields, $O(N)$',
155)
156ax2.set_title('right-hand side evaluations', fontsize=10)
157ax2.legend(frameon=False, fontsize=8)
158fig.tight_layout()
160# %% [markdown]
161# ## Iterations, work and time
162#
163# MLSDC takes half as many iterations as SDC, and with `finter` one more. But look at the expensive evaluations on
164# the fine level: MLSDC needs as many as SDC, because each interpolation needs the fine right-hand side again, as
165# we saw in Part C. Only interpolating the right-hand side as well (`finter`) saves some of them.
166#
167# Whether that shows in the run times depends on how expensive the fine force really is compared with everything
168# else in a sweep. For 50 particles the forces are cheap, computed for all pairs at once: SDC tends to be fastest,
169# and `finter` costs more than it saves. Press {octicon}`play` **Run in browser** to see for yourself. Timings depend
170# on the machine and the implementation; the counts above do not.
171#
172# :::{admonition} Important things to note
173# - Again, the number of MLSDC iterations is highly sensitive to the interplay of all the parameters (number of
174# particles, smoothing parameter, number of nodes, residual tolerance, ...). It is by far not trivial to get a
175# speedup at all, although a reasonable setup has been chosen here.
176# - Count the work, not only the iterations: what an iteration costs differs between the methods.
177# - This kind of coarsening can be combined with the generic ones of [Part B](B_multilevel_hierarchy). Yet, using
178# fewer collocation nodes on the coarse level increases the iteration counts here.
179# :::
180#
181# The checks the tests run:
183# %%
184deviation_mlsdc = abs(energy_sdc[-1][1] - energy_mlsdc[-1][1]) / base_energy
185deviation_finter = abs(energy_mlsdc[-1][1] - energy_mlsdc_finter[-1][1]) / base_energy
186assert deviation_mlsdc < 6e-10, f'ERROR: energy deviated too much between SDC and MLSDC, got {deviation_mlsdc}'
187assert deviation_finter < 8e-10, f'ERROR: energy deviated too much after using finter, got {deviation_finter}'