Coverage for pySDC/projects/Hamiltonian/fput.py: 100%

1import os 

2from collections import defaultdict 

3 

4import dill 

5import numpy as np 

6 

7import pySDC.helpers.plot_helper as plt_helper 

8from pySDC.helpers.stats_helper import get_sorted, filter_stats 

9 

10from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

11from pySDC.implementations.problem_classes.FermiPastaUlamTsingou import fermi_pasta_ulam_tsingou 

12from pySDC.implementations.sweeper_classes.verlet import verlet 

13from pySDC.implementations.transfer_classes.TransferParticles_NoCoarse import particles_to_particles 

14from pySDC.projects.Hamiltonian.hamiltonian_and_energy_output import hamiltonian_and_energy_output 

15 

16 

17def setup_fput(): 

18 """ 

19 Helper routine for setting up everything for the Fermi-Pasta-Ulam-Tsingou problem 

20 

21 Returns: 

22 description (dict): description of the controller 

23 controller_params (dict): controller parameters 

24 """ 

25 

26 # initialize level parameters 

27 level_params = dict() 

28 level_params['restol'] = 1e-12 

29 level_params['dt'] = 2.0 

30 

31 # initialize sweeper parameters 

32 sweeper_params = dict() 

33 sweeper_params['quad_type'] = 'LOBATTO' 

34 sweeper_params['num_nodes'] = [5, 3] 

35 sweeper_params['initial_guess'] = 'zero' 

36 

37 # initialize problem parameters for the Penning trap 

38 problem_params = dict() 

39 problem_params['npart'] = 2048 

40 problem_params['alpha'] = 0.25 

41 problem_params['k'] = 1.0 

42 problem_params['energy_modes'] = [[1, 2, 3, 4]] 

43 

44 # initialize step parameters 

45 step_params = dict() 

46 step_params['maxiter'] = 50 

47 

48 # initialize controller parameters 

49 controller_params = dict() 

50 controller_params['hook_class'] = hamiltonian_and_energy_output 

51 controller_params['logger_level'] = 30 

52 

53 # Fill description dictionary for easy hierarchy creation 

54 description = dict() 

55 description['problem_class'] = fermi_pasta_ulam_tsingou 

56 description['problem_params'] = problem_params 

57 description['sweeper_class'] = verlet 

58 description['sweeper_params'] = sweeper_params 

59 description['level_params'] = level_params 

60 description['step_params'] = step_params 

61 description['space_transfer_class'] = particles_to_particles 

62 

63 return description, controller_params 

64 

65 

66def run_simulation(): 

67 """ 

68 Routine to run the simulation of a second order problem 

69 

70 """ 

71 

72 description, controller_params = setup_fput() 

73 # set time parameters 

74 t0 = 0.0 

75 # set this to 10000 to reproduce the picture in 

76 # http://www.scholarpedia.org/article/Fermi-Pasta-Ulam_nonlinear_lattice_oscillations 

77 Tend = 1000.0 

78 num_procs = 1 

79 

80 f = open('data/fput_out.txt', 'w') 

81 out = 'Running fput problem with %s processors...' % num_procs 

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

83 print(out) 

84 

85 # instantiate the controller 

86 controller = controller_nonMPI(num_procs=num_procs, controller_params=controller_params, description=description) 

87 

88 # get initial values on finest level 

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

90 uinit = P.u_exact(t=t0) 

91 

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

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

94 

95 # filter statistics by type (number of iterations) 

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

97 

98 # compute and print statistics 

99 # for item in iter_counts: 

100 # out = 'Number of iterations for time %4.2f: %2i' % item 

101 # f.write(out + '\n') 

102 # print(out) 

103 

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

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

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

107 print(out) 

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

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

110 print(out) 

111 out = ' Position of max/min number of iterations: %2i -- %2i' % (int(np.argmax(niters)), int(np.argmin(niters))) 

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

113 print(out) 

114 out = ' Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters))) 

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

116 print(out) 

117 

118 # get runtime 

119 timing_run = get_sorted(stats, type='timing_run', sortby='time')[0][1] 

120 out = '... took %6.4f seconds to run this.' % timing_run 

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

122 print(out) 

