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

1import numpy as np 

2from pathlib import Path 

3 

4from pySDC.helpers.stats_helper import get_sorted 

5 

6from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI 

7from pySDC.implementations.problem_classes.PenningTrap_3D import penningtrap 

8from pySDC.implementations.sweeper_classes.boris_2nd_order import boris_2nd_order 

9from pySDC.tutorial.step_3.HookClass_Particles import particle_hook 

10 

11 

12def main(): 

13 """ 

14 A simple test program to retrieve user-defined statistics from a run 

15 """ 

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

17 

18 err, stats = run_penning_trap_simulation() 

19 

20 # filter statistics type (etot) 

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

22 

23 # get base energy and show difference 

24 base_energy = energy[0][1] 

25 f = open('data/step_3_B_out.txt', 'a') 

26 for item in energy: 

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

28 item[0], 

29 item[1], 

30 abs(base_energy - item[1]), 

31 ) 

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

33 print(out) 

34 f.close() 

35 

36 assert abs(base_energy - energy[-1][1]) < 15, 'ERROR: energy deviated too much, got %s' % ( 

37 base_energy - energy[-1][1] 

38 ) 

39 assert err < 5e-04, "ERROR: solution is not as exact as expected, got %s" % err 

40 

41 

42def run_penning_trap_simulation(): 

43 """ 

44 A simple test program to run IMEX SDC for a single time step of the penning trap example 

45 """ 

46 # initialize level parameters 

47 level_params = dict() 

48 level_params['restol'] = 1e-08 

49 level_params['dt'] = 1.0 / 16 

50 

51 # initialize sweeper parameters 

52 sweeper_params = dict() 

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

54 sweeper_params['num_nodes'] = 3 

55 

56 # initialize problem parameters for the Penning trap 

57 problem_params = dict() 

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

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

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

61 problem_params['nparts'] = 1 # number of particles in the trap 

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

63 

64 # initialize step parameters 

65 step_params = dict() 

66 step_params['maxiter'] = 20 

67 

68 # initialize controller parameters 

69 controller_params = dict() 

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

71 controller_params['log_to_file'] = True 

72 controller_params['fname'] = 'data/step_3_B_out.txt' 

73 

74 # Fill description dictionary for easy hierarchy creation 

75 description = dict() 

76 description['problem_class'] = penningtrap 

77 description['problem_params'] = problem_params 

78 description['sweeper_class'] = boris_2nd_order 

79 description['sweeper_params'] = sweeper_params 

80 description['level_params'] = level_params 

81 description['step_params'] = step_params 

82 

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

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

85 

86 # set time parameters 

87 t0 = 0.0 

88 Tend = level_params['dt'] 

89 

90 # get initial values on finest level 

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

92 uinit = P.u_init() 

93 

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

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

96 

97 # compute error compared to know exact solution for one particle 

98 uex = P.u_exact(Tend) 

99 err = np.linalg.norm(uex.pos - uend.pos, np.inf) / np.linalg.norm(uex.pos, np.inf) 

100 

101 return err, stats 

102 

103 

104if __name__ == "__main__": 

105 main()