Coverage for pySDC/tutorial/step_8/C_iteration_estimator.py: 100%
179 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1import numpy as np
2from pathlib import Path
4from pySDC.helpers.stats_helper import get_sorted
6from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
7from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_forced
8from pySDC.implementations.problem_classes.AdvectionEquation_ND_FD import advectionNd
9from pySDC.implementations.problem_classes.Auzinger_implicit import auzinger
10from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
11from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
12from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
13from pySDC.implementations.transfer_classes.TransferMesh_NoCoarse import mesh_to_mesh as mesh_to_mesh_nc
14from pySDC.implementations.convergence_controller_classes.check_iteration_estimator import CheckIterationEstimatorNonMPI
15from pySDC.tutorial.step_8.HookClass_error_output import error_output
18def setup_diffusion(dt=None, ndim=None, ml=False):
19 # initialize level parameters
20 level_params = dict()
21 level_params['restol'] = 1e-10
22 level_params['dt'] = dt # time-step size
23 level_params['nsweeps'] = 1
25 # initialize sweeper parameters
26 sweeper_params = dict()
27 sweeper_params['quad_type'] = 'RADAU-RIGHT'
28 sweeper_params['num_nodes'] = 3
29 sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
30 # sweeper_params['initial_guess'] = 'zero'
32 # initialize problem parameters
33 problem_params = dict()
34 problem_params['order'] = 8 # order of accuracy for FD discretization in space
35 problem_params['nu'] = 0.1 # diffusion coefficient
36 problem_params['bc'] = 'periodic' # boundary conditions
37 problem_params['freq'] = tuple(2 for _ in range(ndim)) # frequencies
38 if ml:
39 problem_params['nvars'] = [tuple(64 for _ in range(ndim)), tuple(32 for _ in range(ndim))] # number of dofs
40 else:
41 problem_params['nvars'] = tuple(64 for _ in range(ndim)) # number of dofs
42 problem_params['solver_type'] = 'CG' # do CG instead of LU
43 problem_params['liniter'] = 10 # number of CG iterations
45 # initialize step parameters
46 step_params = dict()
47 step_params['maxiter'] = 50
48 step_params['errtol'] = 1e-07
50 # initialize space transfer parameters
51 space_transfer_params = dict()
52 space_transfer_params['rorder'] = 2
53 space_transfer_params['iorder'] = 6
54 space_transfer_params['periodic'] = True
56 # setup the iteration estimator
57 convergence_controllers = dict()
58 convergence_controllers[CheckIterationEstimatorNonMPI] = {'errtol': 1e-7}
60 # initialize controller parameters
61 controller_params = dict()
62 controller_params['logger_level'] = 30
63 controller_params['hook_class'] = error_output
65 # fill description dictionary for easy step instantiation
66 description = dict()
67 description['problem_class'] = heatNd_forced # pass problem class
68 description['problem_params'] = problem_params # pass problem parameters
69 description['sweeper_class'] = imex_1st_order # pass sweeper (see part B)
70 description['sweeper_params'] = sweeper_params # pass sweeper parameters
71 description['level_params'] = level_params # pass level parameters
72 description['step_params'] = step_params # pass step parameters
73 description['convergence_controllers'] = convergence_controllers
74 if ml:
75 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
76 description['space_transfer_params'] = space_transfer_params # pass parameters for spatial transfer
78 return description, controller_params
81def setup_advection(dt=None, ndim=None, ml=False):
82 # initialize level parameters
83 level_params = dict()
84 level_params['restol'] = 1e-10
85 level_params['dt'] = dt # time-step size
86 level_params['nsweeps'] = 1
88 # initialize sweeper parameters
89 sweeper_params = dict()
90 sweeper_params['quad_type'] = 'RADAU-RIGHT'
91 sweeper_params['num_nodes'] = 3
92 sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
93 # sweeper_params['initial_guess'] = 'zero'
95 # initialize problem parameters
96 problem_params = dict()
97 problem_params['order'] = 6 # order of accuracy for FD discretization in space
98 problem_params['stencil_type'] = 'center' # order of accuracy for FD discretization in space
99 problem_params['bc'] = 'periodic' # boundary conditions
100 problem_params['c'] = 0.1 # diffusion coefficient
101 problem_params['freq'] = tuple(2 for _ in range(ndim)) # frequencies
102 if ml:
103 problem_params['nvars'] = [tuple(64 for _ in range(ndim)), tuple(32 for _ in range(ndim))] # number of dofs
104 else:
105 problem_params['nvars'] = tuple(64 for _ in range(ndim)) # number of dofs
106 problem_params['solver_type'] = 'GMRES' # do GMRES instead of LU
107 problem_params['liniter'] = 10 # number of GMRES iterations
109 # initialize step parameters
110 step_params = dict()
111 step_params['maxiter'] = 50
112 step_params['errtol'] = 1e-07
114 # initialize space transfer parameters
115 space_transfer_params = dict()
116 space_transfer_params['rorder'] = 2
117 space_transfer_params['iorder'] = 6
118 space_transfer_params['periodic'] = True
120 # setup the iteration estimator
121 convergence_controllers = dict()
122 convergence_controllers[CheckIterationEstimatorNonMPI] = {'errtol': 1e-7}
124 # initialize controller parameters
125 controller_params = dict()
126 controller_params['logger_level'] = 30
127 controller_params['hook_class'] = error_output
129 # fill description dictionary for easy step instantiation
130 description = dict()
131 description['problem_class'] = advectionNd
132 description['problem_params'] = problem_params # pass problem parameters
133 description['sweeper_class'] = generic_implicit
134 description['sweeper_params'] = sweeper_params # pass sweeper parameters
135 description['level_params'] = level_params # pass level parameters
136 description['step_params'] = step_params # pass step parameters
137 description['convergence_controllers'] = convergence_controllers
138 if ml:
139 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
140 description['space_transfer_params'] = space_transfer_params # pass parameters for spatial transfer
142 return description, controller_params
145def setup_auzinger(dt=None, ml=False):
146 # initialize level parameters
147 level_params = dict()
148 level_params['restol'] = 1e-10
149 level_params['dt'] = dt # time-step size
150 level_params['nsweeps'] = 1
152 # initialize sweeper parameters
153 sweeper_params = dict()
154 sweeper_params['quad_type'] = 'RADAU-RIGHT'
155 if ml:
156 sweeper_params['num_nodes'] = [3, 2]
157 else:
158 sweeper_params['num_nodes'] = 3
159 sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
160 # sweeper_params['initial_guess'] = 'zero'
162 # initialize problem parameters
163 problem_params = dict()
164 problem_params['newton_tol'] = 1e-12
165 problem_params['newton_maxiter'] = 10
167 # initialize step parameters
168 step_params = dict()
169 step_params['maxiter'] = 50
170 step_params['errtol'] = 1e-07
172 # setup the iteration estimator
173 convergence_controllers = dict()
174 convergence_controllers[CheckIterationEstimatorNonMPI] = {'errtol': 1e-7}
176 # initialize controller parameters
177 controller_params = dict()
178 controller_params['logger_level'] = 30
179 controller_params['hook_class'] = error_output
181 # fill description dictionary for easy step instantiation
182 description = dict()
183 description['problem_class'] = auzinger
184 description['problem_params'] = problem_params # pass problem parameters
185 description['sweeper_class'] = generic_implicit
186 description['sweeper_params'] = sweeper_params # pass sweeper parameters
187 description['level_params'] = level_params # pass level parameters
188 description['step_params'] = step_params # pass step parameters
189 description['convergence_controllers'] = convergence_controllers
190 if ml:
191 description['space_transfer_class'] = mesh_to_mesh_nc # pass spatial transfer class
193 return description, controller_params
196def run_simulations(type=None, ndim_list=None, Tend=None, nsteps_list=None, ml=False, nprocs=None):
197 """
198 A simple test program to do SDC runs for the heat equation in various dimensions
199 """
201 t0 = None
202 dt = None
203 description = None
204 controller_params = None
206 Path("data").mkdir(parents=True, exist_ok=True)
207 f = open('data/step_8_C_out.txt', 'a')
209 for ndim in ndim_list:
210 for nsteps in nsteps_list:
211 if type == 'diffusion':
212 # set time parameters
213 t0 = 0.0
214 dt = (Tend - t0) / nsteps
215 description, controller_params = setup_diffusion(dt, ndim, ml)
216 mean_number_of_iterations = 3.00 if ml else 5.75
217 elif type == 'advection':
218 # set time parameters
219 t0 = 0.0
220 dt = (Tend - t0) / nsteps
221 description, controller_params = setup_advection(dt, ndim, ml)
222 mean_number_of_iterations = 2.00 if ml else 4.00
223 elif type == 'auzinger':
224 assert ndim == 1
225 # set time parameters
226 t0 = 0.0
227 dt = (Tend - t0) / nsteps
228 description, controller_params = setup_auzinger(dt, ml)
229 mean_number_of_iterations = 3.62 if ml else 5.62
231 out = f'Running {type} in {ndim} dimensions with time-step size {dt}...\n'
232 f.write(out + '\n')
233 print(out)
235 # Warning: this is black magic used to run an 'exact' collocation solver for each step within the hooks
236 description['step_params']['description'] = description
237 description['step_params']['controller_params'] = controller_params
239 # instantiate controller
240 controller = controller_nonMPI(
241 num_procs=nprocs, controller_params=controller_params, description=description
242 )
244 # get initial values on finest level
245 P = controller.MS[0].levels[0].prob
246 uinit = P.u_exact(t0)
248 # call main function to get things done...
249 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
251 # filter statistics by type (number of iterations)
252 iter_counts = get_sorted(stats, type='niter', sortby='time')
254 niters = np.array([item[1] for item in iter_counts])
255 out = f' Mean number of iterations: {np.mean(niters):4.2f}'
256 f.write(out + '\n')
257 print(out)
259 # filter statistics by type (error after time-step)
260 PDE_errors = get_sorted(stats, type='PDE_error_after_step', sortby='time')
261 coll_errors = get_sorted(stats, type='coll_error_after_step', sortby='time')
262 for iters, PDE_err, coll_err in zip(iter_counts, PDE_errors, coll_errors):
263 assert coll_err[1] < description['step_params']['errtol'], f'Error too high, got {coll_err[1]:8.4e}'
264 out = (
265 f' Errors after step {PDE_err[0]:8.4f} with {iters[1]} iterations: '
266 f'{PDE_err[1]:8.4e} / {coll_err[1]:8.4e}'
267 )
268 f.write(out + '\n')
269 print(out)
270 f.write('\n')
271 print()
273 # filter statistics by type (error after time-step)
274 timing = get_sorted(stats, type='timing_run', sortby='time')
275 out = f'...done, took {timing[0][1]} seconds!'
276 f.write(out + '\n')
277 print(out)
279 print()
280 out = '-----------------------------------------------------------------------------'
281 f.write(out + '\n')
282 print(out)
284 f.close()
285 assert np.isclose(
286 mean_number_of_iterations, np.mean(niters), atol=1e-2
287 ), f' \
288 Expected \
289{mean_number_of_iterations:.2f} mean iterations, but got {np.mean(niters):.2f}'
292def main():
293 run_simulations(type='diffusion', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=False, nprocs=1)
294 run_simulations(type='diffusion', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=True, nprocs=1)
296 run_simulations(type='advection', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=False, nprocs=1)
297 run_simulations(type='advection', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=True, nprocs=1)
299 run_simulations(type='auzinger', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=False, nprocs=1)
300 run_simulations(type='auzinger', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=True, nprocs=1)
303if __name__ == "__main__":
304 main()