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

## 85 statements

, created at 2024-09-09 14:59 +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

10import scipy.sparse as sp

12from pySDC.core.collocation import CollBase

13from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced

15# setup id for gathering the results (will sort by dt)

16ID = namedtuple('ID', 'dt')

19def main():

20 """

21 A simple test program to compute the order of accuracy in time

22 """

24 # initialize problem parameters

25 problem_params = {}

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

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

28 problem_params['nvars'] = 16383 # number of DOFs in space

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

31 # instantiate problem

32 prob = heatNd_unforced(**problem_params)

34 # instantiate collocation class, relative to the time interval [0,1]

37 # assemble list of dt

38 dt_list = [0.1 / 2**p for p in range(0, 4)]

40 # run accuracy test for all dt

41 results = run_accuracy_check(prob=prob, coll=coll, dt_list=dt_list)

43 # get order of accuracy

44 order = get_accuracy_order(results)

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

47 f = open('data/step_1_D_out.txt', 'w')

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

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

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

51 print(out)

52 f.close()

54 # visualize results

55 plot_accuracy(results)

57 assert os.path.isfile('data/step_1_accuracy_test_coll.png')

59 assert all(np.isclose(order, 2 * coll.num_nodes - 1, rtol=0.4)), (

60 "ERROR: did not get order of accuracy as expected, got %s" % order

61 )

64def run_accuracy_check(prob, coll, dt_list):

65 """

66 Routine to build and solve the linear collocation problem

68 Args:

69 prob: a problem instance

70 coll: a collocation instance

71 dt_list: list of time-step sizes

73 Return:

74 the analytic error of the solved collocation problem

75 """

77 results = {}

78 # loop over all nvars

79 for dt in dt_list:

80 # shrink collocation matrix: first line and column deals with initial value, not needed here

81 Q = coll.Qmat[1:, 1:]

83 # build system matrix M of collocation problem

84 M = sp.eye(prob.nvars[0] * coll.num_nodes) - dt * sp.kron(Q, prob.A)

86 # get initial value at t0 = 0

87 u0 = prob.u_exact(t=0)

88 # fill in u0-vector as right-hand side for the collocation problem

89 u0_coll = np.kron(np.ones(coll.num_nodes), u0)

90 # get exact solution at Tend = dt

91 uend = prob.u_exact(t=dt)

93 # solve collocation problem directly

94 u_coll = sp.linalg.spsolve(M, u0_coll)

96 # compute error

97 err = np.linalg.norm(u_coll[-prob.nvars[0] :] - uend, np.inf)

98 # get id for this dt and store error in results

99 id = ID(dt=dt)

100 results[id] = err

102 # add list of dt to results for easier access

103 results['dt_list'] = dt_list

104 return results

107def get_accuracy_order(results):

108 """

109 Routine to compute the order of accuracy in time

111 Args:

112 results: the dictionary containing the errors

114 Returns:

115 the list of orders

116 """

118 # retrieve the list of dt from results

119 assert 'dt_list' in results, 'ERROR: expecting the list of dt in the results dictionary'

120 dt_list = sorted(results['dt_list'], reverse=True)

122 order = []

123 # loop over two consecutive errors/dt pairs

124 for i in range(1, len(dt_list)):

125 # get ids

126 id = ID(dt=dt_list[i])

127 id_prev = ID(dt=dt_list[i - 1])

129 # compute order as log(prev_error/this_error)/log(this_dt/old_dt) <-- depends on the sorting of the list!

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

131 order.append(tmp)

133 return order

136def plot_accuracy(results):

137 """

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

140 Args:

141 results: the dictionary containing the errors

142 """

144 # retrieve the list of nvars from results

145 assert 'dt_list' in results, 'ERROR: expecting the list of dts in the results dictionary'

146 dt_list = sorted(results['dt_list'])

148 # Set up plotting parameters

149 params = {

150 'legend.fontsize': 20,

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

152 'axes.labelsize': 20,

153 'axes.titlesize': 20,

154 'xtick.labelsize': 16,

155 'ytick.labelsize': 16,

156 'lines.linewidth': 3,

157 }

158 plt.rcParams.update(params)

160 # create new figure

161 plt.figure()

162 # take x-axis limits from nvars_list + some spacning left and right

163 plt.xlim([min(dt_list) / 2, max(dt_list) * 2])

164 plt.xlabel('dt')

165 plt.ylabel('abs. error')

166 plt.grid()

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

169 # get error for first entry in nvars_list

170 id = ID(dt=dt_list[0])

171 base_error = results[id]

172 # assemble optimal errors for 5th order method and plot

173 order_guide_space = [base_error * (2 ** (5 * i)) for i in range(0, len(dt_list))]

174 plt.loglog(dt_list, order_guide_space, color='k', ls='--', label='5th order')

176 min_err = 1e99

177 max_err = 0e00

178 err_list = []

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

180 for dt in dt_list:

181 id = ID(dt=dt)

182 err = results[id]

183 min_err = min(err, min_err)

184 max_err = max(err, max_err)

185 err_list.append(err)

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

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

190 plt.legend(loc=2, ncol=1, numpoints=1)

192 # save plot as PDF, beautify

193 fname = 'data/step_1_accuracy_test_coll.png'

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

196 return None

199if __name__ == "__main__":

200 main()