Skip to content

Commit 03fdf93

Browse files
author
shi
committed
The comparison with the analytical solution for the uniformly charged sphere expansion added to the exmaples of 3D Space Charge lattice nodes.
1 parent d446207 commit 03fdf93

File tree

3 files changed

+366
-77
lines changed

3 files changed

+366
-77
lines changed

‎examples/SpaceCharge/sc3D/sc_2.5d_drift_latt_uniform_sphere_bunch.py‎

Lines changed: 125 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,8 @@
77
#---------------------------------------------------------
88
# 2.5D Rick Baartman's Space Charge model should not work
99
# for the linac case here it is just as example of setting
10-
# another SC model. N.B. Model is wrong!
10+
# another SC model.
11+
# N.B. Model is wrong!
1112
#---------------------------------------------------------
1213

1314
import sys
@@ -33,9 +34,64 @@
3334
#------------------------------
3435
random.seed(10)
3536

36-
#-------------------------------------
37-
# Longitudinal distribution functions
38-
#-------------------------------------
37+
#------------------------------
38+
# Auxilary functions
39+
#------------------------------
40+
41+
def func_theory(x):
42+
"""
43+
The function used in the charged sphere expansion dynamics calculation.
44+
See the presentation.
45+
"""
46+
val = math.sqrt(abs(x**2-1))
47+
val = math.log(x+val)/2 + x*val/2
48+
return val
49+
50+
def getR_Theory(s_path,r_init,R,Ntot,bunch):
51+
"""
52+
The root finding in the charged sphere expansion dynamics calculation.
53+
See the presentation.
54+
s_path - path length / drift length []
55+
r_init - initial distance from the sphere center [m]
56+
R - bunch sphere radius [m]
57+
Ntot - total macro-size of the bunch
58+
beta - relativistic beta
59+
r_classic - classical radius of particle
60+
"""
61+
r_classic = bunch.classicalRadius()
62+
gamma = bunch.getSyncParticle().gamma()
63+
beta = bunch.getSyncParticle().beta()
64+
val = math.sqrt((r_classic*Ntot/R**3)/2)*s_path/(gamma*beta)
65+
#---- Now we solve equation val = func_theory(x) for x
66+
x_min = 1.000000001
67+
x_max = 10.
68+
x = (x_min + x_max)/2.
69+
#print ("debug s_path=",s_path," val=",val)
70+
n_iter = 15
71+
count = 0
72+
f_min = func_theory(x_min)
73+
f_max = func_theory(x_max)
74+
f_val = func_theory((x_min + x_max)/2.)
75+
while(count < n_iter):
76+
x = (x_min + x_max)/2.
77+
f_val = func_theory(x)
78+
if(val <= f_val):
79+
f_max = f_val
80+
x_max = x
81+
count += 1
82+
continue
83+
if(val >= f_val):
84+
f_min = f_val
85+
x_min = x
86+
count += 1
87+
continue
88+
r_res = 1000.*r_init*x**2
89+
#print ("debug s_path=",s_path," val=",val," f_val=",f_val," x=",x," r_res=",r_res)
90+
return r_res
91+
92+
#----------------------------------------
93+
# Uniform sphere distribution functions
94+
#----------------------------------------
3995

