Coverage for pySDC/tutorial/step_8/A_adaptive_time_stepping.py: 100%
70 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 A: Adaptive time-stepping
12#
13# So far, all our steps had the same size. For many problems this is wasteful: the solution changes quickly at some
14# times and slowly at others, and a fixed step size has to be small enough for the fastest phase. pySDC can choose
15# the step size on the fly with the `Adaptivity` convergence controller. We try it on the van der Pol oscillator
16#
17# $$
18# u'' - \mu (1 - u^2) u' + u = 0
19# $$
20#
21# with $\mu = 5$, written as a system for $u$ and $u'$ in `vanderpol`. In $[0, 10]$, its solution creeps along, jumps
22# once, and creeps along again. We solve it with implicit SDC and exactly 4 iterations per step. The function below
23# does one run, with the fixed step size `dt` or, if a tolerance `e_tol` is given, adaptively, starting with `dt`.
25# %%
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.convergence_controller_classes.adaptivity import Adaptivity
32from pySDC.implementations.hooks.log_solution import LogSolution
33from pySDC.implementations.hooks.log_work import LogWork
34from pySDC.implementations.problem_classes.Van_der_Pol_implicit import vanderpol
35from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
37mu = 5.0
38Tend = 10.0
41def run_vdp(dt, e_tol=None):
42 description = {
43 'problem_class': vanderpol,
44 'problem_params': {'mu': mu, 'newton_tol': 1e-9, 'relative_tolerance': True},
45 'sweeper_class': generic_implicit,
46 'sweeper_params': {'quad_type': 'RADAU-RIGHT', 'num_nodes': 3, 'QI': 'LU'},
47 'level_params': {'dt': dt}, # no restol: every step does maxiter iterations
48 'step_params': {'maxiter': 4},
49 }
50 if e_tol is not None:
51 description['convergence_controllers'] = {Adaptivity: {'e_tol': e_tol}}
53 controller_params = {
54 'logger_level': 30,
55 'mssdc_jac': False, # Adaptivity refuses the (default) Jacobi mode of multi-step SDC
56 'hook_class': [LogSolution, LogWork],
57 }
58 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
59 P = controller.MS[0].levels[0].prob
60 return controller.run(u0=P.u_exact(0.0), t0=0.0, Tend=Tend)
63# reference solution at Tend, computed once with scipy's solve_ivp at a tight tolerance
64u_ref = vanderpol(mu=mu).u_exact(Tend)
67def summary(uend, stats):
68 """Accepted steps, restarts, the work of all steps, restarted ones included, and the error at the end."""
69 return {
70 't_end': get_sorted(stats, type='u', recomputed=False)[-1][0],
71 'steps': len(get_sorted(stats, type='niter', recomputed=False)),
72 'restarts': sum(me[1] for me in get_sorted(stats, type='restart')),
73 'iterations': sum(me[1] for me in get_sorted(stats, type='niter')),
74 'rhs': sum(me[1] for me in get_sorted(stats, type='work_rhs')),
75 'newton': sum(me[1] for me in get_sorted(stats, type='work_newton')),
76 'error': abs(uend - u_ref),
77 }
80def show(name, me):
81 print(
82 f"{name:>12}: {me['steps']:4d} steps, {me['restarts']:3d} restarts, {me['iterations']:5d} iterations, "
83 f"{me['rhs']:6d} right-hand sides, {me['newton']:5d} Newton iterations, error {me['error']:.2e}"
84 )
87# %% [markdown]
88# `LogSolution` records the solution after each step and `LogWork` the increments of the problem's work counters,
89# here of right-hand side evaluations and Newton iterations. The error is measured against a reference solution
90# from scipy.
91#
92# ## An adaptive run
93#
94# All it takes is the `Adaptivity` convergence controller with a tolerance `e_tol` for the local error. We start with
95# a small step and let it grow:
97# %%
98e_tol = 2e-5
99uend, stats = run_vdp(dt=1e-2, e_tol=e_tol)
100adaptive = summary(uend, stats)
101show('adaptive', adaptive)
103# %% [markdown]
104# `Adaptivity` adds the hooks it needs itself: `LogStepSize` records the step size as `dt`, and the error estimate
105# it relies on comes with `LogEmbeddedErrorEstimate`. The restarts are recorded by `LogRestarts`, which comes with
106# the restarting that `controller_nonMPI` always has. When a step is restarted, its values stay in the statistics;
107# `recomputed=False` in `get_sorted` leaves them out and keeps only the accepted steps.
109# %%
110u = get_sorted(stats, type='u', recomputed=False)
111dt = get_sorted(stats, type='dt', recomputed=False)
112e_em = get_sorted(stats, type='error_embedded_estimate', recomputed=False)
113restart_times = [t for t, restarted in get_sorted(stats, type='restart') if restarted]
115ratio = np.array([e for _, e in e_em]) / e_tol
116print(f'estimate / e_tol of the accepted steps: median {np.median(ratio):.2f}, max {ratio.max():.2f}')
118# %% tags=["hide-input"]
119fig, axs = plt.subplots(3, 1, figsize=(8, 6), sharex=True, constrained_layout=True)
120axs[0].plot([t for t, _ in u], [v[0] for _, v in u], label="$u$")
121axs[0].plot([t for t, _ in u], [v[1] for _, v in u], label="$u'$")
122axs[0].legend(frameon=False)
123axs[0].set_ylabel('solution')
124axs[1].semilogy([t for t, _ in dt], [v for _, v in dt], '.-', color='black')
125for t in restart_times:
126 axs[1].axvline(t, color='grey', ls=':', lw=0.8)
127axs[1].set_ylabel(r'$\Delta t$ (dotted: restarts)')
128axs[2].semilogy([t for t, _ in e_em], [v for _, v in e_em], '.', color='magenta', label='embedded estimate')
129axs[2].axhline(e_tol, color='black', ls='--', label='e_tol')
130axs[2].legend(frameon=False)
131axs[2].set_ylabel('local error')
132axs[2].set_xlabel('time')
133for ax in axs:
134 ax.grid(alpha=0.3)
136# %% [markdown]
137# The step size follows the solution: it grows from the initial $0.01$ to about $1$ within the first time unit,
138# shrinks as the jump approaches, is smallest ($0.015$) at $t \approx 5.3$, where $u$ changes fastest ($u'$ reaches
139# $-7.6$), and grows again to its largest value ($1.28$) after the jump. The estimates of the accepted steps stay below
140# `e_tol`, and all but 5 of the 77 above a tenth of it: their median is $0.48$ `e_tol`.
141#
142# ## How it works
143#
144# - **Error estimate.** Each iteration raises the order of SDC by one, so the solution of the next-to-last iteration
145# is a lower order approximation than that of the last one. `EstimateEmbeddedError` takes the difference of the two
146# at the end of the step as an estimate of the local error of the lower order one: the embedded error estimate.
147# - **New step size.** After the last iteration, the next step size is
148# $\Delta t_\mathrm{new} = \beta \, \Delta t \, (\epsilon_\mathrm{tol} / \epsilon)^{1/k}$, with the estimate
149# $\epsilon$, the tolerance $\epsilon_\mathrm{tol}$ (`e_tol`), the safety factor $\beta$ (`beta`, 0.9 by default),
150# and $k$ the number of iterations, which the code uses as the order of the estimate. Hence, all steps must do the
151# same number of iterations: `Adaptivity` refuses a residual tolerance `restol` of zero or more.
152# - **Restarts.** If $\epsilon$ is not below $\epsilon_\mathrm{tol}$, the step is not accepted but restarted with
153# $\Delta t_\mathrm{new}$, which is then smaller than $\Delta t$. A step that still fails after 10 restarts in a
154# row stops the run with a `ConvergenceError` (parameters `max_restarts` and `crash_after_max_restarts` of
155# `BasicRestarting`).
156#
157# The rule proposes the next step as if the error behaves there as in the current one. Where the error grows from step
158# to step, as on the way to the jump, the proposal is too large and the step is restarted. That is where most of the
159# restarts are: each of them redid one step, and none had to redo a step twice.
160#
161# ## A comparison with fixed steps
162#
163# The same SDC with fixed step sizes, halving them until the error matches that of the adaptive run:
165# %%
166dts_fixed = [1 / 8, 1 / 16, 1 / 32]
167fixed = [summary(*run_vdp(dt=me)) for me in dts_fixed]
168for dt_fixed, me in zip(dts_fixed, fixed, strict=True):
169 show(f'dt = 1/{round(1 / dt_fixed)}', me)
170show('adaptive', adaptive)
172# %% tags=["hide-input"]
173fig, ax = plt.subplots(figsize=(6, 3.5), constrained_layout=True)
174ax.loglog([me['newton'] for me in fixed], [me['error'] for me in fixed], 'o-', label='fixed step size')
175for dt_fixed, me in zip(dts_fixed, fixed, strict=True):
176 ax.annotate(f' 1/{round(1 / dt_fixed)}', (me['newton'], me['error']))
177ax.loglog(adaptive['newton'], adaptive['error'], '*', ms=12, label='adaptive')
178ax.set_xlabel('Newton iterations')
179ax.set_ylabel('error at the end')
180ax.legend(frameon=False)
181ax.grid(alpha=0.3)
183# %% [markdown]
184# With $\Delta t = 1/32$, the error is the same as with adaptivity, $7.3 \cdot 10^{-6}$ against
185# $7.2 \cdot 10^{-6}$. For that, the fixed steps need 320 steps, 1280 iterations, 5120 right-hand side evaluations
186# and 4768 Newton iterations. The adaptive run gets there with 77 steps plus 27 restarts, 416 iterations, 1664
187# right-hand side evaluations and 2072 Newton iterations: about a third of the iterations, although a quarter of them
188# went into restarted steps, and less than half the Newton iterations.
189#
190# :::{admonition} Important things to note
191# - `e_tol` bounds the estimated local error of each step, not the error at the end. Here, a local tolerance of
192# $2 \cdot 10^{-5}$ gave an error of $7.2 \cdot 10^{-6}$ at the end, but the relation depends on the problem.
193# - The step size can be limited with the parameters `dt_min` and `dt_max` of `Adaptivity`, which pass them on to
194# the `StepSizeLimiter` convergence controller.
195# - With `avoid_restarts`, `Adaptivity` estimates from the contraction of the iteration how many more iterations
196# would reach the tolerance, and continues iterating instead of restarting if they are few enough.
197# :::
198#
199# The checks the tests run:
201# %%
202assert all(me['t_end'] == Tend for me in [adaptive] + fixed), 'ERROR: not all runs reached Tend'
203assert ratio.max() < 1, f'ERROR: an accepted step has an estimate above the tolerance: {ratio.max()}'
204assert adaptive['restarts'] > 0, 'ERROR: expected restarts in the adaptive run'
205assert (
206 adaptive['error'] < 1.1 * fixed[-1]['error']
207), f"ERROR: adaptive run is less accurate than dt=1/32: {adaptive['error']:.2e} vs. {fixed[-1]['error']:.2e}"
208assert (
209 adaptive['iterations'] < fixed[-1]['iterations'] / 2
210), f"ERROR: adaptive run needs too many iterations: {adaptive['iterations']} vs. {fixed[-1]['iterations']}"
211assert (
212 adaptive['newton'] < fixed[-1]['newton'] / 2
213), f"ERROR: adaptive run needs too many Newton iterations: {adaptive['newton']} vs. {fixed[-1]['newton']}"