Coverage for pySDC/projects/Second_orderSDC/penningtrap_HookClass.py: 100%

1from pySDC.core.hooks import Hooks 

2import numpy as np 

3 

4 

5class particles_output(Hooks): 

6 def __init__(self): 

7 """ 

8 Initialization of particles output 

9 """ 

10 super(particles_output, self).__init__() 

11 

12 def pre_run(self, step, level_number): 

13 """ 

14 Overwrite default routine called before time-loop starts 

15 Args: 

16 step: the current step 

17 level_number: the current level number 

18 """ 

19 super(particles_output, self).pre_run(step, level_number) 

20 

21 def post_step(self, step, level_number): 

22 """ 

23 Default routine called after each iteration 

24 Args: 

25 step: the current step 

26 level_number: the current level number 

27 """ 

28 

29 super(particles_output, self).post_step(step, level_number) 

30 

31 # some abbreviations 

32 L = step.levels[level_number] 

33 

34 # self.bar_run.update(L.time) 

35 

36 L.sweep.compute_end_point() 

37 part = L.uend 

38 N = L.prob.nparts 

39 

40 # ============================================================================= 

41 # 

42 

43 try: # pragma: no cover 

44 _unused = L.prob.Harmonic_oscillator 

45 

46 # add up kinetic and potntial contributions to total energy 

47 epot = 0 

48 ekin = 0 

49 name = str(L.sweep) 

50 for n in range(N): 

51 epot += 1 / 2.0 * np.dot(part.pos[:, n], part.pos[:, n]) 

52 ekin += 1 / 2.0 * np.dot(part.vel[:, n], part.vel[:, n]) 

53 Energy = epot + ekin 

54 uinit = L.u[0] 

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

56 Ham = abs(Energy - H0) / abs(H0) 

57 if 'RKN' in name: 

58 filename = 'data/Ham_RKN' + '{}.txt'.format(step.params.maxiter) 

59 else: 

60 filename = 'data/Ham_SDC' + '{}{}.txt'.format(L.sweep.coll.num_nodes, step.params.maxiter) 

61 

62 if L.time == 0.0: 

63 file = open(filename, 'w') 

64 file.write(str('time') + " | " + str('Hamiltonian error') + '\n') 

65 else: 

66 file = open(filename, 'a') 

67 

68 file.write(str(L.time) + " | " + str(Ham) + '\n') 

69 

70 file.close() 

71 except AttributeError: 

72 pass 

73 # ============================================================================= 

74 

75 part_exact = L.prob.u_exact(L.time + L.dt) 

76 self.add_to_stats( 

77 process=step.status.slot, 

78 time=L.time, 

79 level=L.level_index, 

80 iter=step.status.iter, 

81 sweep=L.status.sweep, 

82 type='position', 

83 value=part.pos, 

84 ) 

85 self.add_to_stats( 

86 process=step.status.slot, 

87 time=L.time, 

88 level=L.level_index, 

89 iter=step.status.iter, 

90 sweep=L.status.sweep, 

91 type='velocity', 

92 value=part.vel, 

93 ) 

94 self.add_to_stats( 

95 process=step.status.slot, 

96 time=L.time, 

97 level=L.level_index, 

98 iter=step.status.iter, 

99 sweep=L.status.sweep, 

100 type='position_exact', 

101 value=part_exact.pos, 

102 ) 

103 self.add_to_stats( 

104 process=step.status.slot, 

105 time=L.time, 

106 level=L.level_index, 

107 iter=step.status.iter, 

108 sweep=L.status.sweep, 

109 type='velocity_exact', 

110 value=part_exact.vel, 

111 )