# Coverage for pySDC/tutorial/step_1/B_spatial_accuracy_check.py: 100%

## 81 statements

, created at 2024-09-20 17:10 +0000

1from pathlib import Path

2import matplotlib

4matplotlib.use('Agg')

6from collections import namedtuple

7import matplotlib.pylab as plt

8import numpy as np

9import os.path

11from pySDC.implementations.datatype_classes.mesh import mesh

12from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced

14# setup id for gathering the results (will sort by nvars)

15ID = namedtuple('ID', 'nvars')

18def main():

19 """

20 A simple test program to check order of accuracy in space for a simple test problem

21 """

23 # initialize problem parameters

24 problem_params = dict()

25 problem_params['nu'] = 0.1 # diffusion coefficient

26 problem_params['freq'] = 4 # frequency for the test value

27 problem_params['bc'] = 'dirichlet-zero' # boundary conditions

29 # create list of nvars to do the accuracy test with

30 nvars_list = [2**p - 1 for p in range(4, 15)]

32 # run accuracy test for all nvars

33 results = run_accuracy_check(nvars_list=nvars_list, problem_params=problem_params)

35 # compute order of accuracy

36 order = get_accuracy_order(results)

38 Path("data").mkdir(parents=True, exist_ok=True)

39 f = open('data/step_1_B_out.txt', 'w')

40 for l in range(len(order)):

41 out = 'Expected order: %2i -- Computed order %4.3f' % (2, order[l])

42 f.write(out + '\n')

43 print(out)

44 f.close()

46 # visualize results

47 plot_accuracy(results)

49 assert os.path.isfile('data/step_1_accuracy_test_space.png'), 'ERROR: plotting did not create file'

51 assert all(np.isclose(order, 2, rtol=0.06)), "ERROR: spatial order of accuracy is not as expected, got %s" % order

54def run_accuracy_check(nvars_list, problem_params):

55 """

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

58 Args:

59 nvars_list: list of nvars to do the testing with

60 problem_params: dictionary containing the problem-dependent parameters

62 Returns:

63 a dictionary containing the errors and a header (with nvars_list)

64 """

66 results = {}

67 # loop over all nvars

68 for nvars in nvars_list:

69 # setup problem

70 problem_params['nvars'] = nvars

71 prob = heatNd_unforced(**problem_params)

73 # create x values, use only inner points

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

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

77 u = prob.u_exact(t=0)

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

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

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

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

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

86 # get id for this nvars and put error into dictionary

87 id = ID(nvars=nvars)

88 results[id] = err

90 # add nvars_list to dictionary for easier access later on

91 results['nvars_list'] = nvars_list

93 return results

96def get_accuracy_order(results):

97 """

98 Routine to compute the order of accuracy in space

100 Args:

101 results: the dictionary containing the errors

103 Returns:

104 the list of orders

105 """

107 # retrieve the list of nvars from results

108 assert 'nvars_list' in results, 'ERROR: expecting the list of nvars in the results dictionary'

109 nvars_list = sorted(results['nvars_list'])

111 order = []

112 # loop over two consecutive errors/nvars pairs

113 for i in range(1, len(nvars_list)):

114 # get ids

115 id = ID(nvars=nvars_list[i])

116 id_prev = ID(nvars=nvars_list[i - 1])

118 # compute order as log(prev_error/this_error)/log(this_nvars/old_nvars) <-- depends on the sorting of the list!

119 tmp = np.log(results[id_prev] / results[id]) / np.log(nvars_list[i] / nvars_list[i - 1])

120 order.append(tmp)

122 return order

125def plot_accuracy(results):

126 """

127 Routine to visualize the errors as well as the expected errors

129 Args:

130 results: the dictionary containing the errors

131 """

133 # retrieve the list of nvars from results

134 assert 'nvars_list' in results, 'ERROR: expecting the list of nvars in the results dictionary'

135 nvars_list = sorted(results['nvars_list'])

137 # Set up plotting parameters

138 params = {

139 'legend.fontsize': 20,

140 'figure.figsize': (12, 8),

141 'axes.labelsize': 20,

142 'axes.titlesize': 20,

143 'xtick.labelsize': 16,

144 'ytick.labelsize': 16,

145 'lines.linewidth': 3,

146 }

147 plt.rcParams.update(params)

149 # create new figure

150 plt.figure()

151 # take x-axis limits from nvars_list + some spacing left and right

152 plt.xlim([min(nvars_list) / 2, max(nvars_list) * 2])

153 plt.xlabel('nvars')

154 plt.ylabel('abs. error')

155 plt.grid()

157 # get guide for the order of accuracy, i.e. the errors to expect

158 # get error for first entry in nvars_list

159 id = ID(nvars=nvars_list[0])

160 base_error = results[id]

161 # assemble optimal errors for 2nd order method and plot

162 order_guide_space = [base_error / (2 ** (2 * i)) for i in range(0, len(nvars_list))]

163 plt.loglog(nvars_list, order_guide_space, color='k', ls='--', label='2nd order')

165 min_err = 1e99

166 max_err = 0e00

167 err_list = []

168 # loop over nvars, get errors and find min/max error for y-axis limits

169 for nvars in nvars_list:

170 id = ID(nvars=nvars)

171 err = results[id]

172 min_err = min(err, min_err)

173 max_err = max(err, max_err)

174 err_list.append(err)

175 plt.loglog(nvars_list, err_list, ls=' ', marker='o', markersize=10, label='experiment')

178 plt.ylim([min_err / 10, max_err * 10])

179 plt.legend(loc=1, ncol=1, numpoints=1)

181 # save plot as PDF, beautify

182 fname = 'data/step_1_accuracy_test_space.png'

183 plt.savefig(fname, bbox_inches='tight')

185 return None

188if __name__ == "__main__":

189 main()