123 f.close() 

124 

125 # measured 5.828 (2914 iterations in 500 steps, 3 to 10 per step); the tolerance allows about 50 iterations to 

126 # change with round-off, since restol=1e-12 is close to it, but not one more (or less) iteration per step 

127 mean_niters = np.mean(niters) 

128 assert np.isclose(mean_niters, 5.828, atol=0.1), 'Mean number of iterations has changed, got %s' % mean_niters 

129 

130 fname = 'data/fput.dat' 

131 f = open(fname, 'wb') 

132 dill.dump(stats, f) 

133 f.close() 

134 

135 assert os.path.isfile(fname), 'Run for %s did not create stats file' 

136 

137 

138def show_results(cwd=''): 

139 """ 

140 Helper function to plot the error of the Hamiltonian 

141 

142 Args: 

143 cwd (str): current working directory 

144 """ 

145 

146 # read in the dill data 

147 f = open(cwd + 'data/fput.dat', 'rb') 

148 stats = dill.load(f) 

149 f.close() 

150 

151 plt_helper.mpl.style.use('classic') 

152 plt_helper.setup_mpl() 

153 

154 # HAMILTONIAN PLOTTING # 

155 # extract error in hamiltonian and prepare for plotting 

156 extract_stats = filter_stats(stats, type='err_hamiltonian') 

157 result = defaultdict(list) 

158 for k, v in extract_stats.items(): 

159 result[k.iter].append((k.time, v)) 

160 for k, _ in result.items(): 

161 result[k] = sorted(result[k], key=lambda x: x[0]) 

162 

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

164 

165 # Rearrange data for easy plotting 

166 err_ham = 1 

167 for k, v in result.items(): 

168 time = [item[0] for item in v] 

169 ham = [item[1] for item in v] 

170 err_ham = ham[-1] 

171 plt_helper.plt.semilogy(time, ham, '-', lw=1, label='Iter ' + str(k)) 

172 print(err_ham) 

173 # the last error in the highest iteration, measured 1.98e-7, which is also the largest error of any converged step 

174 assert err_ham < 3e-7, 'Error in the Hamiltonian is too large, got %s' % err_ham 

175 

176 plt_helper.plt.xlabel('Time') 

177 plt_helper.plt.ylabel('Error in Hamiltonian') 

178 plt_helper.plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) 

179 

180 fname = 'data/fput_hamiltonian' 

181 plt_helper.savefig(fname) 

182 

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

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

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

186 

187 # ENERGY PLOTTING # 

188 # extract error in hamiltonian and prepare for plotting 

189 result = get_sorted(stats, type='energy_step', sortby='time') 

190 

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

192 

193 # Rearrange data for easy plotting 

194 for mode in result[0][1].keys(): 

195 time = [item[0] for item in result] 

196 energy = [item[1][mode] for item in result] 

197 plt_helper.plt.plot(time, energy, label=str(mode) + 'th mode') 

198 

199 plt_helper.plt.xlabel('Time') 

200 plt_helper.plt.ylabel('Energy') 

201 plt_helper.plt.legend(loc='center left', bbox_to_anchor=(1, 0.5)) 

202 

203 fname = 'data/fput_energy' 

204 plt_helper.savefig(fname) 

205 

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

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

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

209 

210 # POSITION PLOTTING # 

211 # extract positions and prepare for plotting 

212 result = get_sorted(stats, type='position', sortby='time') 

213 

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

215 

216 # Rearrange data for easy plotting 

217 nparts = len(result[0][1]) 

218 nsteps = len(result) 

219 pos = np.zeros((nparts, nsteps)) 

220 time = np.zeros(nsteps) 

221 for idx, item in enumerate(result): 

222 time[idx] = item[0] 

223 for n in range(nparts): 

224 pos[n, idx] = item[1][n] 

225 

226 for n in range(min(nparts, 16)): 

227 plt_helper.plt.plot(time, pos[n, :]) 

228 

229 fname = 'data/fput_positions' 

230 plt_helper.savefig(fname) 

231 

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

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

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

235 

236 

237def main(): 

238 run_simulation() 

239 show_results() 

240 

241 

242if __name__ == "__main__": 

243 main()