Coverage for pySDC/tutorial/step_5/C_advection_and_PFASST.py: 100%

1from pathlib import Path 

2import numpy as np 

3 

4from pySDC.helpers.stats_helper import get_sorted 

5 

6from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

7from pySDC.implementations.problem_classes.AdvectionEquation_ND_FD import advectionNd 

8from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit 

9from pySDC.implementations.transfer_classes.TransferMesh import mesh_to_mesh 

10 

11 

12def main(): 

13 """ 

14 A simple test program to run PFASST for the advection equation in multiple ways... 

15 """ 

16 

17 # initialize level parameters 

18 level_params = dict() 

19 level_params['restol'] = 1e-09 

20 level_params['dt'] = 0.0625 

21 

22 # initialize sweeper parameters 

23 sweeper_params = dict() 

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

25 sweeper_params['num_nodes'] = [3] 

26 

27 # initialize problem parameters 

28 problem_params = dict() 

29 problem_params['c'] = 1 # advection coefficient 

30 problem_params['freq'] = 4 # frequency for the test value 

31 problem_params['nvars'] = [128, 64] # number of degrees of freedom for each level 

32 problem_params['order'] = 4 

33 problem_params['bc'] = 'periodic' 

34 problem_params['stencil_type'] = 'center' 

35 

36 # initialize step parameters 

37 step_params = dict() 

38 step_params['maxiter'] = 50 

39 

40 # initialize space transfer parameters 

41 space_transfer_params = dict() 

42 space_transfer_params['rorder'] = 2 

43 space_transfer_params['iorder'] = 6 

44 space_transfer_params['periodic'] = True 

45 

46 # initialize controller parameters 

47 controller_params = dict() 

48 controller_params['logger_level'] = 30 

49 controller_params['predict_type'] = 'pfasst_burnin' 

50 

51 # fill description dictionary for easy step instantiation 

52 description = dict() 

53 description['problem_class'] = advectionNd # pass problem class 

54 description['problem_params'] = problem_params # pass problem parameters 

55 description['sweeper_class'] = generic_implicit # pass sweeper (see part B) 

56 description['level_params'] = level_params # pass level parameters 

57 description['step_params'] = step_params # pass step parameters 

58 description['space_transfer_class'] = mesh_to_mesh # pass spatial transfer class 

59 description['space_transfer_params'] = space_transfer_params # pass parameters for spatial transfer 

60 

61 # set time parameters 

62 t0 = 0.0 

63 Tend = 1.0 

64 

65 # set up list of parallel time-steps to run PFASST with 

66 nsteps = int(Tend / level_params['dt']) 

67 num_proc_list = [2**i for i in range(int(np.log2(nsteps) + 1))] 

68 

69 # set up list of types of implicit SDC sweepers: LU and implicit Euler here 

70 QI_list = ['LU', 'IE'] 

71 niters_min_all = {} 

72 niters_max_all = {} 

73 

74 Path("data").mkdir(parents=True, exist_ok=True) 

75 f = open('data/step_5_C_out.txt', 'w') 

76 # loop over different types of implicit sweeper types 

77 for QI in QI_list: 

78 # define and set preconditioner for the implicit sweeper 

79 sweeper_params['QI'] = QI 

80 description['sweeper_params'] = sweeper_params # pass sweeper parameters 

81 

82 # init min/max iteration counts 

83 niters_min_all[QI] = 99 

84 niters_max_all[QI] = 0 

85 

86 # loop over different number of processes 

87 for num_proc in num_proc_list: 

88 out = 'Working with QI = %s on %2i processes...' % (QI, num_proc) 

89 f.write(out + '\n') 

90 print(out) 

91 # instantiate controller 

92 controller = controller_nonMPI( 

93 num_procs=num_proc, controller_params=controller_params, description=description 

94 ) 

95 

96 # get initial values on finest level 

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

98 uinit = P.u_exact(t0) 

99 

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

101 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend) 

102 

103 # compute exact solution and compare 

104 uex = P.u_exact(Tend) 

105 err = abs(uex - uend) 

106 

107 # filter statistics by type (number of iterations) 

108 iter_counts = get_sorted(stats, type='niter', sortby='time') 

109 

110 # compute and print statistics 

111 niters = np.array([item[1] for item in iter_counts]) 

112 niters_min_all[QI] = min(np.mean(niters), niters_min_all[QI]) 

113 niters_max_all[QI] = max(np.mean(niters), niters_max_all[QI]) 

114 out = ' Mean number of iterations: %4.2f' % np.mean(niters) 

115 f.write(out + '\n') 

116 print(out) 

117 out = ' Range of values for number of iterations: %2i ' % np.ptp(niters) 

118 f.write(out + '\n') 

119 print(out) 

120 out = ' Position of max/min number of iterations: %2i -- %2i' % ( 

121 int(np.argmax(niters)), 

122 int(np.argmin(niters)), 

123 ) 

124 f.write(out + '\n') 

125 print(out) 

126 out = ' Std and var for number of iterations: %4.2f -- %4.2f' % ( 

127 float(np.std(niters)), 

128 float(np.var(niters)), 

129 ) 

130 f.write(out + '\n') 

131 f.write(out + '\n') 

132 print(out) 

133 

134 f.write('\n') 

135 print() 

136 

137 assert err < 5.1365e-04, "ERROR: error is too high, got %s" % err 

138 

139 out = 'Mean number of iterations went up from %4.2f to %4.2f for QI = %s!' % ( 

140 niters_min_all[QI], 

141 niters_max_all[QI], 

142 QI, 

143 ) 

144 f.write(out + '\n') 

145 print(out) 

146 

147 f.write('\n\n') 

148 print() 

149 print() 

150 

151 f.close() 

152 

153 

154if __name__ == "__main__": 

155 main()