Coverage for pySDC/tutorial/step_1/D_collocation_accuracy_check.py: 100%
85 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +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]
35 coll = CollBase(num_nodes=3, tleft=0, tright=1, node_type='LEGENDRE', quad_type='RADAU-RIGHT')
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')
188 # adjust y-axis limits, add legend
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()