Coverage for pySDC/tutorial/step_4/D_MLSDC_with_particles.py: 100%
78 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
1import time
2from pathlib import Path
4import numpy as np
6from pySDC.helpers.stats_helper import get_sorted
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
16def main():
17 """
18 A simple test program to compare SDC with two flavors of MLSDC for particle dynamics
19 """
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)
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)
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')
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()
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 )
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
90 # initialize sweeper parameters
91 sweeper_params = dict()
92 sweeper_params['quad_type'] = 'RADAU-RIGHT'
93 sweeper_params['num_nodes'] = 5
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
103 # initialize step parameters
104 step_params = dict()
105 step_params['maxiter'] = 20
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
112 transfer_params = dict()
113 transfer_params['finter'] = finter
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
131 # instantiate the controller (no controller parameters used here)
132 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
134 # set time parameters
135 t0 = 0.0
136 Tend = level_params['dt']
138 # get initial values on finest level
139 P = controller.MS[0].levels[0].prob
140 uinit = P.u_init()
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
147 return stats, end_time
150if __name__ == "__main__":
151 main()