4096
def getUniformXYZ(radius,gamma):
4197
"""
@@ -45,9 +101,9 @@ def getUniformXYZ(radius,gamma):
45101
(x,y,z) = (radius,radius,radius)
46102
r = math.sqrt(x**2 + y**2 + z**2)
47103
while(r > radius):
48-
x = radius*0.5*(1.0 - 2*random.random())
49-
y = radius*0.5*(1.0 - 2*random.random())
50-
z = radius*0.5*(1.0 - 2*random.random())
104+
x = radius*(1.0 - 2*random.random())
105+
y = radius*(1.0 - 2*random.random())
106+
z = radius*(1.0 - 2*random.random())
51107
r = math.sqrt(x**2 + y**2 + z**2)
52108
return (x,y,z/gamma)
53109

@@ -60,7 +116,7 @@ def getUniformXYZ(radius,gamma):
60116
bunch_length = 10.0 # m
61117
nParts = 1000000
62118
total_macroSize = 1.0e+11
63-
energy = 1.0 # GeV
119+
energy = 0.4 # GeV
64120

65121
syncPart = b_init.getSyncParticle()
66122
syncPart.kinEnergy(energy)
@@ -122,12 +178,6 @@ def getLattice(lattice_length,n_parts,calc3d,pipe_radius):
122178
#print ("debug n sc nodes=",len(scNodes_arr))
123179
return teapot_lattice
124180

125-
# -------------------------------------------------------------
126-
# set of uniformly charged ellipses Space Charge
127-
# -------------------------------------------------------------
128-
nEllipses = 1
129-
calc3d = SpaceChargeCalcUnifEllipse(nEllipses)
130-
131181
#-------------------------------------------------------------
132182
# set 2.5D Rick Baartman's Space Charge model.
133183
# This model should not work for linac case
@@ -136,10 +186,10 @@ def getLattice(lattice_length,n_parts,calc3d,pipe_radius):
136186
sizeX = 64
137187
sizeY = 64
138188
sizeZ = 64
139-
pipe_radius = 0.040
189+
pipe_radius = 0.100
140190
calc3d = SpaceChargeCalc2p5Drb(sizeX,sizeY,sizeZ)
141191

142-
lattice_length = 10.0 # the length of the drift
192+
lattice_length = 0.5 # the length of the drift
143193
n_parts = 30 # number of parts on what the drift will be chopped, or the number of SC nodes
144194
lattice = getLattice(lattice_length,n_parts,calc3d,pipe_radius)
145195

@@ -152,24 +202,69 @@ def getLattice(lattice_length,n_parts,calc3d,pipe_radius):
152202
bunch = Bunch()
153203
b_init.copyBunchTo(bunch)
154204

155-
#---------------------------------------
156-
#track the bunch through the lattice
157-
#---------------------------------------
158-
lattice.trackBunch(bunch)
205+
#-------------------------------------------------------------------
206+
#---- Here we add three particles - along x,y,z axises at
207+
#---- r_bunch/2 disstance from the center. We will use them to
208+
#---- compare analytical solution with the PyORBIT tracking,
209+
#---- N.B. : adding 3 particles did not change total charge much.
210+
#-------------------------------------------------------------------
211+
ip_arr = [nParticlesGlobal,nParticlesGlobal+1,nParticlesGlobal+2]
212+
r_in = 0.5*bunch_radius
213+
bunch.addParticle(-r_in,0.,0.,0.,0.,0.)
214+
bunch.addParticle(0.,0.,r_in,0.,0.,0.)
215+
bunch.addParticle(0.,0.,0.,0.,r_in/gamma,0.)
216+
r_in *= 1000 # now in [mm]
159217

160-
#----------------------------------------------
161-
#---- Analysis of the final bunch
162-
#----------------------------------------------
163-
twiss_analysis.analyzeBunch(bunch)
218+
#-------------------------------------------------------
219+
#---- Let's track bunch through the drift several times
220+
#-------------------------------------------------------
221+
s_path = 0.
164222

165-
x_rms = math.sqrt(twiss_analysis.getCorrelation(0,0)) * 1000.0
166-
y_rms = math.sqrt(twiss_analysis.getCorrelation(2,2)) * 1000.0
167-
z_rms = math.sqrt(twiss_analysis.getCorrelation(4,4)) * 1000.0
223+
print ("s[m] r_rms[mm] r1[mm] r2[mm] r3[mm] r_theory[mm] ")
168224

225+
for ind in range(20):
226+
#---------------------------------------
227+
#track the bunch through the lattice
228+
#---------------------------------------
229+
lattice.trackBunch(bunch)
230+
231+
s_path += lattice_length
232+
233+
#----------------------------------------------
234+
#---- Analysis of the final bunch
235+
#----------------------------------------------
236+
twiss_analysis.analyzeBunch(bunch)
237+
238+
x_rms = math.sqrt(twiss_analysis.getCorrelation(0,0)) * 1000.0
239+
y_rms = math.sqrt(twiss_analysis.getCorrelation(2,2)) * 1000.0
240+
z_rms = math.sqrt(twiss_analysis.getCorrelation(4,4)) * 1000.0
241+
r_rms = math.sqrt(x_rms**2 + y_rms**2 + (z_rms*gamma)**2)
242+
243+
st = ""
244+
r_out = 0.
245+
for ip in ip_arr:
246+
z = bunch.z(ip)*gamma
247+
x = bunch.x(ip)
248+
y = bunch.y(ip)
249+
xp = bunch.xp(ip)*1000
250+
r_out = 1000.*math.sqrt(x**2 + y**2 + z**2)
251+
st += " %6.3f "%(r_out)
252+
253+
254+
r_out_th = getR_Theory(s_path,0.5*bunch_radius,bunch_radius,total_macroSize,bunch)
255+
256+
#r_classic = bunch.classicalRadius()
257+
#gamma = bunch.getSyncParticle().gamma()
258+
#beta = bunch.getSyncParticle().beta()
259+
260+
#xp_th = r_classic*total_macroSize*s_path*r_in/(bunch_radius**3)/(gamma*beta)**2
261+
262+
print (" %4.3f %6.3f "%(s_path,r_rms),st," %6.3f "%r_out_th)
263+
169264
print ("==========After Traking=======================")
170-
print ("Initial x_rms[mm] = %6.5f x_rms/y_rms = %6.5f "%(x_rms,x_rms/y_rms))
171-
print ("Initial y_rms[mm] = %6.5f x_rms/(z_rms*gamma) = %6.5f "%(y_rms,x_rms/(z_rms*gamma)))
172-
print ("Initial z_rms[mm] = %6.5f x_rms/(z_rms*gamma) = %6.5f "%(z_rms,x_rms/(z_rms*gamma)))
265+
print ("x_rms[mm] = %6.5f x_rms/y_rms = %6.5f "%(x_rms,x_rms/y_rms))
266+
print ("y_rms[mm] = %6.5f y_rms/(z_rms*gamma) = %6.5f "%(y_rms,y_rms/(z_rms*gamma)))
267+
print ("z_rms[mm] = %6.5f x_rms/(z_rms*gamma) = %6.5f "%(z_rms,x_rms/(z_rms*gamma)))
173268
print ("=============================================")
174269

175270
print("Finish.")

‎examples/SpaceCharge/sc3D/sc_3D_fft_drift_latt_uniform_sphere_bunch.py‎

Lines changed: 122 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -29,9 +29,64 @@
2929
#------------------------------
3030
random.seed(10)
3131

32-
#-------------------------------------
33-
# Longitudinal distribution functions
34-
#-------------------------------------
32+
#------------------------------
33+
# Auxilary functions
34+
#------------------------------
35+
36+
def func_theory(x):
37+
"""
38+
The function used in the charged sphere expansion dynamics calculation.
39+
See the presentation.
40+
"""
41+
val = math.sqrt(abs(x**2-1))
42+
val = math.log(x+val)/2 + x*val/2
43+
return val
44+
45+
def getR_Theory(s_path,r_init,R,Ntot,bunch):
46+
"""
47+
The root finding in the charged sphere expansion dynamics calculation.
48+
See the presentation.
49+
s_path - path length / drift length []
50+
r_init - initial distance from the sphere center [m]
51+
R - bunch sphere radius [m]
52+
Ntot - total macro-size of the bunch
53+
beta - relativistic beta
54+
r_classic - classical radius of particle
55+
"""
56+
r_classic = bunch.classicalRadius()
57+
gamma = bunch.getSyncParticle().gamma()
58+
beta = bunch.getSyncParticle().beta()
59+
val = math.sqrt((r_classic*Ntot/R**3)/2)*s_path/(gamma*beta)
60+
#---- Now we solve equation val = func_theory(x) for x
61+
x_min = 1.000000001
62+
x_max = 10.
63+
x = (x_min + x_max)/2.
64+
#print ("debug s_path=",s_path," val=",val)
65+
n_iter = 15
66+
count = 0
67+
f_min = func_theory(x_min)
68+
f_max = func_theory(x_max)
69+
f_val = func_theory((x_min + x_max)/2.)
70+
while(count < n_iter):
71+
x = (x_min + x_max)/2.
72+
f_val = func_theory(x)
73+
if(val <= f_val):
74+
f_max = f_val
75+
x_max = x
76+
count += 1
77+
continue
78+
if(val >= f_val):
79+
f_min = f_val
80+
x_min = x
81+
count += 1
82+
continue
83+
r_res = 1000.*r_init*x**2
84+
#print ("debug s_path=",s_path," val=",val," f_val=",f_val," x=",x," r_res=",r_res)
85+
return r_res
86+
87+
#----------------------------------------
88+
# Uniform sphere distribution functions
89+
#----------------------------------------
3590

3691
def getUniformXYZ(radius,gamma):
3792
"""
@@ -41,9 +96,9 @@ def getUniformXYZ(radius,gamma):
4196
(x,y,z) = (radius,radius,radius)
4297
r = math.sqrt(x**2 + y**2 + z**2)
4398
while(r > radius):
44-
x = radius*0.5*(1.0 - 2*random.random())
45-
y = radius*0.5*(1.0 - 2*random.random())
46-
z = radius*0.5*(1.0 - 2*random.random())
99+
x = radius*(1.0 - 2*random.random())
100+
y = radius*(1.0 - 2*random.random())
101+
z = radius*(1.0 - 2*random.random())
47102
r = math.sqrt(x**2 + y**2 + z**2)
48103
return (x,y,z/gamma)
49104

