Coverage for pySDC/tutorial/step_2/A_step_data_structure.py: 100%
23 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: Step data structure
12#
13# Everything in pySDC is built from **steps**. A step represents one time step with whatever hierarchy we
14# prescribe: it holds one or more **levels** (several of them later, for multilevel SDC), and each level holds the
15# problem class, the sweeper, the values at the collocation nodes and its parameters.
16#
17# :::{note}
18# This part is for demonstration only: users normally never touch these data structures themselves. The controller
19# in [Part C](C_using_pySDCs_frontend) does it for them.
20# :::
21#
22# ## What SDC iterates
23#
24# [Step 1](../step_1) ended with the collocation problem of one time step, solved directly. For a general
25# right-hand side $f$, it reads
26#
27# $$
28# \vec u = \vec u_0 + \Delta t\, Q\, \vec f(\vec u) ,
29# $$
30#
31# with $\vec u$ the values at the $M$ nodes. SDC solves it iteratively: it replaces the full matrix $Q$, which
32# couples all nodes, by a lower-triangular approximation $Q_\Delta$, and corrects for the difference with the last
33# iterate,
34#
35# $$
36# \vec u^{k+1} - \Delta t\, Q_\Delta\, \vec f(\vec u^{k+1}) = \vec u_0 + \Delta t\, (Q - Q_\Delta)\, \vec f(\vec u^{k}) .
37# $$
38#
39# Because $Q_\Delta$ is lower triangular, this is solved node by node, from the first to the last: one **sweep**, in
40# which each node needs one implicit solve the size of the spatial problem, like a step of implicit Euler. Once the
41# iteration has converged, $\vec u^{k+1} = \vec u^k$, the $Q_\Delta$ terms cancel, and the result solves the
42# collocation problem, whatever $Q_\Delta$ is. $Q_\Delta$ is the **preconditioner** of the iteration: it decides how
43# fast SDC converges, not what it converges to. The sweeper parameter `QI` chooses it: `'IE'` is implicit Euler
44# from node to node, `'LU'`, used below, takes it from the LU decomposition of $Q^T$ (the "LU trick"). How far an
45# iterate is from the solution is measured by the **residual**
46# $\|\vec u_0 + \Delta t\, Q\, \vec f(\vec u^{k}) - \vec u^{k}\|$. [Part B](B_my_first_sweeper) runs this iteration by hand.
47#
48# ## The description
49#
50# A step is created from a single dictionary, the **description**. It names the classes to use and gives each of
51# them its parameters, again as dictionaries.
53# %%
54from pySDC.core.step import Step
55from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
56from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
58# initialize level parameters
59level_params = {
60 'restol': 1e-10, # stop iterating when the residual is below this
61 'dt': 0.1, # time-step size
62}
64# initialize sweeper parameters
65sweeper_params = {
66 'quad_type': 'RADAU-RIGHT', # collocation node type
67 'num_nodes': 3, # number of collocation nodes
68 'QI': 'LU', # preconditioner of the implicit sweep
69}
71# initialize problem parameters
72problem_params = {
73 'nu': 0.1, # diffusion coefficient
74 'freq': 4, # frequency for the test value
75 'nvars': 1023, # number of degrees of freedom
76 'bc': 'dirichlet-zero', # boundary conditions
77}
79# initialize step parameters
80step_params = {
81 'maxiter': 20, # maximum number of iterations
82}
84# fill description dictionary for easy step instantiation
85description = {
86 'problem_class': heatNd_unforced,
87 'problem_params': problem_params,
88 'sweeper_class': generic_implicit,
89 'sweeper_params': sweeper_params,
90 'level_params': level_params,
91 'step_params': step_params,
92}
94# %% [markdown]
95# Note that the description contains *classes*, not instances: `heatNd_unforced`, not `heatNd_unforced(...)`. pySDC
96# instantiates them itself, with the parameters from the description. Passing classes around makes life much easier,
97# even if it is not the textbook way to program.
98#
99# ## Creating a step
100#
101# That is all a step needs:
103# %%
104# now the description contains more or less everything we need to create a step
105S = Step(description=description)
107# we only have a single level, make a shortcut
108L = S.levels[0]
110# one of the integral parts of each level is the problem class, make a shortcut
111P = L.prob
113print(f'Step with {len(S.levels)} level(s), at most {S.params.maxiter} iterations')
114print(f' └─ level {L.level_index}: dt = {L.params.dt}, restol = {L.params.restol}')
115print(f' ├─ problem: {type(P).__name__} with {P.nvars[0]} unknowns')
116print(f' ├─ sweeper: {type(L.sweep).__name__}, {L.sweep.coll.num_nodes} {L.sweep.params.quad_type} nodes')
117print(f' └─ u, f: values at the initial time and the {len(L.u) - 1} nodes, still empty: {L.u}')
119# %% [markdown]
120# The level holds a list `u` with one entry per collocation node, plus one at the front for the initial value, and
121# the same for the right-hand side `f`. They are empty until someone starts the iteration, which we do in
122# [Part B](B_my_first_sweeper).
123#
124# ## The problem inside the step
125#
126# The problem class inside the level is an ordinary instance, created from `problem_params`. We can run the same
127# operator test on it as in [Step 1, Part A](../step_1/A_spatial_problem_setup):
129# %%
130import numpy as np
132# fill u with a sine wave, and compare eval_f with the exact Laplacian
133x = np.array([(i + 1) * P.dx for i in range(P.nvars[0])])
134u = P.dtype_u(init=P.init)
135u[:] = np.sin(np.pi * P.freq[0] * x)
136err = abs(P.eval_f(u, 0) + (np.pi * P.freq[0]) ** 2 * P.nu * u)
137print(f'Error of the spatial accuracy test: {err:8.6e}')
139# %% [markdown]
140# Same problem, same parameters, same error: the step built exactly what we built by hand before.
141#
142# ## Summary
143#
144# - The `description` dictionary is the central steering tool of pySDC.
145# - A step holds levels; a level holds the problem, the sweeper, the values at the nodes, parameters and status.
146# - Classes, not instances, go into the description.
147#
148# The check the tests run:
150# %%
151assert err <= 2e-04, f"ERROR: the spatial accuracy is higher than expected, got {err}"