## 55 statements

, created at 2024-09-20 17:10 +0000

1import numpy as np

2from pathlib import Path

4from pySDC.helpers.stats_helper import get_sorted

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

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)

18 err, stats = run_penning_trap_simulation()

20 # filter statistics type (etot)

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

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()

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

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

51 # initialize sweeper parameters

52 sweeper_params = dict()

54 sweeper_params['num_nodes'] = 3

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

64 # initialize step parameters

65 step_params = dict()

66 step_params['maxiter'] = 20

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'

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

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

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

86 # set time parameters

87 t0 = 0.0

88 Tend = level_params['dt']

90 # get initial values on finest level

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

92 uinit = P.u_init()

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

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

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)

101 return err, stats

104if __name__ == "__main__":

105 main()