Coverage for pySDC/projects/parallelSDC/preconditioner_playground.py: 97%

142 statements  

« prev     ^ index     » next       coverage.py v7.6.7, created at 2024-11-16 14:51 +0000

1import os 

2import pickle 

3from collections import namedtuple 

4 

5import numpy as np 

6 

7import pySDC.helpers.plot_helper as plt_helper 

8from pySDC.helpers.stats_helper import get_sorted 

9 

10from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

11from pySDC.implementations.problem_classes.AdvectionEquation_ND_FD import advectionNd 

12from pySDC.implementations.problem_classes.GeneralizedFisher_1D_FD_implicit import generalized_fisher 

13from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced 

14from pySDC.implementations.problem_classes.Van_der_Pol_implicit import vanderpol 

15from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit 

16 

17ID = namedtuple('ID', ['setup', 'qd_type', 'param']) 

18 

19 

20def main(): 

21 # initialize level parameters (part I) 

22 level_params = dict() 

23 level_params['restol'] = 1e-08 

24 

25 # initialize sweeper parameters (part I) 

26 sweeper_params = dict() 

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

28 sweeper_params['num_nodes'] = 3 

29 

30 # initialize step parameters 

31 step_params = dict() 

32 step_params['maxiter'] = 100 

33 

34 # initialize controller parameters 

35 controller_params = dict() 

36 controller_params['logger_level'] = 30 

37 

38 # set up list of Q-delta types and setups 

39 qd_list = ['LU', 'IE', 'IEpar', 'Qpar', 'MIN', 'MIN3', 'MIN_GT'] 

40 setup_list = [ 

41 ('heat', 63, [10.0**i for i in range(-3, 3)]), 

42 ('advection', 64, [10.0**i for i in range(-3, 3)]), 

43 ('vanderpol', 2, [0.1 * 2**i for i in range(0, 10)]), 

44 ('fisher', 63, [2**i for i in range(-2, 3)]), 

45 ] 

46 # setup_list = [('fisher', 63, [2 * i for i in range(1, 6)])] 

47 

48 # pre-fill results with lists of setups 

49 results = dict() 

50 for setup, nvars, param_list in setup_list: 

51 results[setup] = (nvars, param_list) 

52 

53 # loop over all Q-delta matrix types 

54 for qd_type in qd_list: 

55 # assign implicit Q-delta matrix 

56 sweeper_params['QI'] = qd_type 

57 

58 # loop over all setups 

59 for setup, nvars, param_list in setup_list: 

60 # initialize problem parameters (part I) 

61 problem_params = dict() 

62 if setup != 'vanderpol': 

63 problem_params['nvars'] = nvars # number of degrees of freedom for each level 

64 

65 # loop over all parameters 

66 for param in param_list: 

67 # fill description for the controller 

68 description = dict() 

69 description['sweeper_class'] = generic_implicit # pass sweeper 

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

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

72 

73 print('working on: %s - %s - %s' % (qd_type, setup, param)) 

74 

75 # decide which setup to take 

76 if setup == 'heat': 

77 problem_params['nu'] = param 

78 problem_params['freq'] = 2 

79 problem_params['bc'] = 'dirichlet-zero' # boundary conditions 

80 

81 level_params['dt'] = 0.1 

82 

83 description['problem_class'] = heatNd_unforced 

84 description['problem_params'] = problem_params 

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

86 

87 elif setup == 'advection': 

88 problem_params['c'] = param 

89 problem_params['order'] = 2 

90 problem_params['freq'] = 2 

91 problem_params['stencil_type'] = 'center' # boundary conditions 

92 problem_params['bc'] = 'periodic' # boundary conditions 

93 

94 level_params['dt'] = 0.1 

95 

96 description['problem_class'] = advectionNd 

97 description['problem_params'] = problem_params 

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

99 

100 elif setup == 'vanderpol': 

101 problem_params['newton_tol'] = 1e-09 

102 problem_params['newton_maxiter'] = 20 

103 problem_params['mu'] = param 

104 problem_params['u0'] = np.array([2.0, 0]) 

105 

106 level_params['dt'] = 0.1 

107 

108 description['problem_class'] = vanderpol 

109 description['problem_params'] = problem_params 

110 description['level_params'] = level_params 

111 

112 elif setup == 'fisher': 

113 problem_params['nu'] = 1 

114 problem_params['lambda0'] = param 

115 problem_params['newton_maxiter'] = 20 

116 problem_params['newton_tol'] = 1e-10 

117 problem_params['interval'] = (-5, 5) 

118 

119 level_params['dt'] = 0.01 

120 

121 description['problem_class'] = generalized_fisher 

122 description['problem_params'] = problem_params 

123 description['level_params'] = level_params 

124 

125 else: 

