Coverage for pySDC/tutorial/step_4/A_spatial_transfer_operators.py: 100%
51 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: Spatial transfer operators
12#
13# Multilevel SDC works on a hierarchy of levels, for example several spatial resolutions of the same problem, and
14# moves values between them: **restriction** from fine to coarse, **interpolation** (or prolongation) from coarse to
15# fine. How this is done depends on the data type, so the user provides it as a `space_transfer_class`. For meshes
16# of finite differences, pySDC has `mesh_to_mesh`, which interpolates with Lagrange polynomials.
17#
18# :::{note}
19# As for the sweeper and everything else, users do not instantiate the transfer class themselves when they use a
20# controller. We do it here to test it.
21# :::
23# %%
24from collections import namedtuple
26import matplotlib.pyplot as plt
27import numpy as np
29from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
30from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
32# initialize problem parameters
33problem_params = {
34 'nu': 0.1, # diffusion coefficient
35 'freq': 3, # frequency for the test value
36 'bc': 'dirichlet-zero', # boundary conditions
37}
39# initialize transfer parameters
40space_transfer_params = {
41 'rorder': 2, # order of the restriction
42 'iorder': 4, # order of the interpolation
43}
45# %% [markdown]
46# ## One interpolation
47#
48# A fine problem and a coarse one with half the unknowns (`nvars` of the form $2^p - 1$ halve nicely, as the
49# boundary points are not unknowns). `prolong` takes a coarse value to the fine mesh:
51# %%
52problem_params['nvars'] = 15
53Pfine = heatNd_unforced(**problem_params)
54problem_params['nvars'] = 7
55Pcoarse = heatNd_unforced(**problem_params)
56T = mesh_to_mesh(fine_prob=Pfine, coarse_prob=Pcoarse, params=space_transfer_params)
58xc = np.array([(i + 1) * Pcoarse.dx for i in range(Pcoarse.nvars[0])])
59xf = np.array([(i + 1) * Pfine.dx for i in range(Pfine.nvars[0])])
60uc = Pcoarse.dtype_u(Pcoarse.init)
61uc[:] = np.sin(np.pi * Pcoarse.freq[0] * xc)
62uf = T.prolong(uc)
64fig, ax = plt.subplots(figsize=(6, 3))
65x = np.linspace(0, 1, 200)
66ax.plot(x, np.sin(np.pi * problem_params['freq'] * x), color='0.7', label='exact')
67ax.plot(xc, uc, 'o', ms=8, label='coarse, 7 unknowns')
68ax.plot(xf, uf, 'x', ms=8, label='interpolated, 15 unknowns')
69ax.set_xlabel('$x$')
70ax.legend(frameon=False, fontsize=8)
71fig.tight_layout()
73# %% [markdown]
74# ## The order of interpolation
75#
76# As before, we check the order: interpolate the exact coarse values for a sequence of meshes and compare with the
77# exact fine values. We expect 4, the `iorder` we asked for.
79# %%
80nvars_fine_list = [2**p - 1 for p in range(5, 10)]
82# setup id for gathering the results (will sort by nvars)
83ID = namedtuple('ID', 'nvars_fine')
85# set up dictionary to store results (plus lists)
86results = {'nvars_list': nvars_fine_list}
88for nvars_fine in nvars_fine_list:
89 print('Working on nvars_fine = %4i...' % nvars_fine)
91 # instantiate fine problem
92 problem_params['nvars'] = nvars_fine # number of degrees of freedom
93 Pfine = heatNd_unforced(**problem_params)
95 # instantiate coarse problem using half of the DOFs
96 problem_params['nvars'] = int((nvars_fine + 1) / 2.0 - 1)
97 Pcoarse = heatNd_unforced(**problem_params)
99 # instantiate spatial interpolation
100 T = mesh_to_mesh(fine_prob=Pfine, coarse_prob=Pcoarse, params=space_transfer_params)
102 # set exact fine solution to compare with
103 xvalues_fine = np.array([(i + 1) * Pfine.dx for i in range(Pfine.nvars[0])])
104 uexact_fine = Pfine.dtype_u(Pfine.init)
105 uexact_fine[:] = np.sin(np.pi * Pfine.freq[0] * xvalues_fine)
107 # set exact coarse solution as source
108 xvalues_coarse = np.array([(i + 1) * Pcoarse.dx for i in range(Pcoarse.nvars[0])])
109 uexact_coarse = Pfine.dtype_u(Pcoarse.init)
110 uexact_coarse[:] = np.sin(np.pi * Pcoarse.freq[0] * xvalues_coarse)
112 # do the interpolation/prolongation
113 uinter = T.prolong(uexact_coarse)
115 # compute error and store
116 id = ID(nvars_fine=nvars_fine)
117 results[id] = abs(uinter - uexact_fine)
119# %%
120# observed order between consecutive meshes, as in step 1, part B
121nvars = sorted(results['nvars_list'])
122errors = [results[ID(nvars_fine=n)] for n in nvars]
123orders = [np.log(errors[i - 1] / errors[i]) / np.log(nvars[i] / nvars[i - 1]) for i in range(1, len(nvars))]
125for p in range(len(orders)):
126 print(
127 'Expected order %2i, got order %5.2f, deviation of %5.2f%%'
128 % (
129 space_transfer_params['iorder'],
130 orders[p],
131 100 * abs(space_transfer_params['iorder'] - orders[p]) / space_transfer_params['iorder'],
132 )
133 )
135# %% [markdown]
136# :::{admonition} Important things to note
137# - MLSDC and PFASST rely on high-order interpolation in space. With Lagrange-based interpolation, orders of 4 to 6
138# and above are recommended.
139# - Restriction, however, can be of order 2, and is thus not tested here.
140# :::
141#
142# :::{admonition} Try it yourself
143# :class: tip
144# Set `space_transfer_params['iorder'] = 6` and run the order check again.
145# :::
146#
147# :::{dropdown} Answer
148# Sixth order, but approached from above: 6.45, 6.83, 6.39 and finally 5.98 on the finest meshes. The coarse ones
149# are not in the asymptotic regime yet, so the check below, which allows 5% everywhere, fails.
150# :::
151#
152# The check the tests run:
154# %%
155for order in orders:
156 assert (
157 abs(space_transfer_params['iorder'] - order) / space_transfer_params['iorder'] < 0.05
158 ), f'ERROR: did not get expected orders for interpolation, got {order}'