Coverage for pySDC/tutorial/step_4/D_MLSDC_with_particles.py: 100%

1import time 

2from pathlib import Path 

3 

4import numpy as np 

5 

6from pySDC.helpers.stats_helper import get_sorted 

7 

8from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

9from pySDC.implementations.problem_classes.PenningTrap_3D import penningtrap 

10from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order 

11from pySDC.implementations.transfer_classes.TransferParticles_NoCoarse import particles_to_particles 

12from pySDC.tutorial.step_3.HookClass_Particles import particle_hook 

13from pySDC.tutorial.step_4.PenningTrap_3D_coarse import penningtrap_coarse 

14 

15 

16def main(): 

17 """ 

18 A simple test program to compare SDC with two flavors of MLSDC for particle dynamics 

19 """ 

20 

21 # run SDC, MLSDC and MLSDC plus f-interpolation and compare 

22 stats_sdc, time_sdc = run_penning_trap_simulation(mlsdc=False) 

23 stats_mlsdc, time_mlsdc = run_penning_trap_simulation(mlsdc=True) 

24 stats_mlsdc_finter, time_mlsdc_finter = run_penning_trap_simulation(mlsdc=True, finter=True) 

25 

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

27 f = open('data/step_4_D_out.txt', 'w') 

28 out = 'Timings for SDC, MLSDC and MLSDC+finter: %12.8f -- %12.8f -- %12.8f' % ( 

29 time_sdc, 

30 time_mlsdc, 

31 time_mlsdc_finter, 

32 ) 

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

34 print(out) 

35 

36 # sort and convert stats to list, sorted by iteration numbers (only pre- and after-step are present here) 

37 energy_sdc = get_sorted(stats_sdc, type='etot', sortby='iter') 

38 energy_mlsdc = get_sorted(stats_mlsdc, type='etot', sortby='iter') 

39 energy_mlsdc_finter = get_sorted(stats_mlsdc_finter, type='etot', sortby='iter') 

40 

41 # get base energy and show differences 

42 base_energy = energy_sdc[0][1] 

43 for item in energy_sdc: 

44 out = 'Total energy and relative deviation in iteration %2i: %12.10f -- %12.8e' % ( 

45 item[0], 

46 item[1], 

47 abs(base_energy - item[1]) / base_energy, 

48 ) 

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

50 print(out) 

51 for item in energy_mlsdc: 

52 out = 'Total energy and relative deviation in iteration %2i: %12.10f -- %12.8e' % ( 

53 item[0], 

54 item[1], 

55 abs(base_energy - item[1]) / base_energy, 

56 ) 

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

58 print(out) 

59 for item in energy_mlsdc_finter: 

60 out = 'Total energy and relative deviation in iteration %2i: %12.10f -- %12.8e' % ( 

61 item[0], 

62 item[1], 

63 abs(base_energy - item[1]) / base_energy, 

64 ) 

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

66 print(out) 

67 f.close() 

68 

69 assert ( 

70 abs(energy_sdc[-1][1] - energy_mlsdc[-1][1]) / base_energy < 6e-10 

71 ), 'ERROR: energy deviated too much between SDC and MLSDC, got %s' % ( 

72 abs(energy_sdc[-1][1] - energy_mlsdc[-1][1]) / base_energy 

73 ) 

74 assert ( 

75 abs(energy_mlsdc[-1][1] - energy_mlsdc_finter[-1][1]) / base_energy < 8e-10 

76 ), 'ERROR: energy deviated too much after using finter, got %s' % ( 

77 abs(energy_mlsdc[-1][1] - energy_mlsdc_finter[-1][1]) / base_energy 

78 ) 

79 

80 

81def run_penning_trap_simulation(mlsdc, finter=False): 

82 """ 

83 A simple test program to run IMEX SDC for a single time step 

84 """ 

85 # initialize level parameters 

86 level_params = dict() 

87 level_params['restol'] = 1e-07 

88 level_params['dt'] = 1.0 / 8 

89 

90 # initialize sweeper parameters 

91 sweeper_params = dict() 

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

93 sweeper_params['num_nodes'] = 5 

94 

95 # initialize problem parameters for the Penning trap 

96 problem_params = dict() 

97 problem_params['omega_E'] = 4.9 # E-field frequency 

98 problem_params['omega_B'] = 25.0 # B-field frequency 

99 problem_params['u0'] = np.array([[10, 0, 0], [100, 0, 100], [1], [1]], dtype=object) # initial center of positions 

100 problem_params['nparts'] = 50 # number of particles in the trap 

101 problem_params['sig'] = 0.1 # smoothing parameter for the forces 

102 

103 # initialize step parameters 

104 step_params = dict() 

105 step_params['maxiter'] = 20 

106 

107 # initialize controller parameters 

108 controller_params = dict() 

109 controller_params['hook_class'] = particle_hook # specialized hook class for more statistics and output 

110 controller_params['logger_level'] = 30 

111 

112 transfer_params = dict() 

113 transfer_params['finter'] = finter 

114 

115 # Fill description dictionary for easy hierarchy creation 

116 description = dict() 

117 if mlsdc: 

118 # MLSDC: provide list of two problem classes: one for the fine, one for the coarse level 

119 description['problem_class'] = [penningtrap, penningtrap_coarse] 

120 else: 

121 # SDC: provide only one problem class 

122 description['problem_class'] = penningtrap 

123 description['problem_params'] = problem_params 

124 description['sweeper_class'] = boris_2nd_order 

125 description['sweeper_params'] = sweeper_params 

126 description['level_params'] = level_params 

127 description['step_params'] = step_params 

128 description['space_transfer_class'] = particles_to_particles 

129 description['base_transfer_params'] = transfer_params 

130 

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

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

133 

134 # set time parameters 

135 t0 = 0.0 

136 Tend = level_params['dt'] 

137 

138 # get initial values on finest level 

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

140 uinit = P.u_init() 

141 

142 # call and time main function to get things done... 

143 start_time = time.perf_counter() 

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

145 end_time = time.perf_counter() - start_time 

146 

147 return stats, end_time 

148 

149 

150if __name__ == "__main__": 

151 main()