Coverage for pySDC/tutorial/step_3/C_study_collocations.py: 100%

1from pathlib import Path 

2import numpy as np 

3 

4from pySDC.helpers.stats_helper import get_sorted 

5from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

6from pySDC.implementations.problem_classes.PenningTrap_3D import penningtrap 

7from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order 

8from pySDC.tutorial.step_3.HookClass_Particles import particle_hook 

9 

10 

11def main(): 

12 """ 

13 A simple test program to show the energy deviation for different quadrature nodes 

14 """ 

15 stats_dict = run_simulation() 

16 

17 ediff = dict() 

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

19 f = open('data/step_3_C_out.txt', 'w') 

20 for cclass, stats in stats_dict.items(): 

21 # filter and convert/sort statistics by etot and iterations 

22 energy = get_sorted(stats, type='etot', sortby='iter') 

23 # compare base and final energy 

24 base_energy = energy[0][1] 

25 final_energy = energy[-1][1] 

26 ediff[cclass] = abs(base_energy - final_energy) 

27 out = "Energy deviation for %s: %12.8e" % (cclass, ediff[cclass]) 

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

29 print(out) 

30 f.close() 

31 

32 # set expected differences and check 

33 ediff_expect = dict() 

34 ediff_expect['RADAU-RIGHT'] = 15 

35 ediff_expect['LOBATTO'] = 1e-05 

36 ediff_expect['GAUSS'] = 3e-05 

37 for k, v in ediff.items(): 

38 assert v < ediff_expect[k], "ERROR: energy deviated too much, got %s" % ediff[k] 

39 

40 

41def run_simulation(): 

42 """ 

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

44 """ 

45 # initialize level parameters 

46 level_params = dict() 

47 level_params['restol'] = 1e-06 

48 level_params['dt'] = 1.0 / 16 

49 

50 # initialize sweeper parameters 

51 sweeper_params = dict() 

52 sweeper_params['num_nodes'] = 3 

53 

54 # initialize problem parameters 

55 problem_params = dict() 

56 problem_params['omega_E'] = 4.9 

57 problem_params['omega_B'] = 25.0 

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

59 problem_params['nparts'] = 1 

60 problem_params['sig'] = 0.1 

61 

62 # initialize step parameters 

63 step_params = dict() 

64 step_params['maxiter'] = 20 

65 

66 # initialize controller parameters 

67 controller_params = dict() 

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

69 controller_params['logger_level'] = 30 # reduce verbosity of each run 

70 

71 # Fill description dictionary for easy hierarchy creation 

72 description = dict() 

73 description['problem_class'] = penningtrap 

74 description['problem_params'] = problem_params 

75 description['sweeper_class'] = boris_2nd_order 

76 description['level_params'] = level_params 

77 description['step_params'] = step_params 

78 

79 # assemble and loop over list of collocation classes 

80 quad_types = ['RADAU-RIGHT', 'GAUSS', 'LOBATTO'] 

81 stats_dict = dict() 

82 for qtype in quad_types: 

83 sweeper_params['quad_type'] = qtype 

84 description['sweeper_params'] = sweeper_params 

85 

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

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

88 

89 # set time parameters 

90 t0 = 0.0 

91 Tend = level_params['dt'] 

92 

93 # get initial values on finest level 

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

95 uinit = P.u_init() 

96 

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

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

99 

100 # gather stats in dictionary, collocation classes being the keys 

101 stats_dict[qtype] = stats 

102 

103 return stats_dict 

104 

105 

106if __name__ == "__main__": 

107 main()