Coverage for pySDC/tutorial/step_8/C_iteration_estimator.py: 100%
110 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: Iteration estimator
12#
13# When should SDC, MLSDC or PFASST stop iterating? So far, we stopped when the residual was small enough or after a
14# fixed number of iterations. Another option is to estimate how many iterations it takes to get the error with
15# respect to the exact collocation solution below a tolerance.
16#
17# Two consecutive iterates give an estimate of the contraction factor of the iteration,
18# $\tilde L = \|u^{k} - u^{k-1}\| / \|u^{k-1} - u^{k-2}\|$. If the iteration contracts like this, the error after
19# $K$ iterations is at most $\tilde L^K \|u^1 - u^0\| / (1 - \tilde L)$, which is below a tolerance $\varepsilon$
20# once
21#
22# $$
23# K \geq \frac{\log\big(\varepsilon\, (1 - \tilde L) / \|u^1 - u^0\|\big)}{\log \tilde L} .
24# $$
25#
26# Add a few magic safety constants ($\tilde L$ is capped at 0.9, and $K$ gets 5% more) and you can guess when to
27# stop. In pySDC, this is the convergence controller `CheckIterationEstimatorNonMPI`, with the tolerance as
28# `errtol`:
29#
30# ```python
31# description['convergence_controllers'] = {CheckIterationEstimatorNonMPI: {'errtol': 1e-7}}
32# ```
33#
34# ## Checking it
35#
36# To see whether the estimator does its job, we need the error with respect to the exact collocation solution.
37# The hook `error_output`, in `HookClass_error_output.py` next to this tutorial, gets it with some black magic: before
38# each step, it runs pySDC once more inside the hook, without the estimator, and solves the collocation problem of
39# that step to a residual of $10^{-14}$ or for 50 iterations, whichever comes first. After the step, it records the
40# error with respect to that and to the exact solution of the PDE. For this, the description and the controller
41# parameters are handed to the hook through the step parameters.
42#
43# :::{literalinclude} HookClass_error_output.py
44# :pyobject: error_output.pre_step
45# :::
46#
47# ## Three test cases
48#
49# A forced heat equation with CG as the spatial solver, an advection equation with GMRES, and the nonlinear ODE of
50# Auzinger, each with SDC and with MLSDC. The setups are long, but the new part is only the entry
51# `convergence_controllers` in each:
53# %% tags=["hide-input"]
54import matplotlib.pyplot as plt
55import numpy as np
57from pySDC.helpers.stats_helper import get_sorted
58from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
59from pySDC.implementations.convergence_controller_classes.check_iteration_estimator import CheckIterationEstimatorNonMPI
60from pySDC.implementations.problem_classes.AdvectionEquation_ND_FD import advectionNd
61from pySDC.implementations.problem_classes.Auzinger_implicit import auzinger
62from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_forced
63from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
64from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
65from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
66from pySDC.implementations.transfer_classes.TransferMesh_NoCoarse import mesh_to_mesh as mesh_to_mesh_nc
67from pySDC.tutorial.step_8.HookClass_error_output import error_output
70def setup_diffusion(dt=None, ndim=None, ml=False):
71 # initialize level parameters
72 level_params = {'restol': 1e-10, 'dt': dt, 'nsweeps': 1}
74 # initialize sweeper parameters
75 sweeper_params = {
76 'quad_type': 'RADAU-RIGHT',
77 'num_nodes': 3,
78 'QI': ['LU'], # For the IMEX sweeper, the LU-trick can be activated for the implicit part
79 }
81 # initialize problem parameters
82 problem_params = {
83 'order': 8, # order of accuracy for FD discretization in space
84 'nu': 0.1, # diffusion coefficient
85 'bc': 'periodic', # boundary conditions
86 'freq': tuple(2 for _ in range(ndim)), # frequencies
87 'solver_type': 'CG', # do CG instead of LU
88 'liniter': 10, # number of CG iterations
89 }
90 if ml:
91 problem_params['nvars'] = [tuple(64 for _ in range(ndim)), tuple(32 for _ in range(ndim))] # number of dofs
92 else:
93 problem_params['nvars'] = tuple(64 for _ in range(ndim)) # number of dofs
95 # initialize step parameters
96 step_params = {'maxiter': 50, 'errtol': 1e-07}
98 # initialize space transfer parameters
99 space_transfer_params = {'rorder': 2, 'iorder': 6, 'periodic': True}
101 # setup the iteration estimator
102 convergence_controllers = {CheckIterationEstimatorNonMPI: {'errtol': 1e-7}}
104 # initialize controller parameters
105 controller_params = {'logger_level': 30, 'hook_class': error_output}
107 # fill description dictionary for easy step instantiation
108 description = {
109 'problem_class': heatNd_forced,
110 'problem_params': problem_params,
111 'sweeper_class': imex_1st_order,
112 'sweeper_params': sweeper_params,
113 'level_params': level_params,
114 'step_params': step_params,
115 'convergence_controllers': convergence_controllers,
116 }
117 if ml:
118 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
119 description['space_transfer_params'] = space_transfer_params # pass parameters for spatial transfer
121 return description, controller_params
124def setup_advection(dt=None, ndim=None, ml=False):
125 # initialize level parameters
126 level_params = {'restol': 1e-10, 'dt': dt, 'nsweeps': 1}
128 # initialize sweeper parameters
129 sweeper_params = {'quad_type': 'RADAU-RIGHT', 'num_nodes': 3, 'QI': ['LU']}
131 # initialize problem parameters
132 problem_params = {
133 'order': 6, # order of accuracy for FD discretization in space
134 'stencil_type': 'center', # centered finite differences
135 'bc': 'periodic', # boundary conditions
136 'c': 0.1, # advection speed
137 'freq': tuple(2 for _ in range(ndim)), # frequencies
138 'solver_type': 'GMRES', # do GMRES instead of LU
139 'liniter': 10, # number of GMRES iterations
140 }
141 if ml:
142 problem_params['nvars'] = [tuple(64 for _ in range(ndim)), tuple(32 for _ in range(ndim))] # number of dofs
143 else:
144 problem_params['nvars'] = tuple(64 for _ in range(ndim)) # number of dofs
146 # initialize step parameters
147 step_params = {'maxiter': 50, 'errtol': 1e-07}
149 # initialize space transfer parameters
150 space_transfer_params = {'rorder': 2, 'iorder': 6, 'periodic': True}
152 # setup the iteration estimator
153 convergence_controllers = {CheckIterationEstimatorNonMPI: {'errtol': 1e-7}}
155 # initialize controller parameters
156 controller_params = {'logger_level': 30, 'hook_class': error_output}
158 # fill description dictionary for easy step instantiation
159 description = {
160 'problem_class': advectionNd,
161 'problem_params': problem_params,
162 'sweeper_class': generic_implicit,
163 'sweeper_params': sweeper_params,
164 'level_params': level_params,
165 'step_params': step_params,
166 'convergence_controllers': convergence_controllers,
167 }
168 if ml:
169 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
170 description['space_transfer_params'] = space_transfer_params # pass parameters for spatial transfer
172 return description, controller_params
175def setup_auzinger(dt=None, ml=False):
176 # initialize level parameters
177 level_params = {'restol': 1e-10, 'dt': dt, 'nsweeps': 1}
179 # initialize sweeper parameters
180 sweeper_params = {'quad_type': 'RADAU-RIGHT', 'num_nodes': [3, 2] if ml else 3, 'QI': ['LU']}
182 # initialize problem parameters
183 problem_params = {'newton_tol': 1e-12, 'newton_maxiter': 10}
185 # initialize step parameters
186 step_params = {'maxiter': 50, 'errtol': 1e-07}
188 # setup the iteration estimator
189 convergence_controllers = {CheckIterationEstimatorNonMPI: {'errtol': 1e-7}}
191 # initialize controller parameters
192 controller_params = {'logger_level': 30, 'hook_class': error_output}
194 # fill description dictionary for easy step instantiation
195 description = {
196 'problem_class': auzinger,
197 'problem_params': problem_params,
198 'sweeper_class': generic_implicit,
199 'sweeper_params': sweeper_params,
200 'level_params': level_params,
201 'step_params': step_params,
202 'convergence_controllers': convergence_controllers,
203 }
204 if ml:
205 description['space_transfer_class'] = mesh_to_mesh_nc # pass spatial transfer class
207 return description, controller_params
210# %% [markdown]
211# ## Running them
212#
213# For each run, we print the mean number of iterations, then for each step the error with respect to the PDE and
214# to the collocation solution, and check that the latter is below the tolerance.
217# %%
218def run_simulations(type=None, ndim_list=None, Tend=None, nsteps_list=None, ml=False, nprocs=None):
219 """
220 Run one of the test cases with the iteration estimator, check it and return the collocation errors per step
221 """
223 t0 = None
224 dt = None
225 description = None
226 controller_params = None
228 for ndim in ndim_list:
229 for nsteps in nsteps_list:
230 if type == 'diffusion':
231 # set time parameters
232 t0 = 0.0
233 dt = (Tend - t0) / nsteps
234 description, controller_params = setup_diffusion(dt, ndim, ml)
235 mean_number_of_iterations = 3.00 if ml else 5.75
236 elif type == 'advection':
237 # set time parameters
238 t0 = 0.0
239 dt = (Tend - t0) / nsteps
240 description, controller_params = setup_advection(dt, ndim, ml)
241 mean_number_of_iterations = 2.00 if ml else 4.00
242 elif type == 'auzinger':
243 assert ndim == 1
244 # set time parameters
245 t0 = 0.0
246 dt = (Tend - t0) / nsteps
247 description, controller_params = setup_auzinger(dt, ml)
248 mean_number_of_iterations = 3.62 if ml else 5.62
250 print(f'Running {type} in {ndim} dimensions with time-step size {dt}...\n')
252 # Warning: this is black magic used to run an 'exact' collocation solver for each step within the hooks
253 description['step_params']['description'] = description
254 description['step_params']['controller_params'] = controller_params
256 # instantiate controller
257 controller = controller_nonMPI(
258 num_procs=nprocs, controller_params=controller_params, description=description
259 )
261 # get initial values on finest level
262 P = controller.MS[0].levels[0].prob
263 uinit = P.u_exact(t0)
265 # call main function to get things done...
266 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
268 # filter statistics by type (number of iterations)
269 iter_counts = get_sorted(stats, type='niter', sortby='time')
271 niters = np.array([item[1] for item in iter_counts])
272 print(f' Mean number of iterations: {np.mean(niters):4.2f}')
274 # filter statistics by type (error after time-step)
275 PDE_errors = get_sorted(stats, type='PDE_error_after_step', sortby='time')
276 coll_errors = get_sorted(stats, type='coll_error_after_step', sortby='time')
277 for iters, PDE_err, coll_err in zip(iter_counts, PDE_errors, coll_errors, strict=True):
278 assert coll_err[1] < description['step_params']['errtol'], f'Error too high, got {coll_err[1]:8.4e}'
279 print(
280 f' Errors after step {PDE_err[0]:8.4f} with {iters[1]} iterations: '
281 f'{PDE_err[1]:8.4e} / {coll_err[1]:8.4e}'
282 )
283 print()
285 # filter statistics by type (error after time-step)
286 timing = get_sorted(stats, type='timing_run', sortby='time')
287 print(f'...done, took {timing[0][1]} seconds!')
289 print()
290 print('-----------------------------------------------------------------------------')
292 assert np.isclose(
293 mean_number_of_iterations, np.mean(niters), atol=1e-2
294 ), f'Expected {mean_number_of_iterations:.2f} mean iterations, but got {np.mean(niters):.2f}'
295 return coll_errors
298# %% tags=["scroll-output"]
299coll_errors = {}
300for case in ['diffusion', 'advection', 'auzinger']:
301 for ml in [False, True]:
302 name = f'{case}, {"MLSDC" if ml else "SDC"}'
303 coll_errors[name] = run_simulations(type=case, ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=ml, nprocs=1)
305# %% tags=["hide-input"]
306fig, ax = plt.subplots(figsize=(7, 3.5))
307for name, errors in coll_errors.items():
308 ax.semilogy(*zip(*errors, strict=True), 'o-', label=name)
309ax.axhline(1e-7, color='k', ls='--', label='errtol')
310ax.set_xlabel('time')
311ax.set_ylabel('error w.r.t. collocation solution')
312ax.legend(frameon=False, fontsize=8, ncol=2)
313ax.grid(alpha=0.3)
314fig.tight_layout()
316# %% [markdown]
317# In every run and every step, the estimator stopped the iteration with the error below the tolerance, and without
318# ever knowing the collocation solution. It stops earlier than the residual tolerance of $10^{-10}$ alone would:
319# without the estimator, the mean numbers of iterations are 7.88, 4.62, 5.00, 3.00, 7.50 and 4.50 instead. It is
320# not overly cautious either: the errors end between about $5 \cdot 10^{-11}$ and $3 \cdot 10^{-8}$, below the
321# tolerance of $10^{-7}$ in every step.
322#
323# :::{admonition} Important things to note
324# - The estimator can be used with several parallel steps, too. With the controller parameter `all_to_done`, the
325# estimate of the last step then decides when the whole block stops. It has seen little testing there, though.
326# - It is not available for the parallel `controller_MPI` yet: as its name says, `CheckIterationEstimatorNonMPI`
327# works with the emulated parallelism of `controller_nonMPI` only.
328# :::
329#
330# The checks the tests run are inside `run_simulations`: the collocation error after every step, and the mean
331# number of iterations of each run.