Skip to content

Commit 8114392

Browse files
giadarolgiadarol
committed
results seem off by a factor 2
investigating
1 parent 231875a commit 8114392

File tree

5 files changed

+420
-5
lines changed

5 files changed

+420
-5
lines changed

‎test_physics/001_check_output.py‎

Lines changed: 99 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,99 @@
1+
import numpy as np
2+
import pylab as pl
3+
import myfilemanager as mfm
4+
import mystyle as ms
5+
6+
import time
7+
8+
n_segments = 3
9+
n_part_per_turn = 5000
10+
N_turns = 8
11+
12+
filename = '../../PyECLOUD/testing/tests_PyEC4PyHT/headtail_for_test/test_protons/SPS_Q20_proton_check_dipole_3kicks_20150212_prb.dat'
13+
appo = np.loadtxt(filename)
14+
15+
parid = np.reshape(appo[:,0], (-1, n_part_per_turn))[::n_segments,:]
16+
x = np.reshape(appo[:,1], (-1, n_part_per_turn))[::n_segments,:]
17+
xp = np.reshape(appo[:,2], (-1, n_part_per_turn))[::n_segments,:]
18+
y = np.reshape(appo[:,3], (-1, n_part_per_turn))[::n_segments,:]
19+
yp =np.reshape(appo[:,4], (-1, n_part_per_turn))[::n_segments,:]
20+
z = np.reshape(appo[:,5], (-1, n_part_per_turn))[::n_segments,:]
21+
zp = np.reshape(appo[:,6], (-1, n_part_per_turn))[::n_segments,:]
22+
23+
rms_err_x_list = []
24+
rms_err_y_list = []
25+
for i_turn in xrange(N_turns-1):
26+
27+
ob = mfm.object_with_arrays_and_scalar_from_h5('particles_at_turn_%d.h5'%i_turn)
28+
29+
# sort id and momenta after track
30+
31+
id_after = ob.id_after
32+
xp_after = ob.xp_after
33+
yp_after = ob.yp_after
34+
z_after = ob.z_after
35+
id_before = ob.id_before
36+
xp_before = ob.xp_before
37+
yp_before = ob.yp_before
38+
39+
fig = pl.figure(100, figsize=(18,10))
40+
pl.clf()
41+
sp1 = pl.subplot(2,3,1)
42+
pl.plot(z_after, (100*np.abs((xp_after-xp_before)-(xp[i_turn+1, :]-xp[0, :]))/np.std(xp[i_turn+1, :]-xp[0, :])), '.r', markersize = 2)
43+
pl.grid('on')
44+
pl.ylabel('''error($\Delta$x')/$\Delta$x' [%]''')
45+
pl.subplot(2,3,4, sharex=sp1)
46+
pl.plot(z_after, (xp[i_turn+1, :]-xp[0, :]), '.r', label='HT', markersize = 2)
47+
pl.plot(z_after, (xp_after-xp_before), '.b', label='PyHT', markersize = 2)
48+
ms.sciy()
49+
pl.grid('on')
50+
pl.ylabel("$\Delta$x'")
51+
pl.xlabel('z[m]')
52+
pl.legend(prop = {'size':14})
53+
pl.subplot(2,3,3)
54+
pl.hist((100*np.abs((xp_after-xp_before)-(xp[i_turn+1, :]-xp[0, :]))/np.std(xp[i_turn+1, :]-xp[0, :])), bins=100, range=(0,100))
55+
pl.xlabel('Error_x [%]')
56+
pl.ylabel('Occurrences')
57+
pl.xlim(0,100)
58+
rms_err_x = np.std(100*np.abs((xp_after-xp_before)-(xp[i_turn+1, :]-xp[0, :]))/np.std(xp[i_turn+1, :]-xp[0, :]))
59+
60+
sp1 = pl.subplot(2,3,2)
61+
pl.plot(z_after, (100*np.abs((yp_after-yp_before)-(yp[i_turn+1, :]-yp[0, :]))/np.std(yp[i_turn+1, :]-yp[0, :])), '.r', markersize = 2)
62+
pl.grid('on')
63+
pl.ylabel('''error($\Delta$y')/$\Delta$y' [%]''')
64+
pl.subplot(2,3,5, sharex=sp1)
65+
pl.plot(z_after, (yp[i_turn+1, :]-yp[0, :]), '.r', label='HT', markersize = 2)
66+
pl.plot(z_after, (yp_after-yp_before), '.b', label='PyHT', markersize = 2)
67+
ms.sciy()
68+
pl.grid('on')
69+
pl.ylabel("$\Delta$y'")
70+
pl.xlabel('z[m]')
71+
pl.legend(prop = {'size':14})
72+
pl.subplot(2,3,6)
73+
pl.hist((100*np.abs((yp_after-yp_before)-(yp[i_turn+1, :]-yp[0, :]))/np.std(yp[i_turn+1, :]-yp[0, :])), bins=100, range=(0,100))
74+
pl.xlim(0,100)
75+
pl.xlabel('Error_y [%]')
76+
pl.ylabel('Occurrences')
77+
rms_err_y = np.std(100*np.abs((yp_after-yp_before)-(yp[i_turn+1, :]-yp[0, :]))/np.std(yp[i_turn+1, :]-yp[0, :]))
78+
79+
80+
pl.suptitle('Turn %d rms_err_x = %e rms_err_y = %e'%(i_turn, rms_err_x, rms_err_y))
81+
82+
pl.savefig(filename.split('_prb.dat')[0]+'_%02d.png'%i_turn, dpi=150)
83+
84+
rms_err_x_list.append(rms_err_x)
85+
rms_err_y_list.append(rms_err_y)
86+
87+
pl.ion()
88+
pl.draw()
89+
time.sleep(1.)
90+
91+
pl.figure(1000)
92+
pl.plot(rms_err_x_list, '.-', markersize = 10, linewidth=2, label='x')
93+
pl.plot(rms_err_y_list, '.-', markersize = 10, linewidth=2, label='y')
94+
pl.grid('on')
95+
pl.ylabel('''Relative r.m.s. error [%]''')
96+
pl.xlabel('Turn')
97+
pl.legend()
98+
pl.savefig(filename.split('_prb.dat')[0]+'_errors.png', dpi=200)
99+
pl.show()

