Skip to content

Commit d0a49dc

Browse files
pancettapancetta
committed
added tutorial and test for iteration estimator
1 parent b64032b commit d0a49dc

File tree

7 files changed

+392
-19
lines changed

7 files changed

+392
-19
lines changed

‎docs/source/tutorial/doc_step_8_C.rst‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
Full code: `pySDC/tutorial/step_8/C_iteration_estimator.py <https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/tutorial/step_8/C_iteration_estimator.py>`_
2+
3+
.. literalinclude:: ../../../pySDC/tutorial/step_8/C_iteration_estimator.py
4+
5+
Results:
6+
7+
.. literalinclude:: ../../../step_8_C_out.txt

‎pySDC/implementations/controller_classes/controller_nonMPI.py‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -73,8 +73,7 @@ def __init__(self, num_procs, controller_params, description):
7373
self.logger.warning('you have specified a predictor type but only a single level.. '
7474
'predictor will be ignored')
7575

76-
@staticmethod
77-
def check_iteration_estimator(MS):
76+
def check_iteration_estimator(self, MS):
7877
"""
7978
Method to check the iteration estimator
8079
@@ -99,8 +98,8 @@ def check_iteration_estimator(MS):
9998
S.status.diff_old_loc = diff_new
10099
alpha = 1 / (1 - Ltilde_loc) * S.status.diff_first_loc
101100
Kest_loc = np.log(S.params.errtol / alpha) / np.log(Ltilde_loc) * 1.05 # Safety factor!
102-
print(f'LOCAL: {L.time:8.4f}, {S.status.iter}: {int(np.ceil(Kest_loc))}, {Ltilde_loc:8.6e}, '
103-
f'{Kest_loc:8.6e}, {Ltilde_loc ** S.status.iter * alpha:8.6e}')
101+
self.logger.debug(f'LOCAL: {L.time:8.4f}, {S.status.iter}: {int(np.ceil(Kest_loc))}, '
102+
f'{Ltilde_loc:8.6e}, {Kest_loc:8.6e}, {Ltilde_loc ** S.status.iter * alpha:8.6e}')
104103
# You should not stop prematurely on earlier steps, since later steps may need more accuracy to reach
105104
# the tolerance themselves. The final Kest_loc is the one that counts.
106105
# if np.ceil(Kest_loc) <= S.status.iter:

‎pySDC/playgrounds/mpifft/grayscott.py‎

Lines changed: 17 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -53,11 +53,13 @@ def run_simulation(spectral=None, ml=None, num_procs=None):
5353
# initialize step parameters
5454
step_params = dict()
5555
step_params['maxiter'] = 500
56+
step_params['errtol'] = 1E-07
5657

5758
# initialize controller parameters
5859
controller_params = dict()
5960
controller_params['logger_level'] = 30 if rank == 0 else 99
6061
# controller_params['predict_type'] = 'fine_only'
62+
controller_params['use_iteration_estimator'] = False
6163

6264
# fill description dictionary for easy step instantiation
6365
description = dict()
@@ -71,7 +73,7 @@ def run_simulation(spectral=None, ml=None, num_procs=None):
7173

7274
# set time parameters
7375
t0 = 0.0
74-
Tend = 20
76+
Tend = 3500
7577

7678
f = None
7779
if rank == 0:
@@ -103,19 +105,19 @@ def run_simulation(spectral=None, ml=None, num_procs=None):
103105
# call main function to get things done...
104106
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
105107

106-
# plt.figure()
107-
# plt.imshow(uend[..., 0])#, vmin=0, vmax=1)
108-
# plt.title('u')
109-
# plt.colorbar()
110-
# plt.figure()
111-
# plt.imshow(uend[..., 1])#, vmin=0, vmax=1)
112-
# plt.title('v')
113-
# plt.colorbar()
108+
plt.figure()
109+
plt.imshow(P.fft.backward(uend[..., 0]))#, vmin=0, vmax=1)
110+
plt.title('u')
111+
plt.colorbar()
112+
plt.figure()
113+
plt.imshow(P.fft.backward(uend[..., 1]))#, vmin=0, vmax=1)
114+
plt.title('v')
115+
plt.colorbar()
114116
# plt.figure()
115117
# plt.imshow(uend[..., 0] + uend[..., 1])
116118
# plt.title('sum')
117119
# plt.colorbar()
118-
# plt.show()
120+
plt.show()
119121
# # exit()
120122

121123
if rank == 0:
@@ -157,12 +159,12 @@ def main():
157159
158160
Note: This can also be run with "mpirun -np 2 python grayscott.py"
159161
"""
160-
run_simulation(spectral=False, ml=False, num_procs=1)
162+
# run_simulation(spectral=False, ml=False, num_procs=1)
161163
run_simulation(spectral=True, ml=False, num_procs=1)
162-
run_simulation(spectral=False, ml=True, num_procs=1)
163-
run_simulation(spectral=True, ml=True, num_procs=1)
164-
run_simulation(spectral=False, ml=True, num_procs=10)
165-
run_simulation(spectral=True, ml=True, num_procs=10)
164+
# run_simulation(spectral=False, ml=True, num_procs=1)
165+
# run_simulation(spectral=True, ml=True, num_procs=1)
166+
# run_simulation(spectral=False, ml=True, num_procs=10)
167+
# run_simulation(spectral=True, ml=True, num_procs=10)
166168

