Coverage for pySDC/tutorial/step_1/B_spatial_accuracy_check.py: 100%
81 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
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')
177 # adjust y-axis limits, add legend
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()