@@ -56,7 +111,7 @@ def getUniformXYZ(radius,gamma):
56111
bunch_length = 10.0 # m
57112
nParts = 1000000
58113
total_macroSize = 1.0e+11
59-
energy = 1.0 # GeV
114+
energy = 0.4 # GeV
60115

61116
syncPart = b_init.getSyncParticle()
62117
syncPart.kinEnergy(energy)
@@ -127,7 +182,7 @@ def getLattice(lattice_length,n_parts,calc3d):
127182
calc3d = SpaceChargeCalc3D(sizeX,sizeY,sizeZ)
128183
print ("SC Solver - SpaceChargeCalc3D, grid sizes (X,Y,Z) = (%3d,%3d,%3d)"%(sizeX,sizeY,sizeZ))
129184

130-
lattice_length = 10.0 # the length of the drift
185+
lattice_length = 0.5 # the length of the drift
131186
n_parts = 30 # number of parts on what the drift will be chopped, or the number of SC nodes
132187
lattice = getLattice(lattice_length,n_parts,calc3d)
133188

@@ -140,24 +195,69 @@ def getLattice(lattice_length,n_parts,calc3d):
140195
bunch = Bunch()
141196
b_init.copyBunchTo(bunch)
142197

143-
#---------------------------------------
144-
#track the bunch through the lattice
145-
#---------------------------------------
146-
lattice.trackBunch(bunch)
198+
#-------------------------------------------------------------------
199+
#---- Here we add three particles - along x,y,z axises at
200+
#---- r_bunch/2 disstance from the center. We will use them to
201+
#---- compare analytical solution with the PyORBIT tracking,
202+
#---- N.B. : adding 3 particles did not change total charge much.
203+
#-------------------------------------------------------------------
204+
ip_arr = [nParticlesGlobal,nParticlesGlobal+1,nParticlesGlobal+2]
205+
r_in = 0.5*bunch_radius
206+
bunch.addParticle(-r_in,0.,0.,0.,0.,0.)
207+
bunch.addParticle(0.,0.,r_in,0.,0.,0.)
208+
bunch.addParticle(0.,0.,0.,0.,r_in/gamma,0.)
209+
r_in *= 1000 # now in [mm]
147210

