Coverage for pySDC/projects/SDC_showdown/SDC_timing_GrayScott.py: 79%
146 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
1import os
2import pickle
4import numpy as np
5from petsc4py import PETSc
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.GrayScott_2D_PETSc_periodic import (
12 petsc_grayscott_multiimplicit,
13 petsc_grayscott_fullyimplicit,
14 petsc_grayscott_semiimplicit,
15)
16from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
17from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
18from pySDC.implementations.sweeper_classes.multi_implicit import multi_implicit
21def setup_parameters():
22 """
23 Helper routine to fill in all relevant parameters
25 Note that this file will be used for all versions of SDC, containing more than necessary for each individual run
27 Returns:
28 description (dict)
29 controller_params (dict)
30 """
32 # initialize level parameters
33 level_params = dict()
34 level_params['restol'] = 1e-08
35 level_params['dt'] = 1.0
36 level_params['nsweeps'] = [1]
38 # initialize sweeper parameters
39 sweeper_params = dict()
40 sweeper_params['quad_type'] = 'RADAU-RIGHT'
41 sweeper_params['num_nodes'] = [3]
42 sweeper_params['Q1'] = ['LU']
43 sweeper_params['Q2'] = ['LU']
44 sweeper_params['QI'] = ['LU']
45 sweeper_params['initial_guess'] = 'zero'
47 # initialize problem parameters
48 problem_params = dict()
49 problem_params['Du'] = 1.0
50 problem_params['Dv'] = 0.01
51 problem_params['A'] = 0.09
52 problem_params['B'] = 0.086
53 problem_params['nvars'] = [(128, 128)]
54 problem_params['nlsol_tol'] = 1e-10
55 problem_params['nlsol_maxiter'] = 100
56 problem_params['lsol_tol'] = 1e-10
57 problem_params['lsol_maxiter'] = 100
59 # initialize step parameters
60 step_params = dict()
61 step_params['maxiter'] = 50
63 # initialize space transfer parameters
64 # space_transfer_params = dict()
65 # space_transfer_params['finter'] = True
67 # initialize controller parameters
68 controller_params = dict()
69 controller_params['logger_level'] = 30
71 # fill description dictionary for easy step instantiation
72 description = dict()
73 description['problem_class'] = None # pass problem class
74 description['problem_params'] = problem_params # pass problem parameters
75 description['sweeper_class'] = None # pass sweeper (see part B)
76 description['sweeper_params'] = sweeper_params # pass sweeper parameters
77 description['level_params'] = level_params # pass level parameters
78 description['step_params'] = step_params # pass step parameters
79 # description['space_transfer_class'] = mesh_to_mesh_petsc_dmda # pass spatial transfer class
80 # description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer
82 return description, controller_params
85def run_SDC_variant(variant=None, inexact=False, cwd=''):
86 """
87 Routine to run particular SDC variant
89 Args:
90 variant (str): string describing the variant
91 inexact (bool): flag to use inexact nonlinear solve (or nor)
92 cwd (str): current working directory
94 Returns:
95 timing (float)
96 niter (float)
97 """
99 # load (incomplete) default parameters
100 description, controller_params = setup_parameters()
102 # add stuff based on variant
103 if variant == 'fully-implicit':
104 description['problem_class'] = petsc_grayscott_fullyimplicit
105 description['sweeper_class'] = generic_implicit
106 elif variant == 'semi-implicit':
107 description['problem_class'] = petsc_grayscott_semiimplicit
108 description['sweeper_class'] = imex_1st_order
109 elif variant == 'multi-implicit':
110 description['problem_class'] = petsc_grayscott_multiimplicit
111 description['sweeper_class'] = multi_implicit
112 else:
113 raise NotImplementedError('Wrong variant specified, got %s' % variant)
115 if inexact:
116 description['problem_params']['lsol_maxiter'] = 2
117 description['problem_params']['nlsol_maxiter'] = 1
118 out = 'Working on inexact %s variant...' % variant
119 else:
120 out = 'Working on exact %s variant...' % variant
121 print(out)
123 # set time parameters
124 t0 = 0.0
125 Tend = 1.0
127 # instantiate controller
128 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
130 # get initial values on finest level
131 P = controller.MS[0].levels[0].prob
132 uinit = P.u_exact(t0)
134 # call main function to get things done...
135 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
137 # load reference solution to compare with
138 fname = cwd + 'data/GS_reference.dat'
139 viewer = PETSc.Viewer().createBinary(fname, 'r')
140 uex = P.u_exact(t0)
141 uex[:] = PETSc.Vec().load(viewer)
142 err = abs(uex - uend)
144 # filter statistics by variant (number of iterations)
145 iter_counts = get_sorted(stats, type='niter', sortby='time')
147 # compute and print statistics
148 niters = np.array([item[1] for item in iter_counts])
149 out = ' Mean number of iterations: %4.2f' % np.mean(niters)
150 print(out)
151 out = ' Range of values for number of iterations: %2i ' % np.ptp(niters)
152 print(out)
153 out = ' Position of max/min number of iterations: %2i -- %2i' % (int(np.argmax(niters)), int(np.argmin(niters)))
154 print(out)
155 out = ' Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters)))
156 print(out)
158 print('Iteration count (nonlinear/linear): %i / %i' % (P.snes_itercount, P.ksp_itercount))
159 print(
160 'Mean Iteration count per call: %4.2f / %4.2f'
161 % (P.snes_itercount / max(P.snes_ncalls, 1), P.ksp_itercount / max(P.ksp_ncalls, 1))
162 )
164 timing = get_sorted(stats, type='timing_run', sortby='time')
166 print('Time to solution: %6.4f sec.' % timing[0][1])
167 print('Error vs. reference solution: %6.4e' % err)
168 print()
170 assert err < 3e-06, 'ERROR: variant %s did not match error tolerance, got %s' % (variant, err)
171 assert np.mean(niters) <= 10, 'ERROR: number of iterations is too high, got %s' % np.mean(niters)
173 return timing[0][1], np.mean(niters)
176def show_results(fname):
177 """
178 Plotting routine
180 Args:
181 fname: file name to read in and name plots
182 """
184 file = open(fname + '.pkl', 'rb')
185 results = pickle.load(file)
186 file.close()
188 plt_helper.mpl.style.use('classic')
189 plt_helper.setup_mpl()
191 plt_helper.newfig(textwidth=238.96, scale=1.0)
193 xcoords = list(range(len(results)))
194 sorted_data = sorted([(key, results[key][0]) for key in results], reverse=True, key=lambda tup: tup[1])
195 heights = [item[1] for item in sorted_data]
196 keys = [(item[0][1] + ' ' + item[0][0]).replace('-', '\n') for item in sorted_data]
198 plt_helper.plt.bar(xcoords, heights, align='center')
200 plt_helper.plt.xticks(xcoords, keys, rotation=90)
201 plt_helper.plt.ylabel('time (sec)')
203 # save plot, beautify
204 plt_helper.savefig(fname)
206 assert os.path.isfile(fname + '.pdf'), 'ERROR: plotting did not create PDF file'
207 # assert os.path.isfile(fname + '.pgf'), 'ERROR: plotting did not create PGF file'
208 assert os.path.isfile(fname + '.png'), 'ERROR: plotting did not create PNG file'
210 return None
213def run_reference():
214 """
215 Helper routine to create a reference solution using very high order SDC and small time-steps
216 """
218 description, controller_params = setup_parameters()
220 description['problem_class'] = petsc_grayscott_semiimplicit
221 description['sweeper_class'] = imex_1st_order
222 description['sweeper_params']['num_nodes'] = 9
223 description['level_params']['dt'] = 0.01
225 # set time parameters
226 t0 = 0.0
227 Tend = 1.0
229 # instantiate controller
230 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
232 # get initial values on finest level
233 P = controller.MS[0].levels[0].prob
234 uinit = P.u_exact(t0)
236 # call main function to get things done...
237 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
239 # filter statistics by variant (number of iterations)
240 iter_counts = get_sorted(stats, type='niter', sortby='time')
242 # compute and print statistics
243 niters = np.array([item[1] for item in iter_counts])
244 out = ' Mean number of iterations: %4.2f' % np.mean(niters)
245 print(out)
246 out = ' Range of values for number of iterations: %2i ' % np.ptp(niters)
247 print(out)
248 out = ' Position of max/min number of iterations: %2i -- %2i' % (int(np.argmax(niters)), int(np.argmin(niters)))
249 print(out)
250 out = ' Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters)))
251 print(out)
253 print('Iteration count (nonlinear/linear): %i / %i' % (P.snes_itercount, P.ksp_itercount))
254 print(
255 'Mean Iteration count per call: %4.2f / %4.2f'
256 % (P.snes_itercount / max(P.snes_ncalls, 1), P.ksp_itercount / max(P.ksp_ncalls, 1))
257 )
259 timing = get_sorted(stats, type='timing_run', sortby='time')
261 print('Time to solution: %6.4f sec.' % timing[0][1])
263 fname = 'data/GS_reference.dat'
264 viewer = PETSc.Viewer().createBinary(fname, 'w')
265 viewer.view(uend)
267 assert os.path.isfile(fname), 'ERROR: PETSc did not create file'
269 return None
272def main(cwd=''):
273 """
274 Main driver
276 Args:
277 cwd (str): current working directory (need this for testing)
278 """
280 # Loop over variants, exact and inexact solves
281 results = {}
282 for variant in ['fully-implicit', 'multi-implicit', 'semi-implicit']:
283 results[(variant, 'exact')] = run_SDC_variant(variant=variant, inexact=False, cwd=cwd)
284 results[(variant, 'inexact')] = run_SDC_variant(variant=variant, inexact=True, cwd=cwd)
286 # dump result
287 fname = 'data/timings_SDC_variants_GrayScott'
288 file = open(fname + '.pkl', 'wb')
289 pickle.dump(results, file)
290 file.close()
291 assert os.path.isfile(fname + '.pkl'), 'ERROR: pickle did not create file'
293 # visualize
294 show_results(fname)
297if __name__ == "__main__":
298 # run_reference()
299 main()