Coverage for pySDC/tutorial/step_1/B_spatial_accuracy_check.py: 100%
46 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: Spatial accuracy check
12#
13# One error value says little. A second-order discretisation should make the error four times smaller whenever
14# the mesh width halves, so we repeat the test of [Part A](A_spatial_problem_setup) on a sequence of meshes and
15# measure the order.
16#
17# Along the way, we set up the problem from a parameter dictionary instead of keyword arguments, which is how
18# pySDC is configured from now on.
20# %%
21from collections import namedtuple
23import matplotlib.pyplot as plt
24import numpy as np
26from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
28problem_params = {
29 'nu': 0.1, # diffusion coefficient
30 'freq': 4, # frequency for the test value
31 'bc': 'dirichlet-zero', # boundary conditions
32}
34# create list of nvars to do the accuracy test with
35nvars_list = [2**p - 1 for p in range(4, 15)]
36print(nvars_list)
38# %% [markdown]
39# ## Collecting results
40#
41# We store each error in a dictionary, keyed by an `ID` that says which run it belongs to, and keep the list of
42# `nvars` in the dictionary too. That makes the results self-describing: the functions that analyse them later
43# need nothing else. The pattern pays off once there are several parameters to vary.
45# %%
46# setup id for gathering the results (will sort by nvars)
47ID = namedtuple('ID', 'nvars')
50def run_accuracy_check(nvars_list, problem_params):
51 """
52 Routine to check the error of the Laplacian vs. its FD discretization
54 Args:
55 nvars_list: list of nvars to do the testing with
56 problem_params: dictionary containing the problem-dependent parameters
58 Returns:
59 a dictionary containing the errors and a header (with nvars_list)
60 """
62 results = {}
63 # loop over all nvars
64 for nvars in nvars_list:
65 # setup problem
66 problem_params['nvars'] = nvars
67 prob = heatNd_unforced(**problem_params)
69 # create x values, use only inner points
70 xvalues = np.array([(i + 1) * prob.dx for i in range(prob.nvars[0])])
72 # create a mesh instance and fill it with a sine wave
73 u = prob.u_exact(t=0)
75 # create a mesh instance and fill it with the Laplacian of the sine wave
76 u_lap = prob.dtype_u(init=prob.init)
77 u_lap[:] = -((np.pi * prob.freq[0]) ** 2) * prob.nu * np.sin(np.pi * prob.freq[0] * xvalues)
79 # compare analytic and computed solution using the eval_f routine of the problem class
80 err = abs(prob.eval_f(u, 0) - u_lap)
82 # get id for this nvars and put error into dictionary
83 id = ID(nvars=nvars)
84 results[id] = err
86 # add nvars_list to dictionary for easier access later on
87 results['nvars_list'] = nvars_list
89 return results
92results = run_accuracy_check(nvars_list=nvars_list, problem_params=problem_params)
94# %% [markdown]
95# ## Measuring the order
96#
97# For two consecutive meshes with $N_{i-1}$ and $N_i$ unknowns and errors $e_{i-1}$, $e_i$, the observed order is
98#
99# $$
100# p_i = \frac{\log(e_{i-1} / e_i)}{\log(N_i / N_{i-1})} .
101# $$
104# %%
105def get_accuracy_order(results):
106 """
107 Routine to compute the order of accuracy in space
109 Args:
110 results: the dictionary containing the errors
112 Returns:
113 the list of orders
114 """
116 # retrieve the list of nvars from results
117 assert 'nvars_list' in results, 'ERROR: expecting the list of nvars in the results dictionary'
118 nvars_list = sorted(results['nvars_list'])
120 order = []
121 # loop over two consecutive errors/nvars pairs
122 for i in range(1, len(nvars_list)):
123 # get ids
124 id = ID(nvars=nvars_list[i])
125 id_prev = ID(nvars=nvars_list[i - 1])
127 # compute order as log(prev_error/this_error)/log(this_nvars/old_nvars) <-- depends on the sorting of the list!
128 tmp = np.log(results[id_prev] / results[id]) / np.log(nvars_list[i] / nvars_list[i - 1])
129 order.append(tmp)
131 return order
134order = get_accuracy_order(results)
135for nvars, p in zip(nvars_list[1:], order, strict=True):
136 print(f'nvars = {nvars:5d}: computed order {p:4.3f} (expected 2)')
138# %% [markdown]
139# The numbers hover around 2. On a log-log plot, second order is a line with slope $-2$:
141# %% tags=["hide-input"]
142errors = [results[ID(nvars=n)] for n in nvars_list]
143fig, ax = plt.subplots(figsize=(6, 4))
144ax.loglog(nvars_list, errors, 'o', label='experiment')
145ax.loglog(nvars_list, [errors[0] * (nvars_list[0] / n) ** 2 for n in nvars_list], 'k--', label='2nd order')
146ax.set_xlabel('nvars')
147ax.set_ylabel('abs. error')
148ax.grid(alpha=0.3)
149ax.legend(frameon=False)
150fig.tight_layout()
152# %% [markdown]
153# :::{warning}
154# Test your operators with care: push `nvars` beyond $2^{15}$ and the error grows again. The truncation error is
155# then smaller than the round-off error of the stencil, which divides by $h^2$.
156# :::
157#
158# :::{admonition} Try it yourself
159# :class: tip
160# Replace the finite-difference stencil by a fourth-order one with `problem_params['order'] = 4`. Which order do
161# you measure now, and from which `nvars` on does round-off take over?
162# :::
163#
164# :::{dropdown} Answer
165# Fourth order (4.1 to 4.9 on the coarse meshes, then 3.99), but only up to 2047 unknowns, where the error bottoms
166# out near $10^{-9}$; from 4095 on it grows again. That is much earlier than with the second-order stencil, because
167# the error is so much smaller. The check in the last cell then fails, of course: it expects second order.
168# :::
169#
170# ## Summary
171#
172# - Collect results in a dictionary with IDs for the runs and a header with the metadata.
173# - Measure orders of accuracy, don't eyeball them.
174#
175# The check the tests run:
177# %%
178assert all(np.isclose(order, 2, rtol=0.06)), f"ERROR: spatial order of accuracy is not as expected, got {order}"