Coverage for pySDC/projects/Hamiltonian/fput.py: 100%
134 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +0000
1import os
2from collections import defaultdict
4import dill
5import numpy as np
7import pySDC.helpers.plot_helper as plt_helper
8from pySDC.helpers.stats_helper import get_sorted, filter_stats
10from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
11from pySDC.implementations.problem_classes.FermiPastaUlamTsingou import fermi_pasta_ulam_tsingou
12from pySDC.implementations.sweeper_classes.verlet import verlet
13from pySDC.implementations.transfer_classes.TransferParticles_NoCoarse import particles_to_particles
14from pySDC.projects.Hamiltonian.hamiltonian_and_energy_output import hamiltonian_and_energy_output
17def setup_fput():
18 """
19 Helper routine for setting up everything for the Fermi-Pasta-Ulam-Tsingou problem
21 Returns:
22 description (dict): description of the controller
23 controller_params (dict): controller parameters
24 """
26 # initialize level parameters
27 level_params = dict()
28 level_params['restol'] = 1e-12
29 level_params['dt'] = 2.0
31 # initialize sweeper parameters
32 sweeper_params = dict()
33 sweeper_params['quad_type'] = 'LOBATTO'
34 sweeper_params['num_nodes'] = [5, 3]
35 sweeper_params['initial_guess'] = 'zero'
37 # initialize problem parameters for the Penning trap
38 problem_params = dict()
39 problem_params['npart'] = 2048
40 problem_params['alpha'] = 0.25
41 problem_params['k'] = 1.0
42 problem_params['energy_modes'] = [[1, 2, 3, 4]]
44 # initialize step parameters
45 step_params = dict()
46 step_params['maxiter'] = 50
48 # initialize controller parameters
49 controller_params = dict()
50 controller_params['hook_class'] = hamiltonian_and_energy_output
51 controller_params['logger_level'] = 30
53 # Fill description dictionary for easy hierarchy creation
54 description = dict()
55 description['problem_class'] = fermi_pasta_ulam_tsingou
56 description['problem_params'] = problem_params
57 description['sweeper_class'] = verlet
58 description['sweeper_params'] = sweeper_params
59 description['level_params'] = level_params
60 description['step_params'] = step_params
61 description['space_transfer_class'] = particles_to_particles
63 return description, controller_params
66def run_simulation():
67 """
68 Routine to run the simulation of a second order problem
70 """
72 description, controller_params = setup_fput()
73 # set time parameters
74 t0 = 0.0
75 # set this to 10000 to reproduce the picture in
76 # http://www.scholarpedia.org/article/Fermi-Pasta-Ulam_nonlinear_lattice_oscillations
77 Tend = 1000.0
78 num_procs = 1
80 f = open('data/fput_out.txt', 'w')
81 out = 'Running fput problem with %s processors...' % num_procs
82 f.write(out + '\n')
83 print(out)
85 # instantiate the controller
86 controller = controller_nonMPI(num_procs=num_procs, controller_params=controller_params, description=description)
88 # get initial values on finest level
89 P = controller.MS[0].levels[0].prob
90 uinit = P.u_exact(t=t0)
92 # call main function to get things done...
93 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
95 # filter statistics by type (number of iterations)
96 iter_counts = get_sorted(stats, type='niter', sortby='time')
98 # compute and print statistics
99 # for item in iter_counts:
100 # out = 'Number of iterations for time %4.2f: %2i' % item
101 # f.write(out + '\n')
102 # print(out)
104 niters = np.array([item[1] for item in iter_counts])
105 out = ' Mean number of iterations: %4.2f' % np.mean(niters)
106 f.write(out + '\n')
107 print(out)
108 out = ' Range of values for number of iterations: %2i ' % np.ptp(niters)
109 f.write(out + '\n')
110 print(out)
111 out = ' Position of max/min number of iterations: %2i -- %2i' % (int(np.argmax(niters)), int(np.argmin(niters)))
112 f.write(out + '\n')
113 print(out)
114 out = ' Std and var for number of iterations: %4.2f -- %4.2f' % (float(np.std(niters)), float(np.var(niters)))
115 f.write(out + '\n')
116 print(out)
118 # get runtime
119 timing_run = get_sorted(stats, type='timing_run', sortby='time')[0][1]
120 out = '... took %6.4f seconds to run this.' % timing_run
121 f.write(out + '\n')
122 print(out)
123 f.close()
125 # assert np.mean(niters) <= 3.46, 'Mean number of iterations is too high, got %s' % np.mean(niters)
127 fname = 'data/fput.dat'
128 f = open(fname, 'wb')
129 dill.dump(stats, f)
130 f.close()
132 assert os.path.isfile(fname), 'Run for %s did not create stats file'
135def show_results(cwd=''):
136 """
137 Helper function to plot the error of the Hamiltonian
139 Args:
140 cwd (str): current working directory
141 """
143 # read in the dill data
144 f = open(cwd + 'data/fput.dat', 'rb')
145 stats = dill.load(f)
146 f.close()
148 plt_helper.mpl.style.use('classic')
149 plt_helper.setup_mpl()
151 # HAMILTONIAN PLOTTING #
152 # extract error in hamiltonian and prepare for plotting
153 extract_stats = filter_stats(stats, type='err_hamiltonian')
154 result = defaultdict(list)
155 for k, v in extract_stats.items():
156 result[k.iter].append((k.time, v))
157 for k, _ in result.items():
158 result[k] = sorted(result[k], key=lambda x: x[0])
160 plt_helper.newfig(textwidth=238.96, scale=0.89)
162 # Rearrange data for easy plotting
163 err_ham = 1
164 for k, v in result.items():
165 time = [item[0] for item in v]
166 ham = [item[1] for item in v]
167 err_ham = ham[-1]
168 plt_helper.plt.semilogy(time, ham, '-', lw=1, label='Iter ' + str(k))
169 print(err_ham)
170 # assert err_ham < 6E-10, 'Error in the Hamiltonian is too large, got %s' % err_ham
172 plt_helper.plt.xlabel('Time')
173 plt_helper.plt.ylabel('Error in Hamiltonian')
174 plt_helper.plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
176 fname = 'data/fput_hamiltonian'
177 plt_helper.savefig(fname)
179 assert os.path.isfile(fname + '.pdf'), 'ERROR: plotting did not create PDF file'
180 # assert os.path.isfile(fname + '.pgf'), 'ERROR: plotting did not create PGF file'
181 assert os.path.isfile(fname + '.png'), 'ERROR: plotting did not create PNG file'
183 # ENERGY PLOTTING #
184 # extract error in hamiltonian and prepare for plotting
185 result = get_sorted(stats, type='energy_step', sortby='time')
187 plt_helper.newfig(textwidth=238.96, scale=0.89)
189 # Rearrange data for easy plotting
190 for mode in result[0][1].keys():
191 time = [item[0] for item in result]
192 energy = [item[1][mode] for item in result]
193 plt_helper.plt.plot(time, energy, label=str(mode) + 'th mode')
195 plt_helper.plt.xlabel('Time')
196 plt_helper.plt.ylabel('Energy')
197 plt_helper.plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
199 fname = 'data/fput_energy'
200 plt_helper.savefig(fname)
202 assert os.path.isfile(fname + '.pdf'), 'ERROR: plotting did not create PDF file'
203 # assert os.path.isfile(fname + '.pgf'), 'ERROR: plotting did not create PGF file'
204 assert os.path.isfile(fname + '.png'), 'ERROR: plotting did not create PNG file'
206 # POSITION PLOTTING #
207 # extract positions and prepare for plotting
208 result = get_sorted(stats, type='position', sortby='time')
210 plt_helper.newfig(textwidth=238.96, scale=0.89)
212 # Rearrange data for easy plotting
213 nparts = len(result[0][1])
214 nsteps = len(result)
215 pos = np.zeros((nparts, nsteps))
216 time = np.zeros(nsteps)
217 for idx, item in enumerate(result):
218 time[idx] = item[0]
219 for n in range(nparts):
220 pos[n, idx] = item[1][n]
222 for n in range(min(nparts, 16)):
223 plt_helper.plt.plot(time, pos[n, :])
225 fname = 'data/fput_positions'
226 plt_helper.savefig(fname)
228 assert os.path.isfile(fname + '.pdf'), 'ERROR: plotting did not create PDF file'
229 # assert os.path.isfile(fname + '.pgf'), 'ERROR: plotting did not create PGF file'
230 assert os.path.isfile(fname + '.png'), 'ERROR: plotting did not create PNG file'
233def main():
234 run_simulation()
235 show_results()
238if __name__ == "__main__":
239 main()