126 print('Setup not implemented..', setup) 

127 exit() 

128 

129 # instantiate controller 

130 controller = controller_nonMPI( 

131 num_procs=1, controller_params=controller_params, description=description 

132 ) 

133 

134 # get initial values on finest level 

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

136 uinit = P.u_exact(0) 

137 

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

139 uend, stats = controller.run(u0=uinit, t0=0, Tend=level_params['dt']) 

140 

141 # filter statistics by type (number of iterations) 

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

143 

144 # just one time-step, grep number of iteration and store 

145 niter = iter_counts[0][1] 

146 id = ID(setup=setup, qd_type=qd_type, param=param) 

147 results[id] = niter 

148 

149 assert len(results) == (6 + 6 + 10 + 5) * 7 + 4, 'ERROR: did not get all results, got %s' % len(results) 

150 

151 # write out for later visualization 

152 file = open('data/parallelSDC_iterations_precond.pkl', 'wb') 

153 pickle.dump(results, file) 

154 

155 assert os.path.isfile('data/parallelSDC_iterations_precond.pkl'), 'ERROR: pickle did not create file' 

156 

157 

158def plot_iterations(): 

159 """ 

160 Helper routine to plot iteration counts 

161 """ 

162 

163 file = open('data/parallelSDC_iterations_precond.pkl', 'rb') 

164 results = pickle.load(file) 

165 

166 # find the lists/header required for plotting 

167 qd_type_list = [] 

168 setup_list = [] 

169 for key in results.keys(): 

170 if isinstance(key, ID): 

171 if key.qd_type not in qd_type_list: 

172 qd_type_list.append(key.qd_type) 

173 elif isinstance(key, str): 

174 setup_list.append(key) 

175 print('Found these type of preconditioners:', qd_type_list) 

176 print('Found these setups:', setup_list) 

177 

178 assert len(qd_type_list) == 7, 'ERROR did not find five preconditioners, got %s' % qd_type_list 

179 assert len(setup_list) == 4, 'ERROR: did not find three setup, got %s' % setup_list 

180 

181 qd_type_list = ['LU', 'IE', 'IEpar', 'Qpar', 'MIN', 'MIN3', 'MIN_GT'] 

182 marker_list = [None, None, 's', 'o', '^', 'd', 'x'] 

183 color_list = ['k', 'k', 'r', 'g', 'b', 'c', 'm'] 

184 

185 plt_helper.setup_mpl() 

186 

187 # loop over setups and Q-delta types: one figure per setup, all Qds in one plot 

188 for setup in setup_list: 

189 plt_helper.newfig(textwidth=238.96, scale=0.89) 

190 

191 for qd_type, marker, color in zip(qd_type_list, marker_list, color_list): 

192 niter = np.zeros(len(results[setup][1])) 

193 for key in results.keys(): 

194 if isinstance(key, ID): 

195 if key.setup == setup and key.qd_type == qd_type: 

196 xvalue = results[setup][1].index(key.param) 

197 niter[xvalue] = results[key] 

198 if qd_type == 'LU': 

199 ls = '--' 

200 lw = 0.5 

201 elif qd_type == 'IE': 

202 ls = '-.' 

203 lw = 0.5 

204 else: 

205 ls = '-' 

206 lw = 1 

207 plt_helper.plt.semilogx( 

208 results[setup][1], 

209 niter, 

210 label=qd_type, 

211 lw=lw, 

212 linestyle=ls, 

213 color=color, 

214 marker=marker, 

215 markeredgecolor='k', 

216 ) 

217 

218 if setup == 'heat': 

219 xlabel = r'$\nu$' 

220 elif setup == 'advection': 

221 xlabel = r'$c$' 

222 elif setup == 'fisher': 

223 xlabel = r'$\lambda_0$' 

224 elif setup == 'vanderpol': 

225 xlabel = r'$\mu$' 

226 else: 

227 print('Setup not implemented..', setup) 

228 exit() 

229 

230 plt_helper.plt.ylim([0, 60]) 

231 plt_helper.plt.legend(loc=2, ncol=1) 

232 plt_helper.plt.ylabel('number of iterations') 

233 plt_helper.plt.xlabel(xlabel) 

234 plt_helper.plt.grid() 

235 

236 # save plot as PDF and PGF 

237 fname = 'data/parallelSDC_preconditioner_' + setup 

238 plt_helper.savefig(fname) 

239 

240 assert os.path.isfile(fname + '.pdf'), 'ERROR: plotting did not create PDF file' 

241 # assert os.path.isfile(fname + '.pgf'), 'ERROR: plotting did not create PGF file' 

242 assert os.path.isfile(fname + '.png'), 'ERROR: plotting did not create PNG file' 

243 

244 

245if __name__ == "__main__": 

246 main() 

247 plot_iterations()