148-
#----------------------------------------------
149-
#---- Analysis of the final bunch
150-
#----------------------------------------------
151-
twiss_analysis.analyzeBunch(bunch)
211+
#-------------------------------------------------------
212+
#---- Let's track bunch through the drift several times
213+
#-------------------------------------------------------
214+
s_path = 0.
152215

153-
x_rms = math.sqrt(twiss_analysis.getCorrelation(0,0)) * 1000.0
154-
y_rms = math.sqrt(twiss_analysis.getCorrelation(2,2)) * 1000.0
155-
z_rms = math.sqrt(twiss_analysis.getCorrelation(4,4)) * 1000.0
216+
print ("s[m] r_rms[mm] r1[mm] r2[mm] r3[mm] r_theory[mm] ")
156217

218+
for ind in range(20):
219+
#---------------------------------------
220+
#track the bunch through the lattice
221+
#---------------------------------------
222+
lattice.trackBunch(bunch)
223+
224+
s_path += lattice_length
225+
226+
#----------------------------------------------
227+
#---- Analysis of the final bunch
228+
#----------------------------------------------
229+
twiss_analysis.analyzeBunch(bunch)
230+
231+
x_rms = math.sqrt(twiss_analysis.getCorrelation(0,0)) * 1000.0
232+
y_rms = math.sqrt(twiss_analysis.getCorrelation(2,2)) * 1000.0
233+
z_rms = math.sqrt(twiss_analysis.getCorrelation(4,4)) * 1000.0
234+
r_rms = math.sqrt(x_rms**2 + y_rms**2 + (z_rms*gamma)**2)
235+
236+
st = ""
237+
r_out = 0.
238+
for ip in ip_arr:
239+
z = bunch.z(ip)*gamma
240+
x = bunch.x(ip)
241+
y = bunch.y(ip)
242+
xp = bunch.xp(ip)*1000
243+
r_out = 1000.*math.sqrt(x**2 + y**2 + z**2)
244+
st += " %6.3f "%(r_out)
245+
246+
247+
r_out_th = getR_Theory(s_path,0.5*bunch_radius,bunch_radius,total_macroSize,bunch)
248+
249+
#r_classic = bunch.classicalRadius()
250+
#gamma = bunch.getSyncParticle().gamma()
251+
#beta = bunch.getSyncParticle().beta()
252+
253+
#xp_th = r_classic*total_macroSize*s_path*r_in/(bunch_radius**3)/(gamma*beta)**2
254+
255+
print (" %4.3f %6.3f "%(s_path,r_rms),st," %6.3f "%r_out_th)
256+
157257
print ("==========After Traking=======================")
158-
print ("Initial x_rms[mm] = %6.5f x_rms/y_rms = %6.5f "%(x_rms,x_rms/y_rms))
159-
print ("Initial y_rms[mm] = %6.5f x_rms/(z_rms*gamma) = %6.5f "%(y_rms,x_rms/(z_rms*gamma)))
160-
print ("Initial z_rms[mm] = %6.5f x_rms/(z_rms*gamma) = %6.5f "%(z_rms,x_rms/(z_rms*gamma)))
258+
print ("x_rms[mm] = %6.5f x_rms/y_rms = %6.5f "%(x_rms,x_rms/y_rms))
259+
print ("y_rms[mm] = %6.5f y_rms/(z_rms*gamma) = %6.5f "%(y_rms,y_rms/(z_rms*gamma)))
260+
print ("z_rms[mm] = %6.5f x_rms/(z_rms*gamma) = %6.5f "%(z_rms,x_rms/(z_rms*gamma)))
161261
print ("=============================================")
162262

163263
print("Finish.")

0 commit comments

Comments
 (0)