‎test_physics/Simulation_with_eclouds.py‎

Lines changed: 14 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -7,21 +7,22 @@
77
epsn_x = 2.5e-6
88
epsn_y = 2.5e-6
99
B_multip = [0.5]
10+
n_slices = 64
11+
n_segments = 3
1012

1113
filename = '../../PyECLOUD/testing/tests_PyEC4PyHT/headtail_for_test/test_protons/SPS_Q20_proton_check_dipole_3kicks_20150212_prb.dat'
1214
B_multip = [0.5]
15+
N_turns_to_simulate = 8
1316

1417

1518
class Simulation(object):
1619
def __init__(self):
17-
self.N_turns = None
20+
self.N_turns = N_turns_to_simulate
1821

1922
def init_all(self):
2023

21-
n_slices = 64
24+
2225
self.n_slices = n_slices
23-
24-
n_segments = 3
2526
self.n_segments = n_segments
2627

2728
from machines_for_testing import SPS
@@ -107,7 +108,6 @@ def init_master(self):
107108
yp =np.reshape(appo[:,4], (-1, self.n_part_per_turn))[::self.n_segments,:]
108109
z = np.reshape(appo[:,5], (-1, self.n_part_per_turn))[::self.n_segments,:]
109110
zp = np.reshape(appo[:,6], (-1, self.n_part_per_turn))[::self.n_segments,:]
110-
self.N_turns = len(x[:,0])
111111

112112
# replace first particles with HEADTAIL ones
113113
bunch.x[:self.n_part_per_turn] = x[0,:]
@@ -133,6 +133,8 @@ def init_master(self):
133133
slice_obj_list = bunch.extract_slices(self.slicer)
134134

135135
pieces_to_be_treated = slice_obj_list
136+
137+
print 'N_turns', self.N_turns
136138

137139
return pieces_to_be_treated
138140

