Coverage for pySDC/tutorial/step_2/B_my_first_sweeper.py: 100%
42 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 B: My first sweeper
12#
13# Now we run SDC for one time step, by hand. The problem is the heat equation with a forcing term,
14# `heatNd_forced`, and the sweeper `imex_1st_order` treats the two parts differently: diffusion implicitly,
15# forcing explicitly. This is **IMEX SDC**, and it needs the right-hand side in two parts, which the data type
16# `imex_mesh` provides as `.impl` and `.expl`.
17#
18# :::{note}
19# Again, this is for demonstration only. Users normally let a controller do all of this, see
20# [Part C](C_using_pySDCs_frontend).
21# :::
23# %%
24import matplotlib.pyplot as plt
26from pySDC.core.step import Step
27from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_forced
28from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
30# initialize level parameters
31level_params = {'restol': 1e-10, 'dt': 0.1}
33# initialize sweeper parameters
34sweeper_params = {'quad_type': 'RADAU-RIGHT', 'num_nodes': 3}
36# initialize problem parameters
37problem_params = {
38 'nu': 0.1, # diffusion coefficient
39 'freq': 4, # frequency for the test value
40 'nvars': 1023, # number of degrees of freedom
41 'bc': 'dirichlet-zero', # boundary conditions
42}
44# initialize step parameters
45step_params = {'maxiter': 20}
47# Fill description dictionary for easy hierarchy creation
48description = {
49 'problem_class': heatNd_forced,
50 'problem_params': problem_params,
51 'sweeper_class': imex_1st_order,
52 'sweeper_params': sweeper_params,
53 'level_params': level_params,
54 'step_params': step_params,
55}
57# instantiate the step we are going to work on, and make shortcuts for the level and the problem
58S = Step(description=description)
59L = S.levels[0]
60P = L.prob
62# %% [markdown]
63# ## Parameters and status
64#
65# Steps and levels carry two kinds of attributes. **Parameters** are what we asked for in the description, such as
66# `S.params.maxiter` or `L.params.restol`. The **status** reflects where the computation currently is, such as
67# `S.status.iter`, `L.status.time` or `L.status.residual`. We set the status of the level to the start of the step
68# and put the initial value at index 0:
70# %%
71# set initial time in the status of the level
72L.status.time = 0.1
73# compute initial value (using the exact function here)
74L.u[0] = P.u_exact(L.time)
76# %% [markdown]
77# ## Predict, sweep, check
78#
79# The sweeper does the actual work, with three routines. `predict` fills the nodes with a first guess; without it,
80# they stay empty. `compute_residual` measures how far the values at the nodes are from solving the collocation
81# problem. `update_nodes` is one SDC sweep.
83# %%
84# access the sweeper's predict routine to get things started
85L.sweep.predict()
86# compute the residual (we may be done already!)
87L.sweep.compute_residual()
89print(f'right-hand side at node 1: {type(L.f[1]).__name__} with parts .impl and .expl')
90print(f'residual of the initial guess: {L.status.residual:12.8e}')
92# %% [markdown]
93# Now we sweep until the residual is below `restol` or we run out of iterations:
95# %%
96# reset iteration counter
97S.status.iter = 0
98residuals = [L.status.residual]
99# run the SDC iteration until either the maximum number of iterations is reached or the residual is small enough
100while S.status.iter < S.params.maxiter and L.status.residual > L.params.restol:
101 # this is where the nodes are actually updated according to the SDC formulas
102 L.sweep.update_nodes()
103 # compute/update the residual
104 L.sweep.compute_residual()
105 # increment the iteration counter
106 S.status.iter += 1
107 residuals.append(L.status.residual)
108 print(
109 f'Time {L.time:4.2f} of {L.level_index} -- Iteration: {S.status.iter:2d} -- Residual: {L.status.residual:12.8e}'
110 )
112# %% tags=["hide-input"]
113fig, ax = plt.subplots(figsize=(6, 3.5))
114ax.semilogy(residuals, 'o-')
115ax.axhline(L.params.restol, color='k', ls='--', label='restol')
116ax.set_xlabel('iteration')
117ax.set_ylabel('residual')
118ax.grid(alpha=0.3)
119ax.legend(frameon=False)
120fig.tight_layout()
122# %% [markdown]
123# Most sweeps reduce the residual by a factor of five to ten, with a slower patch around the sixth. Once converged,
124# the values at the nodes solve a collocation problem like the one in
125# [Step 1, Part C](../step_1/C_collocation_problem_setup), without ever assembling it.
126#
127# ## Finishing the step
128#
129# The solution at the end of the step is computed from the values at the nodes by `compute_end_point`, and only
130# there. Then the time moves on and we compare with the exact solution:
132# %%
133# compute the interval's endpoint: this (and only this) will set uend, depending on the collocation nodes
134L.sweep.compute_end_point()
135# update the simulation time
136L.status.time += L.dt
138# compute exact solution and compare
139err = abs(P.u_exact(L.status.time) - L.uend)
140res, niter = L.status.residual, S.status.iter
141print(f'Error and residual: {err:12.8e} -- {res:12.8e}')
143# %% [markdown]
144# :::{admonition} Try it yourself
145# :class: tip
146# The sweeper preconditions each sweep with a lower-triangular matrix $Q_\Delta$, chosen with the sweeper parameter
147# `QI`. The default is `'IE'` (implicit Euler). Set `sweeper_params['QI'] = 'LU'` in the first cell, or
148# `'MIN-SR-S'`, and run everything again. How many iterations do you need now?
149# :::
150#
151# :::{dropdown} Answer
152# `'LU'` needs 10 iterations instead of 12, `'MIN-SR-S'` needs 11. The error stays the same: all of them converge
153# to the same collocation solution, only the path differs. `'MIN-SR-S'` is diagonal, so its nodes could be solved in
154# parallel, which is what the project [parallel SDC](../../projects/parallelSDC) is about.
155# :::
156#
157# ## Summary
158#
159# - One SDC step is: set the initial value, `predict`, then `update_nodes` and `compute_residual` until converged,
160# then `compute_end_point`.
161# - Parameters are what you asked for; status is where the computation is.
162# - This logic is simple but tedious, and it gets much worse with several steps, levels or processes. Hence
163# [Part C](C_using_pySDCs_frontend).
164#
165# The checks the tests run:
167# %%
168assert err <= 1e-5, f"ERROR: IMEX SDC iteration did not reduce the error enough, got {err}"
169assert res <= level_params['restol'], f"ERROR: IMEX SDC iteration did not reduce the residual enough, got {res}"
170assert niter <= 12, f"ERROR: IMEX SDC took too many iterations, got {niter}"