Coverage for pySDC/projects/Second_orderSDC/penningtrap_run_Hamiltonian_error.py: 0%

1# This code checks if the "data" folder exists or not. 

2exec(open("check_data_folder.py").read()) 

3 

4import matplotlib.pyplot as plt 

5import numpy as np 

6 

7from pySDC.helpers.stats_helper import get_sorted 

8 

9from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

10from pySDC.implementations.problem_classes.PenningTrap_3D import penningtrap 

11from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order 

12from pySDC.projects.Second_orderSDC.penningtrap_HookClass import particles_output 

13from pySDC.implementations.sweeper_classes.Runge_Kutta_Nystrom import RKN 

14from pySDC.projects.Second_orderSDC.plot_helper import set_fixed_plot_params 

15 

16 

17def main(dt, tend, maxiter, M, sweeper): # pragma: no cover 

18 """ 

19 Implementation of Hamiltonian error for Harmonic oscillator problem 

20 mu=0 

21 kappa=1 

22 omega=1 

23 Args: 

24 dt: time step 

25 tend: final time 

26 maxiter: maximal iteration 

27 M: Number of quadrature nodes 

28 sweeper: sweeper class 

29 returns: 

30 Ham_err 

31 """ 

32 # initialize level parameters 

33 level_params = dict() 

34 level_params['restol'] = 1e-16 

35 level_params['dt'] = dt 

36 # initialize sweeper parameters 

37 sweeper_params = dict() 

38 sweeper_params['quad_type'] = 'GAUSS' 

39 sweeper_params['num_nodes'] = M 

40 

41 # initialize problem parameters for the Penning trap 

42 problem_params = dict() 

43 problem_params['omega_E'] = 1 

44 problem_params['omega_B'] = 0 

45 problem_params['u0'] = np.array([[0, 0, 0], [0, 0, 1], [1], [1]], dtype=object) 

46 problem_params['nparts'] = 1 

47 problem_params['sig'] = 0.1 

48 # problem_params['Tend'] = 16.0 

49 

50 # initialize step parameters 

51 step_params = dict() 

52 step_params['maxiter'] = maxiter 

53 

54 # initialize controller parameters 

55 controller_params = dict() 

56 controller_params['hook_class'] = particles_output # specialized hook class for more statistics and output 

57 controller_params['logger_level'] = 30 

58 penningtrap.Harmonic_oscillator = True 

59 # Fill description dictionary for easy hierarchy creation 

60 description = dict() 

61 description['problem_class'] = penningtrap 

62 description['problem_params'] = problem_params 

63 description['sweeper_class'] = sweeper 

64 description['sweeper_params'] = sweeper_params 

65 description['level_params'] = level_params 

66 # description['space_transfer_class'] = particles_to_particles # this is only needed for more than 2 levels 

67 description['step_params'] = step_params 

68 

69 # instantiate the controller (no controller parameters used here) 

70 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) 

71 

72 # set time parameters 

73 t0 = 0.0 

74 Tend = tend 

75 

76 # get initial values on finest level 

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

78 uinit = P.u_init() 

79 

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

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

82 

83 sortedlist_stats = get_sorted(stats, type='etot', sortby='time') 

84 

85 # energy = [entry[1] for entry in sortedlist_stats] 

86 H0 = 1 / 2 * (np.dot(uinit.vel[:].T, uinit.vel[:]) + np.dot(uinit.pos[:].T, uinit.pos[:])) 

87 

88 Ham_err = np.ravel([abs(entry[1] - H0) / H0 for entry in sortedlist_stats]) 

89 return Ham_err 

90 

91 

92def plot_Hamiltonian_error(K, M, dt): # pragma: no cover 

93 """ 

94 Plot Hamiltonian Error 

95 Args: 

96 K: list of maxiter 

97 M: number of quadrature nodes 

98 dt: time step 

99 """ 

100 set_fixed_plot_params() 

101 # Define final time 

102 time = 1e6 

103 tn = dt 

104 # Find time nodes 

105 t = np.arange(0, time + tn, tn) 

106 # Get saved data 

107 t_len = len(t) 

108 RKN1_ = np.loadtxt('data/Ham_RKN1.csv', delimiter='*') 

109 SDC2_ = np.loadtxt('data/Ham_SDC{}{}.csv'.format(M, K[0]), delimiter='*') 

110 SDC3_ = np.loadtxt('data/Ham_SDC{}{}.csv'.format(M, K[1]), delimiter='*') 

111 SDC4_ = np.loadtxt('data/Ham_SDC{}{}.csv'.format(M, K[2]), delimiter='*') 

112 # Only save Hamiltonian error 

113 RKN1 = RKN1_[:, 1] 

114 SDC2 = SDC2_[:, 1] 

115 SDC3 = SDC3_[:, 1] 

116 SDC4 = SDC4_[:, 1] 

117 step = 3000 

118 # plot Hamiltonian error 

119 plt.loglog(t[:t_len:step], RKN1[:t_len:step], label='RKN-4', marker='.', linestyle=' ') 

120 plt.loglog(t[:t_len:step], SDC2[:t_len:step], label='K={}'.format(K[0]), marker='s', linestyle=' ') 

121 plt.loglog(t[:t_len:step], SDC3[:t_len:step], label='K={}'.format(K[1]), marker='*', linestyle=' ') 

122 plt.loglog(t[:t_len:step], SDC4[:t_len:step], label='K={}'.format(K[2]), marker='H', linestyle=' ') 

123 plt.xlabel('$\omega \cdot t$') 

124 plt.ylabel('$\Delta H^{\mathrm{(rel)}}$') 

125 plt.ylim(1e-11, 1e-3 + 0.001) 

126 plt.legend(fontsize=15) 

127 plt.tight_layout() 

128 

129 

130if __name__ == "__main__": 

131 """ 

132 Important: 

133 * Every simulation needs to run individually otherwise it may crash at some point. 

134 I don't know why 

135 * All of the data saved in /data folder 

136 """ 

137 K = (2, 3, 4) 

138 M = 3 

139 dt = 2 * np.pi / 10 

140 tend = 2 * np.pi * 1e6 

141 

142 Ham_SDC1 = main(dt, tend, K[0], M, boris_2nd_order) 

143 

144 # Ham_SDC2=main(dt, tend, K[1], M, boris_2nd_order) 

145 

146 # Ham_SDC3=main(dt, tend, K[2], M, boris_2nd_order) 

147 

148 # Ham_RKN=main(dt, tend, 1, M, RKN) 

149 

150 # plot_Hamiltonian_error(K, M, dt)