Coverage for pySDC/projects/Second_orderSDC/plot_helper.py: 100%
4 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
1# import matplotlib
3# matplotlib.use('Agg')
4# import os
6import numpy as np
7import matplotlib.pyplot as plt
9FONT_SIZE = 16
10FIG_SIZE = (7.44, 6.74)
13def set_fixed_plot_params(): # pragma: no cover
14 """
15 Set fixed parameters for all plots
16 """
17 plt.rcParams['figure.figsize'] = FIG_SIZE
18 plt.rcParams['pgf.rcfonts'] = False
20 plt.rcParams['lines.linewidth'] = 2.5
21 plt.rcParams['axes.titlesize'] = FONT_SIZE + 5
22 plt.rcParams['axes.labelsize'] = FONT_SIZE + 5
23 plt.rcParams['xtick.labelsize'] = FONT_SIZE
24 plt.rcParams['ytick.labelsize'] = FONT_SIZE
25 plt.rcParams['xtick.major.pad'] = 5
26 plt.rcParams['ytick.major.pad'] = 5
27 plt.rcParams['axes.labelpad'] = 6
28 plt.rcParams['lines.markersize'] = FONT_SIZE - 2
29 plt.rcParams['lines.markeredgewidth'] = 1
30 plt.rcParams['mathtext.fontset'] = 'cm'
31 plt.rcParams['mathtext.rm'] = 'serif'
32 plt.rc('font', size=FONT_SIZE)
35class PlotManager(object): # pragma: no cover
36 """
37 This class generates all of the plots of the Second-order SDC plots.
38 """
40 def __init__(self, controller_params, description, time_iter=3, K_iter=(1, 2, 3), Tend=2, axes=(1,), cwd=''):
41 self.controller_params = controller_params
42 self.description = description
43 self.time_iter = time_iter
44 self.K_iter = K_iter
45 self.Tend = Tend
46 self.axes = axes
47 self.cwd = cwd
48 self.quad_type = self.description['sweeper_params']['quad_type']
49 self.num_nodes = self.description['sweeper_params']['num_nodes']
50 self.error_type = 'local'
52 def plot_convergence(self):
53 """
54 Plot convergence order plots for the position and velocity
55 If you change parameters of the values you need set y_lim values need to set manually
56 """
57 set_fixed_plot_params()
58 [N, time_data, error_data, order_data, convline] = self.organize_data(
59 filename='data/dt_vs_{}_errorSDC.csv'.format(self.error_type)
60 )
62 color = ['r', 'brown', 'g', 'blue']
63 shape = ['o', 'd', 's', 'x']
65 fig1, ax1 = plt.subplots()
66 fig2, ax2 = plt.subplots()
67 value = self.axes[0]
68 for ii in range(0, N):
69 ax1.loglog(time_data[ii, :], convline['pos'][value, ii, :], color='black')
70 ax1.loglog(
71 time_data[ii, :],
72 error_data['pos'][value, ii, :],
73 ' ',
74 color=color[ii],
75 marker=shape[ii],
76 label='k={}'.format(int(self.K_iter[ii])),
77 )
78 if value == 2:
79 ax1.text(
80 time_data[ii, 1],
81 0.25 * convline['pos'][value, ii, 1],
82 r"$\mathcal{O}(\Delta t^{%d})$" % (order_data['pos'][ii, 0, 1]),
83 size=18,
84 )
85 else:
86 ax1.text(
87 time_data[ii, 1],
88 0.25 * convline['pos'][value, ii, 1],
89 r"$\mathcal{O}(\Delta t^{%d})$" % (order_data['pos'][ii, 0, 0]),
90 size=18,
91 )
93 if self.error_type == 'Local':
94 ax1.set_ylabel(r'$\Delta x^{\mathrm{(abs)}}_{%d}$' % (value + 1))
95 else:
96 ax1.set_ylabel(r'$\Delta x^{\mathrm{(rel)}}_{%d}$' % (value + 1))
97 ax1.set_title('{} order of convergence, $M={}$'.format(self.error_type, self.num_nodes))
98 ax1.set_xlabel(r'$\omega_{B} \cdot \Delta t$')
100 ax1.legend(loc='best')
101 fig1.tight_layout()
102 fig1.savefig(self.cwd + 'data/{}_conv_plot_pos{}.pdf'.format(self.error_type, value + 1))
104 for ii in range(0, N):
105 ax2.loglog(time_data[ii, :], convline['vel'][value, ii, :], color='black')
106 ax2.loglog(
107 time_data[ii, :],
108 error_data['vel'][value, ii, :],
109 ' ',
110 color=color[ii],
111 marker=shape[ii],
112 label='k={}'.format(int(self.K_iter[ii])),
113 )
115 if value == 2:
116 ax2.text(
117 time_data[ii, 1],
118 0.25 * convline['vel'][value, ii, 1],
119 r"$\mathcal{O}(\Delta t^{%d})$" % (order_data['vel'][ii, 0, 1]),
120 size=18,
121 )
122 else:
123 ax2.text(
124 time_data[ii, 1],
125 0.25 * convline['vel'][value, ii, 1],
126 r"$\mathcal{O}(\Delta t^{%d})$" % (order_data['vel'][ii, 0, 0]),
127 size=18,
128 )
130 if self.error_type == 'Local':
131 ax2.set_ylabel(r'$\Delta v^{\mathrm{(abs)}}_{%d}$' % (value + 1))
132 else:
133 ax2.set_ylabel(r'$\Delta v^{\mathrm{(rel)}}_{%d}$' % (value + 1))
134 ax2.set_title(r'{} order of convergence, $M={}$'.format(self.error_type, self.num_nodes))
135 ax2.set_xlabel(r'$\omega_{B} \cdot \Delta t$')
136 # =============================================================================
137 # Setting y axis min and max values
138 # =============================================================================
139 if self.error_type == 'global':
140 ax2.set_ylim(1e-14, 1e1)
141 ax1.set_ylim(1e-14, 1e1)
142 else:
143 ax2.set_ylim(np.min(ax1.get_ylim()), np.max(ax2.get_ylim()))
144 ax1.set_ylim(np.min(ax1.get_ylim()), np.max(ax2.get_ylim()))
145 ax2.legend(loc='best')
146 fig2.tight_layout()
147 plt.show()
148 fig2.savefig(self.cwd + 'data/{}_conv_plot_vel{}.pdf'.format(self.error_type, value + 1))
150 def format_number(self, data_value, indx):
151 """
152 Change format of the x axis for the work precision plots
153 """
154 if data_value >= 1_000_000:
155 formatter = "{:1.1f}M".format(data_value * 0.000_001)
156 else:
157 formatter = "{:1.0f}K".format(data_value * 0.001)
158 return formatter
160 def plot_work_precision(self):
161 """
162 Generate work precision plots
163 """
164 set_fixed_plot_params()
165 [N, func_eval_SDC, error_SDC, *_] = self.organize_data(
166 filename=self.cwd + 'data/rhs_eval_vs_global_errorSDC.csv',
167 time_iter=self.time_iter,
168 )
170 [N, func_eval_picard, error_picard, *_] = self.organize_data(
171 filename=self.cwd + 'data/rhs_eval_vs_global_errorPicard.csv',
172 time_iter=self.time_iter,
173 )
175 color = ['r', 'brown', 'g', 'blue']
176 shape = ['o', 'd', 's', 'x']
177 fig1, ax1 = plt.subplots()
178 fig2, ax2 = plt.subplots()
179 value = self.axes[0]
181 if self.RK:
182 [N, func_eval_RKN, error_RKN, *_] = self.organize_data(
183 filename=self.cwd + 'data/rhs_eval_vs_global_errorRKN.csv',
184 time_iter=self.time_iter,
185 )
187 ax1.loglog(
188 func_eval_RKN[0],
189 error_RKN['pos'][value,][0][:],
190 ls='dashdot',
191 color='purple',
192 marker='p',
193 label='RKN-4',
194 )
195 ax2.loglog(
196 func_eval_RKN[0],
197 error_RKN['vel'][value,][0][:],
198 ls='dashdot',
199 color='purple',
200 marker='p',
201 label='RKN-4',
202 )
203 if self.VV:
204 [N, func_eval_VV, error_VV, *_] = self.organize_data(
205 filename=self.cwd + 'data/rhs_eval_vs_global_errorVV.csv',
206 time_iter=self.time_iter,
207 )
209 ax1.loglog(
210 func_eval_VV[0],
211 error_VV['pos'][value,][0][:],
212 ls='dashdot',
213 color='blue',
214 marker='H',
215 label='Velocity-Verlet',
216 )
217 ax2.loglog(
218 func_eval_VV[0],
219 error_VV['vel'][value,][0][:],
220 ls='dashdot',
221 color='blue',
222 marker='H',
223 label='Velocity-Verlet',
224 )
226 for ii, jj in enumerate(self.K_iter):
227 # =============================================================================
228 # # If you want to get exactly the same picture like in paper uncomment this only for vertical axis
229 # if ii==0 or ii==1:
230 # ax1.loglog(func_eval_SDC[ii, :][1:], error_SDC['pos'][value, ii, :][1:], ls='solid', color=color[ii], marker=shape[ii], label='k={}'.format(jj))
231 # ax1.loglog(func_eval_picard[ii,:][1:], error_picard['pos'][value, ii, :][1:], ls='--', color=color[ii], marker=shape[ii])
233 # ax2.loglog(func_eval_SDC[ii, :][1:], error_SDC['vel'][value, ii, :][1:], ls='solid', color=color[ii], marker=shape[ii], label='k={}'.format(jj))
234 # ax2.loglog(func_eval_picard[ii,:][1:], error_picard['vel'][value, ii, :][1:], ls='--', color=color[ii], marker=shape[ii])
235 # else:
237 # ax1.loglog(func_eval_SDC[ii, :][:-1], error_SDC['pos'][value, ii, :][:-1], ls='solid', color=color[ii], marker=shape[ii], label='k={}'.format(jj))
238 # ax1.loglog(func_eval_picard[ii,:][:-1], error_picard['pos'][value, ii, :][:-1], ls='--', color=color[ii], marker=shape[ii])
240 # ax2.loglog(func_eval_SDC[ii, :][:-1], error_SDC['vel'][value, ii, :][:-1], ls='solid', color=color[ii], marker=shape[ii], label='k={}'.format(jj))
241 # ax2.loglog(func_eval_picard[ii,:][:-1], error_picard['vel'][value, ii, :][:-1], ls='--', color=color[ii], marker=shape[ii])
242 #
243 # =============================================================================
244 ax1.loglog(
245 func_eval_SDC[ii, :],
246 error_SDC['pos'][value, ii, :],
247 ls='solid',
248 color=color[ii],
249 marker=shape[ii],
250 label='k={}'.format(jj),
251 )
252 ax1.loglog(
253 func_eval_picard[ii, :], error_picard['pos'][value, ii, :], ls='--', color=color[ii], marker=shape[ii]
254 )
256 ax2.loglog(
257 func_eval_SDC[ii, :],
258 error_SDC['vel'][value, ii, :],
259 ls='solid',
260 color=color[ii],
261 marker=shape[ii],
262 label='k={}'.format(jj),
263 )
264 ax2.loglog(
265 func_eval_picard[ii, :], error_picard['vel'][value, ii, :], ls='--', color=color[ii], marker=shape[ii]
266 )
268 xmin = np.min(ax1.get_xlim())
269 xmax = np.max(ax2.get_xlim())
270 xmin = round(xmin, -3)
271 xmax = round(xmax, -3)
273 xx = np.linspace(np.log(xmin), np.log(xmax), 5)
274 xx = 3**xx
275 xx = xx[np.where(xx < xmax)]
276 # xx=[2*1e+3,4*1e+3, 8*1e+3]
277 ax1.grid(True)
279 ax1.set_title("$M={}$".format(self.num_nodes))
280 ax1.set_xlabel("Number of RHS evaluations")
281 ax1.set_ylabel(r'$\Delta x^{\mathrm{(rel)}}_{%d}$' % (value + 1))
282 ax1.loglog([], [], color="black", ls="--", label="Picard iteration")
283 ax1.loglog([], [], color="black", ls="solid", label="Boris-SDC iteration")
285 ax1.set_xticks(xx)
286 ax1.xaxis.set_major_formatter(self.format_number)
287 ax1.set_ylim(np.min(ax1.get_ylim()), np.max(ax2.get_ylim()))
288 # ax1.set_ylim(1e-14, 1e+0)
290 ax1.legend(loc="best", fontsize=12)
291 fig1.tight_layout()
292 fig1.savefig(self.cwd + "data/f_eval_pos_{}_M={}.pdf".format(value, self.num_nodes))
294 ax2.grid(True)
295 ax2.xaxis.set_major_formatter(self.format_number)
296 ax2.set_title("$M={}$".format(self.num_nodes))
297 ax2.set_xlabel("Number of RHS evaluations")
298 ax2.set_ylabel(r'$\Delta v^{\mathrm{(rel)}}_{%d}$' % (value + 1))
299 ax2.loglog([], [], color="black", ls="--", label="Picard iteration")
300 ax2.loglog([], [], color="black", ls="solid", label="Boris-SDC iteration")
301 ax2.set_xticks(xx)
302 ax2.xaxis.set_major_formatter(self.format_number)
303 ax2.set_ylim(np.min(ax1.get_ylim()), np.max(ax2.get_ylim()))
304 # ax2.set_ylim(1e-14, 1e+0)
305 ax2.legend(loc="best", fontsize=12)
306 fig2.tight_layout()
307 fig2.savefig(self.cwd + "data/f_eval_vel_{}_M={}.pdf".format(value, self.num_nodes))
308 plt.show()
310 def organize_data(self, filename='data/dt_vs_local_errorSDC.csv', time_iter=None):
311 """
312 Organize data according to plot
313 Args:
314 filename (string): data to find approximate order
315 time_iter : in case it you used different time iterations
316 """
317 if time_iter is None:
318 time_iter = self.time_iter
320 items = np.genfromtxt(filename, delimiter=',', skip_header=1)
321 time = items[:, 0]
322 N = int(np.size(time) / time_iter)
324 error_data = {'pos': np.zeros([3, N, time_iter]), 'vel': np.zeros([3, N, time_iter])}
325 order_data = {'pos': np.zeros([N, time_iter, 2]), 'vel': np.zeros([N, time_iter, 2])}
326 time_data = np.zeros([N, time_iter])
327 convline = {'pos': np.zeros([3, N, time_iter]), 'vel': np.zeros([3, N, time_iter])}
329 time_data = time.reshape([N, time_iter])
331 order_data['pos'][:, :, 0] = items[:, 1].reshape([N, time_iter])
332 order_data['pos'][:, :, 1] = items[:, 2].reshape([N, time_iter])
333 order_data['vel'][:, :, 0] = items[:, 6].reshape([N, time_iter])
334 order_data['vel'][:, :, 1] = items[:, 7].reshape([N, time_iter])
336 for ii in range(0, 3):
337 error_data['pos'][ii, :, :] = items[:, ii + 3].reshape([N, time_iter])
338 error_data['vel'][ii, :, :] = items[:, ii + 8].reshape([N, time_iter])
340 for jj in range(0, 3):
341 if jj == 2:
342 convline['pos'][jj, :, :] = (
343 (time_data / time_data[0, 0]).T ** order_data['pos'][:, jj, 1]
344 ).T * error_data['pos'][jj, :, 0][:, None]
345 convline['vel'][jj, :, :] = (
346 (time_data / time_data[0, 0]).T ** order_data['vel'][:, jj, 1]
347 ).T * error_data['vel'][jj, :, 0][:, None]
348 else:
349 convline['pos'][jj, :, :] = (
350 (time_data / time_data[0, 0]).T ** order_data['pos'][:, jj, 0]
351 ).T * error_data['pos'][jj, :, 0][:, None]
352 convline['vel'][jj, :, :] = (
353 (time_data / time_data[0, 0]).T ** order_data['vel'][:, jj, 0]
354 ).T * error_data['vel'][jj, :, 0][:, None]
356 return [N, time_data, error_data, order_data, convline]
358 # find approximate order
359 def find_approximate_order(self, filename='data/dt_vs_local_errorSDC.csv'):
360 """
361 This function finds approximate convergence rate and saves in the data folder
362 Args:
363 filename: given data
364 return:
365 None
366 """
367 [N, time_data, error_data, order_data, convline] = self.organize_data(self.cwd + filename)
368 approx_order = {'pos': np.zeros([1, N]), 'vel': np.zeros([1, N])}
370 for jj in range(0, 3):
371 if jj == 0:
372 file = open(self.cwd + 'data/{}_order_vs_approx_order.csv'.format(self.error_type), 'w')
374 else:
375 file = open(self.cwd + 'data/{}_order_vs_approx_order.csv'.format(self.error_type), 'a')
377 for ii in range(0, N):
378 approx_order['pos'][0, ii] = np.polyfit(
379 np.log(time_data[ii, :]), np.log(error_data['pos'][jj, ii, :]), 1
380 )[0].real
381 approx_order['vel'][0, ii] = np.polyfit(
382 np.log(time_data[ii, :]), np.log(error_data['vel'][jj, ii, :]), 1
383 )[0].real
384 if jj == 2:
385 file.write(
386 str(order_data['pos'][:, jj, 1])
387 + ' | '
388 + str(approx_order['pos'][0])
389 + ' | '
390 + str(order_data['vel'][:, jj, 1])
391 + ' | '
392 + str(approx_order['vel'][0])
393 + '\n'
394 )
395 else:
396 file.write(
397 str(order_data['pos'][:, jj, 0])
398 + ' | '
399 + str(approx_order['pos'][0])
400 + ' | '
401 + str(order_data['vel'][:, jj, 0])
402 + ' | '
403 + str(approx_order['vel'][0])
404 + '\n'
405 )
406 file.close()