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

## 63 statements

, created at 2024-09-20 16:55 +0000

1from pathlib import Path

2import numpy as np

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

11def main():

12 """

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

14 """

15 stats_dict = run_simulation()

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

32 # set expected differences and check

33 ediff_expect = dict()

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]

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

50 # initialize sweeper parameters

51 sweeper_params = dict()

52 sweeper_params['num_nodes'] = 3

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

62 # initialize step parameters

63 step_params = dict()

64 step_params['maxiter'] = 20

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

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

79 # assemble and loop over list of collocation classes

81 stats_dict = dict()

84 description['sweeper_params'] = sweeper_params

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

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

89 # set time parameters

90 t0 = 0.0

91 Tend = level_params['dt']

93 # get initial values on finest level

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

95 uinit = P.u_init()

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

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

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

101 stats_dict[qtype] = stats

103 return stats_dict

106if __name__ == "__main__":

107 main()