Coverage for pySDC/tutorial/step_3/C_study_collocations.py: 100%
42 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-09-29 12:50 +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 C: Studying collocation node types
12#
13# In [Part B](B_adding_statistics), the energy of the particle changed quite a bit in a single step. Here we test
14# whether the collocation nodes are to blame, and on the way show how to set up a parameter study with pySDC:
15# describe the whole setup except the parameter to vary, then loop over its values.
17# %%
18import matplotlib.pyplot as plt
19import numpy as np
21from pySDC.helpers.stats_helper import get_sorted
22from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
23from pySDC.implementations.problem_classes.PenningTrap_3D import penningtrap
24from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order
25from pySDC.tutorial.step_3.HookClass_Particles import particle_hook
27# initialize level parameters
28level_params = {'restol': 1e-06, 'dt': 1.0 / 16}
30# initialize sweeper parameters, the node type is set in the loop below
31sweeper_params = {'num_nodes': 3}
33# initialize problem parameters
34problem_params = {
35 'omega_E': 4.9,
36 'omega_B': 25.0,
37 'u0': np.array([[10, 0, 0], [100, 0, 100], [1], [1]], dtype=object),
38 'nparts': 1,
39 'sig': 0.1,
40}
42# initialize step parameters
43step_params = {'maxiter': 20}
45# initialize controller parameters
46controller_params = {
47 'hook_class': particle_hook, # specialized hook class for more statistics and output
48 'logger_level': 30, # reduce verbosity of each run
49}
51# Fill description dictionary for easy hierarchy creation
52description = {
53 'problem_class': penningtrap,
54 'problem_params': problem_params,
55 'sweeper_class': boris_2nd_order,
56 'level_params': level_params,
57 'step_params': step_params,
58}
60# %% [markdown]
61# ## The study
62#
63# For each node type we build a new controller. That is slightly inefficient, but it makes sure that all variables
64# and statistics start afresh. The `stats` of each run go into one dictionary, keyed by the node type.
66# %%
67# assemble and loop over list of collocation classes
68quad_types = ['RADAU-RIGHT', 'GAUSS', 'LOBATTO']
69stats_dict = {}
70for qtype in quad_types:
71 sweeper_params['quad_type'] = qtype
72 description['sweeper_params'] = sweeper_params
74 # instantiate the controller
75 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
77 # set time parameters
78 t0 = 0.0
79 Tend = level_params['dt']
81 # get initial values on finest level
82 P = controller.MS[0].levels[0].prob
83 uinit = P.u_init()
85 # call main function to get things done...
86 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
88 # gather stats in dictionary, collocation classes being the keys
89 stats_dict[qtype] = stats
91# %% [markdown]
92# Then we compare the energy before and after the step, for each node type:
94# %%
95ediff = {}
96for cclass, stats in stats_dict.items():
97 # filter and convert/sort statistics by etot and iterations
98 energy = get_sorted(stats, type='etot', sortby='iter')
99 # compare base and final energy
100 base_energy = energy[0][1]
101 final_energy = energy[-1][1]
102 ediff[cclass] = abs(base_energy - final_energy)
103 print("Energy deviation for %s: %12.8e" % (cclass, ediff[cclass]))
105# %% tags=["hide-input"]
106fig, ax = plt.subplots(figsize=(5, 3))
107ax.bar(list(ediff), list(ediff.values()), color=['#e8743b', '#4f6bed', '#19a979'])
108ax.set_yscale('log')
109ax.set_ylabel('energy deviation after one step')
110ax.axhline(level_params['restol'], color='k', ls='--', label='restol')
111ax.legend(frameon=False)
112fig.tight_layout()
114# %% [markdown]
115# Gauss-Radau loses energy, Gauss-Legendre and Gauss-Lobatto hardly any, about $10^{-5}$. The difference is
116# symmetry. Both Gauss-Legendre and Gauss-Lobatto nodes are symmetric within the step, Gauss-Radau
117# nodes are not.
118#
119# :::{admonition} Try it yourself
120# :class: tip
121# Lower the residual tolerance, e.g. to `level_params['restol'] = 1e-10`, and run the study again. What happens to
122# the energy deviation of each node type?
123# :::
124#
125# :::{dropdown} Answer
126# For the symmetric nodes it drops with the tolerance, to about $6 \cdot 10^{-11}$ (Gauss) and
127# $8 \cdot 10^{-10}$ (Lobatto), at the price of a few more iterations. Rule of thumb: energy is conserved to within
128# an order of magnitude or so of the residual tolerance. Gauss-Radau stays at 14.5 whatever the tolerance: its error is the one of the
129# collocation method itself, not of the iteration. The checks below still pass.
130# :::
131#
132# ## Summary
133#
134# - A parameter study: set up everything but the parameter, then loop, with a new controller for each value.
135# - Working with several `stats` dictionaries is not straightforward, but a meta-dictionary like `stats_dict`
136# helps. Alternatively, process each `stats` right after its run and keep only what you need.
137# - Symmetric collocation nodes conserve the energy of this problem to within an order of magnitude or so of the
138# residual tolerance.
139#
140# The checks the tests run:
142# %%
143# set expected differences and check
144ediff_expect = {'RADAU-RIGHT': 15, 'LOBATTO': 1e-05, 'GAUSS': 3e-05}
145for k, v in ediff.items():
146 assert v < ediff_expect[k], f"ERROR: energy deviated too much, got {ediff[k]}"