Coverage for pySDC/projects/DAE/run/run_iteration_test.py: 100%

1import numpy as np 

2import statistics 

3import pickle 

4 

5from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

6from pySDC.projects.DAE.problems.simpleDAE import SimpleDAE 

7from pySDC.projects.DAE.sweepers.fullyImplicitDAE import FullyImplicitDAE 

8from pySDC.projects.DAE.misc.hooksDAE import LogGlobalErrorPostStepDifferentialVariable 

9from pySDC.helpers.stats_helper import get_sorted 

10from pySDC.helpers.stats_helper import filter_stats 

11 

12 

13def setup(): 

14 """ 

15 Routine to initialise iteration test parameters 

16 """ 

17 

18 # initialize level parameters 

19 level_params = dict() 

20 level_params['restol'] = -1 

21 # Used for generating the first set of plots. Chose this because in the convergence plots the three collocation methods investigated had converged. Maybe too big? -> actually looked at results for different step sizes. There was no real difference. 

22 # level_params['dt'] = 1e-3 

23 level_params['dt'] = 1e-4 

24 

25 # This comes as read-in for the sweeper class 

26 sweeper_params = dict() 

27 sweeper_params['quad_type'] = 'RADAU-RIGHT' 

28 

29 # This comes as read-in for the problem class 

30 problem_params = dict() 

31 # Absolute termination tollerance for implicit solver 

32 # Exactly how this is used can be adjusted in update_nodes() in the fully implicit sweeper 

33 problem_params['newton_tol'] = 1e-7 

34 

35 # This comes as read-in for the step class 

36 step_params = dict() 

37 

38 # initialize controller parameters 

39 controller_params = dict() 

40 controller_params['logger_level'] = 30 

41 controller_params['hook_class'] = LogGlobalErrorPostStepDifferentialVariable 

42 

43 # Fill description dictionary for easy hierarchy creation 

44 description = dict() 

45 description['problem_class'] = SimpleDAE 

46 description['problem_params'] = problem_params 

47 description['sweeper_class'] = FullyImplicitDAE 

48 description['sweeper_params'] = sweeper_params 

49 description['level_params'] = level_params 

50 description['step_params'] = step_params 

51 

52 # set simulation parameters 

53 run_params = dict() 

54 run_params['t0'] = 0.0 

55 run_params['tend'] = 0.1 

56 max_iter_low = 4 

57 max_iter_high = 6 

58 run_params['max_iter_list'] = list(range(max_iter_low, max_iter_high)) 

59 run_params['qd_list'] = ['IE', 'LU'] 

60 run_params['num_nodes_list'] = [3] 

61 

62 return description, controller_params, run_params 

63 

64 

65def run(description, controller_params, run_params): 

66 """ 

67 Routine to run simulation 

68 """ 

69 conv_data = dict() 

70 

71 for qd_type in run_params['qd_list']: 

72 description['sweeper_params']['QI'] = qd_type 

73 conv_data[qd_type] = dict() 

74 

75 for num_nodes in run_params['num_nodes_list']: 

76 description['sweeper_params']['num_nodes'] = num_nodes 

77 conv_data[qd_type][num_nodes] = dict() 

78 conv_data[qd_type][num_nodes]['error'] = np.zeros_like(run_params['max_iter_list'], dtype=float) 

79 conv_data[qd_type][num_nodes]['residual'] = np.zeros_like(run_params['max_iter_list'], dtype=float) 

80 conv_data[qd_type][num_nodes]['niter'] = np.zeros_like(run_params['max_iter_list'], dtype='int') 

81 conv_data[qd_type][num_nodes]['max_iter'] = run_params['max_iter_list'] 

82 

83 for i, max_iter in enumerate(run_params['max_iter_list']): 

84 print('Working on Qdelta=%s -- num. nodes=%i -- max. iter.=%i' % (qd_type, num_nodes, max_iter)) 

85 description['step_params']['maxiter'] = max_iter 

86 

87 # instantiate the controller 

88 controller = controller_nonMPI( 

89 num_procs=1, controller_params=controller_params, description=description 

90 ) 

91 # get initial values 

92 P = controller.MS[0].levels[0].prob 

93 uinit = P.u_exact(run_params['t0']) 

94 

95 # call main function to get things done... 

96 uend, stats = controller.run(u0=uinit, t0=run_params['t0'], Tend=run_params['tend']) 

97 

98 # compute exact solution and compare 

99 err = get_sorted(stats, type='e_global_differential_post_step', sortby='time') 

100 residual = get_sorted(stats, type='residual_post_step', sortby='time') 

101 niter = filter_stats(stats, type='niter') 

102 

103 conv_data[qd_type][num_nodes]['error'][i] = np.linalg.norm([err[j][1] for j in range(len(err))], np.inf) 

104 conv_data[qd_type][num_nodes]['residual'][i] = np.linalg.norm( 

105 [residual[j][1] for j in range(len(residual))], np.inf 

106 ) 

107 conv_data[qd_type][num_nodes]['niter'][i] = round(statistics.mean(niter.values())) 

108 print( 

109 "Error=", 

110 conv_data[qd_type][num_nodes]['error'][i], 

111 " Residual=", 

112 conv_data[qd_type][num_nodes]['residual'][i], 

113 ) 

114 

115 return conv_data 

116 

117 

118if __name__ == "__main__": 

119 """ 

120 Routine to run simple differential-algebraic-equation example with various max iters, preconditioners and collocation node counts 

121 In contrast to run_convergence_test.py, in which max iters is set large enough to not be the limiting factor, max iters is varied for a fixed time step and the improvement in the error is measured 

122 Error data is stored in a dictionary and then pickled for use with the loglog_plot.py routine 

123 """ 

124 

125 description, controller_params, run_params = setup() 

126 conv_data = run(description, controller_params, run_params) 

127 pickle.dump(conv_data, open("data/dae_iter_data.p", 'wb')) 

128 print("Done")