@@ -168,6 +170,13 @@ def finalize_turn_on_master(self, pieces_treated):
168170
# save results
169171
import myfilemanager as mfm
170172
mfm.save_dict_to_h5('particles_at_turn_%d.h5'%self.ring_of_CPUs.i_turn,{\
173+
'id_after': id_after,
174+
'xp_after': xp_after,
175+
'yp_after': yp_after,
176+
'z_after': z_after,
177+
'id_before':self.id_before,
178+
'xp_before':self.xp_before,
179+
'yp_before':self.yp_before})
171180

172181

173182
# prepare next turn (re-slice)

‎test_physics/a001_compare_kick_against_headtail_serial.py‎

Lines changed: 204 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,204 @@
1+
import sys, os
2+
BIN=os.path.expanduser('../../')
3+
sys.path.append(BIN)
4+
5+
6+
import numpy as np
7+
import pylab as pl
8+
import mystyle as ms
9+
import time
10+
11+
show_movie = False
12+
13+
14+
filename = '../../PyECLOUD/testing/tests_PyEC4PyHT/headtail_for_test/test_protons/SPS_Q20_proton_check_dipole_3kicks_20150212_prb.dat'
15+
B_multip = [0.5]
16+
N_kicks = 3
17+
18+
n_part_per_turn = 5000
19+
20+
appo = np.loadtxt(filename)
21+
22+
parid = np.reshape(appo[:,0], (-1, n_part_per_turn))[::N_kicks,:]
23+
x = np.reshape(appo[:,1], (-1, n_part_per_turn))[::N_kicks,:]
24+
xp = np.reshape(appo[:,2], (-1, n_part_per_turn))[::N_kicks,:]
25+
y = np.reshape(appo[:,3], (-1, n_part_per_turn))[::N_kicks,:]
26+
yp =np.reshape(appo[:,4], (-1, n_part_per_turn))[::N_kicks,:]
27+
z = np.reshape(appo[:,5], (-1, n_part_per_turn))[::N_kicks,:]
28+
zp = np.reshape(appo[:,6], (-1, n_part_per_turn))[::N_kicks,:]
29+
N_turns = len(x[:,0])
30+
31+
pl.close('all')
32+
ms.mystyle(fontsz=14)
33+
34+
35+
# define machine for PyHEADTAIL
36+
from PyHEADTAIL.particles.slicing import UniformBinSlicer
37+
from machines_for_testing import SPS
38+
machine = SPS(n_segments = N_kicks, machine_configuration = 'Q20-injection', accQ_x=20., accQ_y=20.)
39+
#machine.one_turn_map.remove(machine.longitudinal_map)
40+
41+
42+
# compute sigma x and y
43+
epsn_x = 2.5e-6
44+
epsn_y = 2.5e-6
45+
46+
inj_optics = machine.transverse_map.get_injection_optics()
47+
sigma_x = np.sqrt(inj_optics['beta_x']*epsn_x/machine.betagamma)
48+
sigma_y = np.sqrt(inj_optics['beta_y']*epsn_y/machine.betagamma)
49+
50+
# define apertures and Dh_sc to simulate headtail conditions
51+
x_aper = 20*sigma_x
52+
y_aper = 20*sigma_y
53+
Dh_sc = 2*x_aper/128/2
54+
if show_movie: Dh_sc*=2
55+
56+
# ecloud
57+
import PyECLOUD.PyEC4PyHT as PyEC4PyHT
58+
from PyHEADTAIL.particles.slicing import UniformBinSlicer
59+
slicer = UniformBinSlicer(n_slices = 64, n_sigma_z = 3.)
60+
61+
init_unif_edens_flag=1
62+
init_unif_edens=2e11
63+
N_MP_ele_init = 100000
64+
N_mp_max = N_MP_ele_init*4.
65+
66+
nel_mp_ref_0 = init_unif_edens*4*x_aper*y_aper/N_MP_ele_init
67+
68+
69+
ecloud = PyEC4PyHT.Ecloud(L_ecloud=machine.circumference/N_kicks, slicer=slicer,
70+
Dt_ref=25e-12, pyecl_input_folder='../../PyECLOUD/testing/tests_PyEC4PyHT/drift_sim/',
71+
x_aper=x_aper, y_aper=y_aper, Dh_sc=Dh_sc,
72+
init_unif_edens_flag=init_unif_edens_flag,
73+
init_unif_edens=init_unif_edens,
74+
N_MP_ele_init=N_MP_ele_init,
75+
N_mp_max=N_mp_max,
76+
nel_mp_ref_0=nel_mp_ref_0,
77+
B_multip=B_multip)
78+
machine.install_after_each_transverse_segment(ecloud)
79+
80+
if show_movie:
81+
ecloud.save_ele_distributions_last_track = True
82+
ecloud.save_ele_potential_and_field = True
83+
84+
# generate a bunch
85+
bunch = machine.generate_6D_Gaussian_bunch(n_macroparticles=300000, intensity=1.15e11, epsn_x=epsn_x, epsn_y=epsn_y, sigma_z=0.2)
86+
87+
88+
# replace first particles with HEADTAIL ones
89+
bunch.x[:n_part_per_turn] = x[0,:]
90+
bunch.xp[:n_part_per_turn] = xp[0,:]
91+
bunch.y[:n_part_per_turn] = y[0,:]
92+
bunch.yp[:n_part_per_turn] = yp[0,:]
93+
bunch.z[:n_part_per_turn] = z[0,:]
94+
bunch.dp[:n_part_per_turn] =zp[0,:]
95+
96+
# save id and momenta before track
97+
id_before = bunch.id[bunch.id<=n_part_per_turn]
98+
xp_before = bunch.xp[bunch.id<=n_part_per_turn]
99+
yp_before = bunch.yp[bunch.id<=n_part_per_turn]
100+
101+
rms_err_x_list = []
102+
rms_err_y_list = []
103+
for ii in xrange(N_turns-1):
104+
# track
105+
machine.track(bunch)#, verbose = True)
106+
print 'Turn', ii
107+
108+
# id and momenta after track
109+
id_after = bunch.id[bunch.id<=n_part_per_turn]
110+
xp_after = bunch.xp[bunch.id<=n_part_per_turn]
111+
z_after = bunch.z[bunch.id<=n_part_per_turn]
112+
yp_after = bunch.yp[bunch.id<=n_part_per_turn]
113+
114+
# sort id and momenta after track
115+
indsort = np.argsort(id_after)
116+
id_after = np.take(id_after, indsort)
117+
xp_after = np.take(xp_after, indsort)
118+
yp_after = np.take(yp_after, indsort)
119+
z_after = np.take(z_after, indsort)
120+
121+
fig = pl.figure(100, figsize=(18,10))
122+
pl.clf()
123+
sp1 = pl.subplot(2,3,1)
124+
pl.plot(z_after, (100*np.abs((xp_after-xp_before)-(xp[ii+1, :]-xp[0, :]))/np.std(xp[ii+1, :]-xp[0, :])), '.r', markersize = 2)
125+
pl.grid('on')
126+
pl.ylabel('''error($\Delta$x')/$\Delta$x' [%]''')
127+
pl.subplot(2,3,4, sharex=sp1)
128+
pl.plot(z_after, (xp[ii+1, :]-xp[0, :]), '.r', label='HT', markersize = 2)
129+
pl.plot(z_after, (xp_after-xp_before), '.b', label='PyHT', markersize = 2)
130+
ms.sciy()
131+
pl.grid('on')
132+
pl.ylabel("$\Delta$x'")
133+
pl.xlabel('z[m]')
134+
pl.legend(prop = {'size':14})
135+
pl.subplot(2,3,3)
136+
pl.hist((100*np.abs((xp_after-xp_before)-(xp[ii+1, :]-xp[0, :]))/np.std(xp[ii+1, :]-xp[0, :])), bins=100, range=(0,100))
137+
pl.xlabel('Error_x [%]')
138+
pl.ylabel('Occurrences')
139+
pl.xlim(0,100)
140+
rms_err_x = np.std(100*np.abs((xp_after-xp_before)-(xp[ii+1, :]-xp[0, :]))/np.std(xp[ii+1, :]-xp[0, :]))
141+
142+
sp1 = pl.subplot(2,3,2)
143+
pl.plot(z_after, (100*np.abs((yp_after-yp_before)-(yp[ii+1, :]-yp[0, :]))/np.std(yp[ii+1, :]-yp[0, :])), '.r', markersize = 2)
144+
pl.grid('on')
145+
pl.ylabel('''error($\Delta$y')/$\Delta$y' [%]''')
146+
pl.subplot(2,3,5, sharex=sp1)
147+
pl.plot(z_after, (yp[ii+1, :]-yp[0, :]), '.r', label='HT', markersize = 2)
148+
pl.plot(z_after, (yp_after-yp_before), '.b', label='PyHT', markersize = 2)
149+
ms.sciy()
150+
pl.grid('on')
151+
pl.ylabel("$\Delta$y'")
152+
pl.xlabel('z[m]')
153+
pl.legend(prop = {'size':14})
154+
pl.subplot(2,3,6)
155+
pl.hist((100*np.abs((yp_after-yp_before)-(yp[ii+1, :]-yp[0, :]))/np.std(yp[ii+1, :]-yp[0, :])), bins=100, range=(0,100))
156+
pl.xlim(0,100)
157+
pl.xlabel('Error_y [%]')
158+
pl.ylabel('Occurrences')
159+
rms_err_y = np.std(100*np.abs((yp_after-yp_before)-(yp[ii+1, :]-yp[0, :]))/np.std(yp[ii+1, :]-yp[0, :]))
160+
161+
162+
pl.suptitle('Turn %d rms_err_x = %e rms_err_y = %e'%(ii, rms_err_x, rms_err_y))
163+
164+
pl.savefig(filename.split('_prb.dat')[0]+'_%02d.png'%ii, dpi=150)
165+
166+
rms_err_x_list.append(rms_err_x)
167+
rms_err_y_list.append(rms_err_y)
168+
169+
pl.ion()
170+
pl.draw()
171+
time.sleep(1.)
172+
173+
pl.figure(1000)
174+
pl.plot(rms_err_x_list, '.-', markersize = 10, linewidth=2, label='x')
175+
pl.plot(rms_err_y_list, '.-', markersize = 10, linewidth=2, label='y')
176+
pl.grid('on')
177+
pl.ylabel('''Relative r.m.s. error [%]''')
178+
pl.xlabel('Turn')
179+
pl.legend()
180+
pl.savefig(filename.split('_prb.dat')[0]+'_errors.png', dpi=200)
181+
pl.show()
182+
183+
slices = bunch.get_slices(ecloud.slicer)
184+
if show_movie:
185+
n_frames = ecloud.rho_ele_last_track.shape[0]
186+
pl.close(1)
187+
pl.figure(1, figsize=(8,8))
188+
vmax = np.max(ecloud.rho_ele_last_track[:])
189+
vmin = np.min(ecloud.rho_ele_last_track[:])
190+
for ii in xrange(n_frames-1, 0, -1):
191+
pl.subplot2grid((10,1),(0,0), rowspan=3)
192+
pl.plot(slices.z_centers, np.float_(slices.n_macroparticles_per_slice)/np.max(slices.n_macroparticles_per_slice))
193+
pl.xlabel('z [m]');pl.ylabel('Long. profile')
194+
pl.axvline(x = slices.z_centers[ii], color='r')
195+
pl.subplot2grid((10,1),(4,0), rowspan=6)
196+
pl.pcolormesh(ecloud.spacech_ele.xg*1e3, ecloud.spacech_ele.yg*1e3, ecloud.rho_ele_last_track[ii,:,:].T, vmax=vmax, vmin=vmin)
197+
pl.xlabel('x [mm]');pl.ylabel('y [mm]')
198+
pl.axis('equal')
199+
pl.subplots_adjust(hspace=.5)
200+
pl.draw()
201+
time.sleep(.2)
202+
pl.show()
203+
204+

0 commit comments

Comments
 (0)