Coverage for pySDC/projects/TOMS/AllenCahn_contracting_circle.py: 95%
185 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 20:28 +0000
1import os
3import dill
4import matplotlib.ticker as ticker
5import numpy as np
7import pySDC.helpers.plot_helper as plt_helper
8from pySDC.helpers.stats_helper import get_sorted
10from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
11from pySDC.implementations.problem_classes.AllenCahn_2D_FD import (
12 allencahn_fullyimplicit,
13 allencahn_semiimplicit,
14 allencahn_semiimplicit_v2,
15 allencahn_multiimplicit,
16 allencahn_multiimplicit_v2,
17)
18from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
19from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
20from pySDC.implementations.sweeper_classes.multi_implicit import multi_implicit
21from pySDC.implementations.hooks.AllenCahn_monitor import AllenCahnMonitor
23# http://www.personal.psu.edu/qud2/Res/Pre/dz09sisc.pdf
26def setup_parameters():
27 """
28 Helper routine to fill in all relevant parameters
30 Note that this file will be used for all versions of SDC, containing more than necessary for each individual run
32 Returns:
33 description (dict)
34 controller_params (dict)
35 """
37 # initialize level parameters
38 level_params = dict()
39 level_params['restol'] = 1e-08
40 level_params['dt'] = 1e-03
41 level_params['nsweeps'] = [1]
43 # initialize sweeper parameters
44 sweeper_params = dict()
45 sweeper_params['quad_type'] = 'RADAU-RIGHT'
46 sweeper_params['num_nodes'] = [3]
47 sweeper_params['Q1'] = ['LU']
48 sweeper_params['Q2'] = ['LU']
49 sweeper_params['QI'] = ['LU']
50 sweeper_params['QE'] = ['EE']
51 sweeper_params['initial_guess'] = 'zero'
53 # This comes as read-in for the problem class
54 problem_params = dict()
55 problem_params['nvars'] = [(128, 128)]
56 problem_params['eps'] = [0.04]
57 problem_params['newton_maxiter'] = 100
58 problem_params['newton_tol'] = 1e-09
59 problem_params['lin_tol'] = 1e-10
60 problem_params['lin_maxiter'] = 100
61 problem_params['radius'] = 0.25
63 # initialize step parameters
64 step_params = dict()
65 step_params['maxiter'] = 50
67 # initialize controller parameters
68 controller_params = dict()
69 controller_params['logger_level'] = 30
70 controller_params['hook_class'] = AllenCahnMonitor
72 # fill description dictionary for easy step instantiation
73 description = dict()
74 description['problem_class'] = None # pass problem class
75 description['problem_params'] = problem_params # pass problem parameters
76 description['sweeper_class'] = None # pass sweeper (see part B)
77 description['sweeper_params'] = sweeper_params # pass sweeper parameters
78 description['level_params'] = level_params # pass level parameters
79 description['step_params'] = step_params # pass step parameters
81 return description, controller_params
84def run_SDC_variant(variant=None, inexact=False):
85 """
86 Routine to run particular SDC variant
88 Args:
89 variant (str): string describing the variant
90 inexact (bool): flag to use inexact nonlinear solve (or nor)
92 Returns:
93 timing (float)
94 niter (float)
95 """
97 # load (incomplete) default parameters
98 description, controller_params = setup_parameters()
100 # add stuff based on variant
101 if variant == 'fully-implicit':
102 description['problem_class'] = allencahn_fullyimplicit
103 description['sweeper_class'] = generic_implicit
104 if inexact:
105 description['problem_params']['newton_maxiter'] = 1
106 elif variant == 'semi-implicit':
107 description['problem_class'] = allencahn_semiimplicit
108 description['sweeper_class'] = imex_1st_order
109 if inexact:
110 description['problem_params']['lin_maxiter'] = 10
111 elif variant == 'semi-implicit_v2':
112 description['problem_class'] = allencahn_semiimplicit_v2
113 description['sweeper_class'] = imex_1st_order
114 if inexact:
115 description['problem_params']['newton_maxiter'] = 1
116 elif variant == 'multi-implicit':
117 description['problem_class'] = allencahn_multiimplicit
118 description['sweeper_class'] = multi_implicit
119 if inexact:
120 description['problem_params']['newton_maxiter'] = 1
121 description['problem_params']['lin_maxiter'] = 10
122 elif variant == 'multi-implicit_v2':
123 description['problem_class'] = allencahn_multiimplicit_v2
124 description['sweeper_class'] = multi_implicit
125 if inexact:
126 description['problem_params']['newton_maxiter'] = 1
127 else:
128 raise NotImplementedError('Wrong variant specified, got %s' % variant)
130 if inexact:
131 out = 'Working on inexact %s variant...' % variant
132 else:
133 out = 'Working on exact %s variant...' % variant
134 print(out)
136 # setup parameters "in time"
137 t0 = 0
138 Tend = 0.032
140 # instantiate controller
141 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
143 # get initial values on finest level
144 P = controller.MS[0].levels[0].prob
145 uinit = P.u_exact(t0)
147 # call main function to get things done...
148 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
150 # filter statistics by variant (number of iterations)
151 iter_counts = get_sorted(stats, type='niter', sortby='time')
153 # compute and print statistics
154 niters = np.array([item[1] for item in iter_counts])
155 out = ' Mean number of iterations: %4.2f' % np.mean(niters)
156 print(out)
157 out = ' Range of values for number of iterations: %2i ' % np.ptp(niters)
158 print(out)
159 out = ' Position of max/min number of iterations: %2i -- %2i' % (int(np.argmax(niters)), int(np.argmin(niters)))
160 print(out)
161 out = ' Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters)))
162 print(out)
164 newton_iters = P.work_counters['newton'].niter
165 lin_iters = P.work_counters['linear'].niter
166 print(' Iteration count (nonlinear/linear): %i / %i' % (newton_iters, lin_iters))
167 print(
168 ' Mean Iteration count per call: %4.2f / %4.2f'
169 % (newton_iters / max(P.newton_ncalls, 1), lin_iters / max(P.lin_ncalls, 1))
170 )
172 timing = get_sorted(stats, type='timing_run', sortby='time')
174 print('Time to solution: %6.4f sec.' % timing[0][1])
175 print()
177 return stats
180def show_results(fname, cwd=''):
181 """
182 Plotting routine
184 Args:
185 fname (str): file name to read in and name plots
186 cwd (str): current working directory
187 """
189 file = open(cwd + fname + '.pkl', 'rb')
190 results = dill.load(file)
191 file.close()
193 # plt_helper.mpl.style.use('classic')
194 plt_helper.setup_mpl()
196 # set up plot for timings
197 fig, ax1 = plt_helper.newfig(textwidth=238.96, scale=1.5, ratio=0.4)
199 timings = {}
200 niters = {}
201 for key, item in results.items():
202 timings[key] = get_sorted(item, type='timing_run', sortby='time')[0][1]
203 iter_counts = get_sorted(item, type='niter', sortby='time')
204 niters[key] = np.mean(np.array([item[1] for item in iter_counts]))
206 xcoords = list(range(len(timings)))
207 sorted_timings = sorted([(key, timings[key]) for key in timings], reverse=True, key=lambda tup: tup[1])
208 sorted_niters = [(k, niters[k]) for k in [key[0] for key in sorted_timings]]
209 heights_timings = [item[1] for item in sorted_timings]
210 heights_niters = [item[1] for item in sorted_niters]
211 keys = [(item[0][1] + ' ' + item[0][0]).replace('-', '\n').replace('_v2', ' mod.') for item in sorted_timings]
213 ax1.bar(xcoords, heights_timings, align='edge', width=-0.3, label='timings (left axis)')
214 ax1.set_ylabel('time (sec)')
216 ax2 = ax1.twinx()
217 ax2.bar(xcoords, heights_niters, color='lightcoral', align='edge', width=0.3, label='iterations (right axis)')
218 ax2.set_ylabel('mean number of iterations')
220 ax1.set_xticks(xcoords)
221 ax1.set_xticklabels(keys, rotation=90, ha='center')
223 # ask matplotlib for the plotted objects and their labels
224 lines, labels = ax1.get_legend_handles_labels()
225 lines2, labels2 = ax2.get_legend_handles_labels()
226 ax2.legend(lines + lines2, labels + labels2, loc=0)
228 # save plot, beautify
229 f = fname + '_timings'
230 plt_helper.savefig(f)
232 assert os.path.isfile(f + '.pdf'), 'ERROR: plotting did not create PDF file'
233 # assert os.path.isfile(f + '.pgf'), 'ERROR: plotting did not create PGF file'
234 assert os.path.isfile(f + '.png'), 'ERROR: plotting did not create PNG file'
236 # set up plot for radii
237 fig, ax = plt_helper.newfig(textwidth=238.96, scale=1.0)
239 exact_radii = []
240 for key, item in results.items():
241 computed_radii = get_sorted(item, type='computed_radius', sortby='time')
243 xcoords = [item0[0] for item0 in computed_radii]
244 radii = [item0[1] for item0 in computed_radii]
245 if key[0] + ' ' + key[1] == 'fully-implicit exact':
246 ax.plot(xcoords, radii, label=(key[0] + ' ' + key[1]).replace('_v2', ' mod.'))
248 exact_radii = get_sorted(item, type='exact_radius', sortby='time')
250 diff = np.array([abs(item0[1] - item1[1]) for item0, item1 in zip(exact_radii, computed_radii, strict=True)])
251 max_pos = int(np.argmax(diff))
252 assert max(diff) < 0.07, 'ERROR: computed radius is too far away from exact radius, got %s' % max(diff)
253 assert 0.028 < computed_radii[max_pos][0] < 0.03, (
254 'ERROR: largest difference is at wrong time, got %s' % computed_radii[max_pos][0]
255 )
257 xcoords = [item[0] for item in exact_radii]
258 radii = [item[1] for item in exact_radii]
259 ax.plot(xcoords, radii, color='k', linestyle='--', linewidth=1, label='exact')
261 ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%1.2f'))
262 ax.set_ylabel('radius')
263 ax.set_xlabel('time')
264 ax.grid()
265 ax.legend(loc=3)
267 # save plot, beautify
268 f = fname + '_radii'
269 plt_helper.savefig(f)
271 assert os.path.isfile(f + '.pdf'), 'ERROR: plotting did not create PDF file'
272 # assert os.path.isfile(f + '.pgf'), 'ERROR: plotting did not create PGF file'
273 assert os.path.isfile(f + '.png'), 'ERROR: plotting did not create PNG file'
275 # set up plot for interface width
276 fig, ax = plt_helper.newfig(textwidth=238.96, scale=1.0)
278 interface_width = []
279 for key, item in results.items():
280 interface_width = get_sorted(item, type='interface_width', sortby='time')
281 xcoords = [item[0] for item in interface_width]
282 width = [item[1] for item in interface_width]
283 if key[0] + ' ' + key[1] == 'fully-implicit exact':
284 ax.plot(xcoords, width, label=key[0] + ' ' + key[1])
286 xcoords = [item[0] for item in interface_width]
287 init_width = [interface_width[0][1]] * len(xcoords)
288 ax.plot(xcoords, init_width, color='k', linestyle='--', linewidth=1, label='exact')
290 ax.yaxis.set_major_formatter(ticker.FormatStrFormatter('%1.2f'))
291 ax.set_ylabel(r'interface width ($\epsilon$)')
292 ax.set_xlabel('time')
293 ax.grid()
294 ax.legend(loc=3)
296 # save plot, beautify
297 f = fname + '_interface'
298 plt_helper.savefig(f)
300 assert os.path.isfile(f + '.pdf'), 'ERROR: plotting did not create PDF file'
301 # assert os.path.isfile(f + '.pgf'), 'ERROR: plotting did not create PGF file'
302 assert os.path.isfile(f + '.png'), 'ERROR: plotting did not create PNG file'
304 return None
307def main(cwd=''):
308 """
309 Main driver
311 Args:
312 cwd (str): current working directory (need this for testing)
313 """
315 # Loop over variants, exact and inexact solves
316 results = {}
317 for variant in ['multi-implicit', 'semi-implicit', 'fully-implicit', 'semi-implicit_v2', 'multi-implicit_v2']:
318 results[(variant, 'exact')] = run_SDC_variant(variant=variant, inexact=False)
319 results[(variant, 'inexact')] = run_SDC_variant(variant=variant, inexact=True)
321 # dump result
322 fname = 'data/results_SDC_variants_AllenCahn_1E-03'
323 file = open(cwd + fname + '.pkl', 'wb')
324 dill.dump(results, file)
325 file.close()
326 assert os.path.isfile(cwd + fname + '.pkl'), 'ERROR: dill did not create file'
328 # visualize
329 show_results(fname, cwd=cwd)
332if __name__ == "__main__":
333 main()