Coverage for pySDC/tutorial/step_4/C_SDC_vs_MLSDC.py: 100%
58 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 C: SDC vs. MLSDC
12#
13# Now we run multilevel SDC and compare it with SDC, for one step of the unforced heat equation. Two descriptions,
14# two controllers: SDC on one level, MLSDC on three, coarsened in space and in the number of collocation nodes.
16# %%
17import matplotlib.pyplot as plt
19from pySDC.helpers.stats_helper import get_sorted
20from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
21from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
22from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
23from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
25# initialize level parameters
26level_params = {'restol': 1e-09, 'dt': 0.1}
28# initialize sweeper parameters
29sweeper_params_sdc = {'node_type': 'LEGENDRE', 'quad_type': 'RADAU-RIGHT', 'num_nodes': 5, 'QI': 'LU'}
30sweeper_params_mlsdc = {'node_type': 'LEGENDRE', 'quad_type': 'RADAU-RIGHT', 'num_nodes': [5, 3, 2], 'QI': 'LU'}
32# initialize problem parameters
33problem_params_sdc = {
34 'nu': 0.1, # diffusion coefficient
35 'freq': 4, # frequency for the test value
36 'nvars': 1023, # number of degrees of freedom
37 'bc': 'dirichlet-zero', # boundary conditions
38}
39problem_params_mlsdc = {
40 'nu': 0.1, # diffusion coefficient
41 'freq': 4, # frequency for the test value
42 'nvars': [1023, 511, 255], # number of degrees of freedom for each level
43 'bc': 'dirichlet-zero', # boundary conditions
44}
46# initialize step parameters
47step_params = {'maxiter': 20}
49# initialize space transfer parameters
50space_transfer_params = {'rorder': 2, 'iorder': 6}
52# initialize controller parameters
53controller_params = {'logger_level': 30}
55# fill description dictionary for SDC
56description_sdc = {
57 'problem_class': heatNd_unforced,
58 'problem_params': problem_params_sdc,
59 'sweeper_class': generic_implicit,
60 'sweeper_params': sweeper_params_sdc,
61 'level_params': level_params,
62 'step_params': step_params,
63}
65# fill description dictionary for MLSDC
66description_mlsdc = {
67 'problem_class': heatNd_unforced,
68 'problem_params': problem_params_mlsdc,
69 'sweeper_class': generic_implicit,
70 'sweeper_params': sweeper_params_mlsdc,
71 'level_params': level_params,
72 'step_params': step_params,
73 'space_transfer_class': mesh_to_mesh,
74 'space_transfer_params': space_transfer_params,
75}
77# instantiate the controllers
78controller_sdc = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description_sdc)
79controller_mlsdc = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description_mlsdc)
81# %% [markdown]
82# To see what the two cost, we also count how often each level evaluates the right-hand side, by wrapping the
83# `eval_f` of each level's problem in a counter:
86# %%
87def count_rhs_evaluations(controller):
88 """Wrap eval_f of each level's problem, and return one counter per level"""
89 counters = []
90 for L in controller.MS[0].levels:
91 counter, eval_f = [0], L.prob.eval_f
93 def counting_eval_f(*args, eval_f=eval_f, counter=counter, **kwargs):
94 counter[0] += 1
95 return eval_f(*args, **kwargs)
97 L.prob.eval_f = counting_eval_f
98 counters.append(counter)
99 return counters
102rhs_sdc = count_rhs_evaluations(controller_sdc)
103rhs_mlsdc = count_rhs_evaluations(controller_mlsdc)
105# %% [markdown]
106# ## Running both
108# %%
109# set time parameters
110t0 = 0.0
111Tend = 0.1
113# get initial values on finest level
114P = controller_sdc.MS[0].levels[0].prob
115uinit = P.u_exact(t0)
117# call main functions to get things done...
118uend_sdc, stats_sdc = controller_sdc.run(u0=uinit, t0=t0, Tend=Tend)
119uend_mlsdc, stats_mlsdc = controller_mlsdc.run(u0=uinit, t0=t0, Tend=Tend)
121# get number of iterations for both
122niter_sdc = get_sorted(stats_sdc, type='niter', sortby='time')[0][1]
123niter_mlsdc = get_sorted(stats_mlsdc, type='niter', sortby='time')[0][1]
125# compute exact solution and compare both
126uex = P.u_exact(Tend)
127err_sdc = abs(uex - uend_sdc)
128err_mlsdc = abs(uex - uend_mlsdc)
129diff = abs(uend_mlsdc - uend_sdc)
131print('Error SDC and MLSDC: %12.8e -- %12.8e' % (err_sdc, err_mlsdc))
132print('Difference SDC vs. MLSDC: %12.8e' % diff)
133print('Number of iterations SDC and MLSDC: %2i -- %2i' % (niter_sdc, niter_mlsdc))
134print('Right-hand side evaluations per level, SDC and MLSDC:', [c[0] for c in rhs_sdc], [c[0] for c in rhs_mlsdc])
136# %% tags=["hide-input"]
137fig, ax = plt.subplots(figsize=(6, 3.5))
138for name, stats in [('SDC', stats_sdc), ('MLSDC', stats_mlsdc)]:
139 residual = get_sorted(stats, type='residual_post_iteration', sortby='iter')
140 ax.semilogy(*zip(*residual, strict=True), 'o-', label=name)
141ax.axhline(level_params['restol'], color='k', ls='--', label='restol')
142ax.set_xlabel('iteration')
143ax.set_ylabel('residual')
144ax.legend(frameon=False)
145ax.grid(alpha=0.3)
146fig.tight_layout()
148# %% [markdown]
149# The same result, in half the iterations. This is the best case, and in many situations it cannot be achieved;
150# in particular, the interpolation order is crucial.
151#
152# The count tells the other half of the story: MLSDC evaluates the right-hand side on the finest level exactly as
153# often as SDC, plus the evaluations on the coarse levels. After each coarse-level correction, the interpolated
154# values need their right-hand side on the fine level again, before the next fine sweep. What halves with the
155# iterations are the fine sweeps, and with them the implicit solves on the finest level; the fine right-hand side
156# work does not drop at all.
157#
158# :::{admonition} Try it yourself
159# :class: tip
160# 1. Interpolate the right-hand side too, instead of evaluating it again: add
161# `description_mlsdc['base_transfer_params'] = {'finter': True}` before the controllers are created, and run
162# everything again.
163# 2. Switch on a predictor with the controller parameter `predict_type`: the default `None` has none,
164# `'pfasst_burnin'` sweeps on the coarsest level first. How many iterations does MLSDC need then? (The SDC
165# controller warns that it ignores the predictor: it has only one level.)
166# :::
167#
168# :::{dropdown} Answers
169# 1. The iterations and the result stay the same, but the finest level needs 36 evaluations instead of 66.
170# 2. 4 instead of 6. The result stays within the tolerance of SDC's: the difference grows from $8 \cdot 10^{-11}$
171# to $2 \cdot 10^{-10}$. The predictor's own sweeps do not count as iterations, though.
172# :::
173#
174# The checks the tests run:
176# %%
177assert diff < 6e-10, f"ERROR: difference between MLSDC and SDC is higher than expected, got {diff}"
178assert niter_sdc - niter_mlsdc >= 6, f"ERROR: MLSDC required more iterations than expected, got {niter_mlsdc}"