167169

168170
if __name__ == "__main__":

‎pySDC/tests/test_tutorials/test_step_8.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,13 @@
11
from pySDC.tutorial.step_8.A_visualize_residuals import main as main_A
22
from pySDC.tutorial.step_8.B_multistep_SDC import main as main_B
3+
from pySDC.tutorial.step_8.C_iteration_estimator import main as main_C
34

45

56
def test_A():
67
main_A()
78

89
def test_B():
910
main_B()
11+
12+
def test_C():
13+
main_C()

‎pySDC/tutorial/step_8/C_iteration_estimator.py‎

Lines changed: 282 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,282 @@
1+
import numpy as np
2+
3+
from pySDC.helpers.stats_helper import filter_stats, sort_stats
4+
from pySDC.implementations.collocation_classes.gauss_radau_right import CollGaussRadau_Right
5+
from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
6+
from pySDC.implementations.problem_classes.HeatEquation_ND_FD_forced_periodic import heatNd_periodic
7+
from pySDC.implementations.problem_classes.AdvectionEquation_ND_FD_periodic import advectionNd_periodic
8+
from pySDC.implementations.problem_classes.Auzinger_implicit import auzinger
9+
from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order
10+
from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
11+
from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh
12+
from pySDC.implementations.transfer_classes.TransferMesh_NoCoarse import mesh_to_mesh as mesh_to_mesh_nc
13+
from pySDC.playgrounds.compression.HookClass_error_output import error_output
14+
15+
16+
def setup_diffusion(dt=None, ndim=None, ml=False):
17+
18+
# initialize level parameters
19+
level_params = dict()
20+
level_params['restol'] = 1E-10
21+
level_params['dt'] = dt # time-step size
22+
level_params['nsweeps'] = 1
23+
24+
# initialize sweeper parameters
25+
sweeper_params = dict()
26+
sweeper_params['collocation_class'] = CollGaussRadau_Right
27+
sweeper_params['num_nodes'] = 3
28+
sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
29+
# sweeper_params['initial_guess'] = 'zero'
30+
31+
# initialize problem parameters
32+
problem_params = dict()
33+
problem_params['ndim'] = ndim # will be iterated over
34+
problem_params['order'] = 8 # order of accuracy for FD discretization in space
35+
problem_params['nu'] = 0.1 # diffusion coefficient
36+
problem_params['freq'] = tuple(2 for _ in range(ndim)) # frequencies
37+
if ml:
38+
problem_params['nvars'] = [tuple(64 for _ in range(ndim)), tuple(32 for _ in range(ndim))] # number of dofs
39+
else:
40+
problem_params['nvars'] = tuple(64 for _ in range(ndim)) # number of dofs
41+
problem_params['direct_solver'] = False # do GMRES instead of LU
42+
problem_params['liniter'] = 10 # number of GMRES iterations
43+
44+
# initialize step parameters
45+
step_params = dict()
46+
step_params['maxiter'] = 50
47+
step_params['errtol'] = 1E-07
48+
49+
# initialize space transfer parameters
50+
space_transfer_params = dict()
51+
space_transfer_params['rorder'] = 2
52+
space_transfer_params['iorder'] = 6
53+
space_transfer_params['periodic'] = True
54+
55+
# initialize controller parameters
56+
controller_params = dict()
57+
controller_params['logger_level'] = 30
58+
controller_params['use_iteration_estimator'] = True
59+
controller_params['hook_class'] = error_output
60+
61+
# fill description dictionary for easy step instantiation
62+
description = dict()
63+
description['problem_class'] = heatNd_periodic # pass problem class
64+
description['problem_params'] = problem_params # pass problem parameters
65+
description['sweeper_class'] = imex_1st_order # pass sweeper (see part B)
66+
description['sweeper_params'] = sweeper_params # pass sweeper parameters
67+
description['level_params'] = level_params # pass level parameters
68+
description['step_params'] = step_params # pass step parameters
69+
if ml:
70+
description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
71+
description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer
72+
73+
return description, controller_params
74+
75+
76+
def setup_advection(dt=None, ndim=None, ml=False):
77+
78+
# initialize level parameters
79+
level_params = dict()
80+
level_params['restol'] = 1E-10
81+
level_params['dt'] = dt # time-step size
82+
level_params['nsweeps'] = 1
83+
84+
# initialize sweeper parameters
85+
sweeper_params = dict()
86+
sweeper_params['collocation_class'] = CollGaussRadau_Right
87+
sweeper_params['num_nodes'] = 3
88+
sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
89+
# sweeper_params['initial_guess'] = 'zero'
90+
91+
# initialize problem parameters
92+
problem_params = dict()
93+
problem_params['ndim'] = ndim # will be iterated over
94+
problem_params['order'] = 6 # order of accuracy for FD discretization in space
95+
problem_params['type'] = 'center' # order of accuracy for FD discretization in space
96+
problem_params['c'] = 0.1 # diffusion coefficient
97+
problem_params['freq'] = tuple(2 for _ in range(ndim)) # frequencies
98+
if ml:
99+
problem_params['nvars'] = [tuple(64 for _ in range(ndim)), tuple(32 for _ in range(ndim))] # number of dofs
100+
else:
101+
problem_params['nvars'] = tuple(64 for _ in range(ndim)) # number of dofs
102+
problem_params['direct_solver'] = False # do GMRES instead of LU
103+
problem_params['liniter'] = 10 # number of GMRES iterations
104+
105+
# initialize step parameters
106+
step_params = dict()
107+
step_params['maxiter'] = 50
108+
step_params['errtol'] = 1E-07
109+
110+
# initialize space transfer parameters
111+
space_transfer_params = dict()
112+
space_transfer_params['rorder'] = 2
113+
space_transfer_params['iorder'] = 6
114+
space_transfer_params['periodic'] = True
115+
116+
# initialize controller parameters
117+
controller_params = dict()
118+
controller_params['logger_level'] = 30
119+
controller_params['use_iteration_estimator'] = True
120+
controller_params['hook_class'] = error_output
121+
122+
# fill description dictionary for easy step instantiation
123+
description = dict()
124+
description['problem_class'] = advectionNd_periodic
125+
description['problem_params'] = problem_params # pass problem parameters
126+
description['sweeper_class'] = generic_implicit
127+
description['sweeper_params'] = sweeper_params # pass sweeper parameters
128+
description['level_params'] = level_params # pass level parameters
129+
description['step_params'] = step_params # pass step parameters
130+
if ml:
131+
description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class
132+
description['space_transfer_params'] = space_transfer_params # pass paramters for spatial transfer
133+
134+
return description, controller_params
135+
136+
137+
def setup_auzinger(dt=None, ml=False):
138+
139+
# initialize level parameters
140+
level_params = dict()
141+
level_params['restol'] = 1E-10
142+
level_params['dt'] = dt # time-step size
143+
level_params['nsweeps'] = 1
144+
145+
# initialize sweeper parameters
146+
sweeper_params = dict()
147+
sweeper_params['collocation_class'] = CollGaussRadau_Right
148+
if ml:
149+
sweeper_params['num_nodes'] = [3, 2]
150+
else:
151+
sweeper_params['num_nodes'] = 3
152+
sweeper_params['QI'] = ['LU'] # For the IMEX sweeper, the LU-trick can be activated for the implicit part
153+
# sweeper_params['initial_guess'] = 'zero'
154+
155+
# initialize problem parameters
156+
problem_params = dict()
157+
problem_params['newton_tol'] = 1E-12
158+
problem_params['newton_maxiter'] = 10
159+
160+
# initialize step parameters
161+
step_params = dict()
162+
step_params['maxiter'] = 50
163+
step_params['errtol'] = 1E-07
164+
165+
# initialize controller parameters
166+
controller_params = dict()
167+
controller_params['logger_level'] = 30
168+
controller_params['use_iteration_estimator'] = True
169+
controller_params['hook_class'] = error_output
170+
171+
# fill description dictionary for easy step instantiation
172+
description = dict()
173+
description['problem_class'] = auzinger
174+
description['problem_params'] = problem_params # pass problem parameters
175+
description['sweeper_class'] = generic_implicit
176+
description['sweeper_params'] = sweeper_params # pass sweeper parameters
177+
description['level_params'] = level_params # pass level parameters
178+
description['step_params'] = step_params # pass step parameters
179+
if ml:
180+
description['space_transfer_class'] = mesh_to_mesh_nc # pass spatial transfer class
181+
182+
return description, controller_params
183+
184+
185+
def run_simulations(type=None, ndim_list=None, Tend=None, nsteps_list=None, ml=False, nprocs=None):
186+
"""
187+
A simple test program to do SDC runs for the heat equation in various dimensions
188+
"""
189+
190+
t0 = None
191+
dt = None
192+
description = None
193+
controller_params = None
194+
195+
f = open('step_8_C_out.txt', 'a')
196+
197+
for ndim in ndim_list:
198+
for nsteps in nsteps_list:
199+
200+
if type == 'diffusion':
201+
# set time parameters
202+
t0 = 0.0
203+
dt = (Tend - t0) / nsteps
204+
description, controller_params = setup_diffusion(dt, ndim, ml)
205+
elif type == 'advection':
206+
# set time parameters
207+
t0 = 0.0
208+
dt = (Tend - t0) / nsteps
209+
description, controller_params = setup_advection(dt, ndim, ml)
210+
elif type == 'auzinger':
211+
assert ndim == 1
212+
# set time parameters
213+
t0 = 0.0
214+
dt = (Tend - t0) / nsteps
215+
description, controller_params = setup_auzinger(dt, ml)
216+
217+
out = f'Running {type} in {ndim} dimensions with time-step size {dt}...\n'
218+
f.write(out + '\n')
219+
print(out)
220+
221+
# Warning: this is black magic used to run an 'exact' collocation solver for each step within the hooks
222+
description['step_params']['description'] = description
223+
description['step_params']['controller_params'] = controller_params
224+
225+
# instantiate controller
226+
controller = controller_nonMPI(num_procs=nprocs, controller_params=controller_params,
227+
description=description)
228+
229+
# get initial values on finest level
230+
P = controller.MS[0].levels[0].prob
231+
uinit = P.u_exact(t0)
232+
233+
# call main function to get things done...
234+
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
235+
236+
# filter statistics by type (number of iterations)
237+
iter_counts = sort_stats(filter_stats(stats, type='niter'), sortby='time')
238+
239+
niters = np.array([item[1] for item in iter_counts])
240+
out = f' Mean number of iterations: {np.mean(niters):4.2f}'
241+
f.write(out + '\n')
242+
print(out)
243+
244+
# filter statistics by type (error after time-step)
245+
PDE_errors = sort_stats(filter_stats(stats, type='PDE_error_after_step'), sortby='time')
246+
coll_errors = sort_stats(filter_stats(stats, type='coll_error_after_step'), sortby='time')
247+
for iters, PDE_err, coll_err in zip(iter_counts, PDE_errors, coll_errors):
248+
assert coll_err[1] < description['step_params']['errtol'], f'Error too high, got {coll_err[1]:8.4e}'
249+
out = f' Errors after step {PDE_err[0]:8.4f} with {iters[1]} iterations: ' \
250+
f'{PDE_err[1]:8.4e} / {coll_err[1]:8.4e}'
251+
f.write(out + '\n')
252+
print(out)
253+
f.write('\n')
254+
print()
255+
256+
# filter statistics by type (error after time-step)
257+
timing = sort_stats(filter_stats(stats, type='timing_run'), sortby='time')
258+
out = f'...done, took {timing[0][1]} seconds!'
259+
f.write(out + '\n')
260+
print(out)
261+
262+
print()
263+
out = '-----------------------------------------------------------------------------'
264+
f.write(out + '\n')
265+
print(out)
266+
267+
f.close()
268+
269+
270+
def main():
271+
run_simulations(type='diffusion', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=False, nprocs=1)
272+
run_simulations(type='diffusion', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=True, nprocs=1)
273+
274+
run_simulations(type='advection', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=False, nprocs=1)
275+
run_simulations(type='advection', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=True, nprocs=1)
276+
277+
run_simulations(type='auzinger', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=False, nprocs=1)
278+
run_simulations(type='auzinger', ndim_list=[1], Tend=1.0, nsteps_list=[8], ml=True, nprocs=1)
279+
280+
281+
if __name__ == "__main__":
282+
main()

0 commit comments

Comments
 (0)