Coverage for pySDC/tutorial/step_1/A_spatial_problem_setup.py: 100%

1# --- 

2# jupyter: 

3# jupytext: 

4# formats: py:percent 

5# kernelspec: 

6# display_name: Python 3 

7# name: python3 

8# --- 

9 

10# %% [markdown] 

11# # Part A: Spatial problem setup 

12# 

13# We start as simple as possible: no SDC, no time stepping, just a spatial problem. pySDC calls this a 

14# **problem class**: it holds the spatial discretisation of an equation and knows how to evaluate its right-hand 

15# side. Here we take `heatNd_unforced`, the heat equation 

16# 

17# $$ 

18# u_t = \nu\, \Delta u \quad \text{on } [0, 1], \qquad u(0, t) = u(1, t) = 0, 

19# $$ 

20# 

21# discretised with second-order finite differences. 

22 

23# %% 

24import numpy as np 

25import matplotlib.pyplot as plt 

26 

27from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced 

28 

29prob = heatNd_unforced( 

30 nvars=1023, # number of degrees of freedom 

31 nu=0.1, # diffusion coefficient 

32 freq=4, # frequency of the test function 

33 bc='dirichlet-zero', # boundary conditions 

34) 

35print(f'{prob.nvars[0]} unknowns, mesh width {prob.dx:.2e}, operator of shape {prob.A.shape}') 

36 

37# %% [markdown] 

38# Parameters go in as keyword arguments. Almost everything in pySDC is configured this way; later, the same 

39# arguments are collected in a dictionary and handed over to the controller. 

40# 

41# ## Data types 

42# 

43# The values of a solution live in a pySDC **data type**, here `mesh`, which is a NumPy array with some extra 

44# meta-information. Every problem class says which data type it uses, as `dtype_u` for the solution and `dtype_f` 

45# for the right-hand side, and every data type is created from the problem's `init` attribute: 

46 

47# %% 

48u = prob.dtype_u(init=prob.init) 

49print(type(u).__name__, u.shape, 'init =', prob.init) 

50 

51# %% [markdown] 

52# :::{note} 

53# `init` and `nvars` look the same here, so why two attributes? `init` is whatever the data type needs to set 

54# itself up. For a finite-difference mesh that is just the number of unknowns, but for finite elements it is 

55# the function space. So: whenever you create a data type, use `init`. 

56# ::: 

57# 

58# ## Testing the operator 

59# 

60# Before building anything on top of this discretisation, we check it. We fill `u` with $\sin(\pi k x)$, whose 

61# exact right-hand side is $-\nu (\pi k)^2 \sin(\pi k x)$, and compare that with what `eval_f` computes. We only 

62# use inner points: the boundary values are fixed at zero and are not part of the unknowns. 

63 

64 

65# %% 

66def run_accuracy_check(prob): 

67 """ 

68 Routine to check the error of the Laplacian vs. its FD discretization 

69 

70 Args: 

71 prob: a problem instance 

72 

73 Returns: 

74 the error between the analytic Laplacian and the computed one of a given function 

75 """ 

76 

77 # create x values, use only inner points 

78 xvalues = np.array([(i + 1) * prob.dx for i in range(prob.nvars[0])]) 

79 

80 # create a mesh instance and fill it with a sine wave 

81 u = prob.dtype_u(init=prob.init) 

82 u[:] = np.sin(np.pi * prob.freq[0] * xvalues) 

83 

84 # create a mesh instance and fill it with the Laplacian of the sine wave 

85 u_lap = prob.dtype_u(init=prob.init) 

86 u_lap[:] = -((np.pi * prob.freq[0]) ** 2) * prob.nu * np.sin(np.pi * prob.freq[0] * xvalues) 

87 

88 # compare analytic and computed solution using the eval_f routine of the problem class 

89 err = abs(prob.eval_f(u, 0) - u_lap) 

90 

91 return err 

92 

93 

94err = run_accuracy_check(prob) 

95print(f'Error of the spatial accuracy test: {err:8.6e}') 

96 

97# %% [markdown] 

98# `abs` of a `mesh` is its maximum norm, so this is the largest pointwise error. Relative to the size of the 

99# right-hand side, about $0.1 \cdot (4\pi)^2 \approx 16$, that is small. Seeing it helps: 

100 

101# %% tags=["hide-input"] 

102x = np.array([(i + 1) * prob.dx for i in range(prob.nvars[0])]) 

103u[:] = np.sin(np.pi * prob.freq[0] * x) 

104f_exact = -((np.pi * prob.freq[0]) ** 2) * prob.nu * u 

105f_disc = prob.eval_f(u, 0) 

106 

107fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3)) 

108ax1.plot(x, f_exact, lw=3, alpha=0.4, label='exact') 

109ax1.plot(x, f_disc, '--', label='eval_f') 

110ax1.set_xlabel('$x$') 

111ax1.set_title(r'$\nu\,\Delta u$') 

112ax1.legend(frameon=False) 

113ax2.plot(x, np.abs(f_disc - f_exact)) 

114ax2.set_xlabel('$x$') 

115ax2.set_title('pointwise error') 

116fig.tight_layout() 

117 

118# %% [markdown] 

119# The error has the shape of the solution itself, as expected from the leading term of a second-order 

120# finite-difference stencil, $\frac{h^2}{12} u''''$. In [Part B](B_spatial_accuracy_check) we check that it 

121# really shrinks like $h^2$. 

122# 

123# ## Summary 

124# 

125# - A problem class bundles the discretised equation; parameters are passed as keyword arguments. 

126# - Values live in data types such as `mesh`, created from the problem's `init` attribute. 

127# - Test your operators before you build a time integrator on top of them. 

128# 

129# The tests check exactly what we saw above: 

130 

131# %% 

132assert err <= 2e-04, f"ERROR: the spatial accuracy is higher than expected, got {err}"