J-PARC Trace3D Upgrades
Christopher K. Allen
Los Alamos National Laboratory

Trace3D “Project”
 Split the big FORTRAN file
 Each SUBROUTINE has file
 Project built using a Makefile (K. Furukawa) Version controlled (K. Furukawa) Repository jkksv01.j-parc.jp:/jk/master/t3d_ckaCan reassemble back into the big FORTRAN
file if desiredGeneral Additions
Debugging and Data I/O
Added output subroutines for debugging and recording Twiss
parameters along beamline
 InputDump – stores input “deck” to disk (file “InputData.txt”) tocheck that data is being properly input

TwissDump – stores Twiss parameters along beamline to disk file
“TwissDump.txt”. Twiss parameters are stored during simulation
since trajectory data is not kept in Trace3D memory.

Convenience Matrix Functions
Multiply, scalar multiply, commutator
Matrix exponent
Matrix logarithm (E. Forest)

General Additions (cont.) Steering Magnets

Added steering magnets to the element library.
 Note that only the centroid tracking is affectedby this element, it has no effect on the secondorder moment dynamics.
 The type identifier is 19 (NT=19).
 The parameters is Δx’, Δy’General Additions (cont.)
 Arbitrarily Oriented Beam Ellipsoids

To compute the electric self forces of the beam
Trace3D performs a coordinate transform to“eigen-coordinates” of ellipsoid In the current version this transform fails when

the beam is skewed arbitrarily off the design axis
 That is, skewed in a direction other that one of thetwo transverse directions.

 This condition was fixed.

Choice of Elliptic Integral RD or Form Factor ξ For the space charge calculations you mayswitch between the use of the form factor anddirect (numerical) evaluation of the ellipticintegral.
 This feature is available using the input variable

iFlgUseRd in the $DATA section of the Trace3Dinput file iFlgUseRd = 0: use form factor ξ (default)
 iFlgUseRd = 1: use elliptic integral RD

Form Factor vs. Elliptic Integral
Original Trace3D uses a “form factor” ξ in the self fieldcalculations.
λThe function ξ(s) is part of an analytic approximation to an
elliptic integral RD, which is defined
R D ( x, y , z ) ≡ξ(s)
1


0.8

3 ∞
dt
2 0 (t + x )1 / 2 (t + y )1 / 2 (t + z )3 / 2∫The form factor is defineds
1−cos −1 s

2s ∞
dt
1 
1−s
ξ (s) ≡=
s2 0 (t + 1) t + s 2 3 / 2 1 − s 2 1 −cosh −1 ss2 − 1

0.6

∫

0.4
0.2
12

345( )

sfor s < 1for s > 1

Form Factor vs. Elliptic Integral (cont.)
The approximations for RD in terms of the form factor are then given as
RD ( Z 2 , Y 2 , X 2 ) = Z 
31  + O(ε ),1−ξXZ X + Y 
 XY RD ( Z 2 , X 2 , Y 2 ) ≈ Z 
31  + O(ε ),1−ξ
YZ X + Y 
 XY RD ( X 2 , Y 2 , Z 2 ) ≈3  Z  + O(ε 2 ),ξ 
XYZ  XY 

Where X, Y, and Z are the semi-axes of the beam ellipsoid andε = (X-Y)/2is related to the eccentricity of the transverse ellipse

0.65
1.2s (m)

2.08

3.51

3.42

3.42

3.3

3.243.09

1.2
1.37

432

108s (m)765

BETA_x
BETA_y
BETA_z

3.51

3.42

3.42

3.3

3.24

3.09

3.03

2.94

2.91

2.88

2.26

2.15

2.15

2.08

2.01

1.95

1.68

1.61

1.51

1.44

1.41

SNS MEBT - SNS Trace (Using Form Factor)

3.03

2.94

2.91

2.88

2.26

2.15

2.15s (m)

2.01

1.95

1.68

1.61

1.51

1.44

1.41

1.37

0.73

SNS MEBT - KEK Trace3D (Using Rd)

0.73

0.65

00.54

00.43

10.65

10.54

20.65

20.28

30.14

30.35

SNS Version of Trace3D (form fac)
40.35

50.43

5060.14

60.28

70.14

7

beta (mm)80

Beta (mm)80.14

Beta (mm)

43.51

3.42

3.42

3.3

3.24

3.09

3.03

2.94

2.91

2.88

2.26

2.15

2.15

2.08

2.01

1.95

1.68

1.61

1.51

1.44

1.41

1.37

1.2

0.73

0.65

0.65

0.54

0.43

0.35

0.28

0.14

0.14

0

Comparison of Trace3D elliptic integral
and form factor simulations
SNS MEBT - KEK Trace3D (Using Form Factor)

BETA_x
BETA_y
BETA_z
BETA_x
BETA_y
BETA_zAdaptive Integration SteppingApproachλλForm a transfer matrix Φ(s;s0) that includes space effects tosecond order (2nd order accurate)
Choose error tolerance ε in the solution (~ 10-5 to 10-7)
Use Φ(s;s0) to propagate σ in steps h whose length isdetermined adaptively to maintain εAdaptive Stepping and Trace3D

Due to Trace3D “architecture”
implementing adaptive
stepping may require a major
rewrite

Brittle – dangerous

Implementation possible if it
can be done inSUBROUTINE TRANSλλ

Compute log(Φ) = ΔsA
Compute exp(hA)
May be too CPU intensive

XAL Architectureϒ

XAL Architecture ismodernλλλϒ

Not coupled
Easier to upgrade
Easier to maintain

XAL Appl. DevelλλλHard part is setting up DB
API is object oriented anddocumented (easy)
New features are(relatively) easily installed

XAL Architecture – Sequence DiagramSummary Adding adaptive stepping to Trace3D may not beworth the effort
Ability to implement it in TRANS
 Compute matrix logarithms and exponentials
 Add space charge to SAD?
Unknown effort – inexperienced with SAD

 Add J-PARC features to XAL?

Is XAL a player?
 Very familiar with XAL4. Space Charge Algorithm (cont.)
Assume that A and B are constantϒThen full transfer matrix Φ(s;s0) = e(s-s0)(A+B)

For practical reasons, we are usually given ΦA and ΦB
λλBeam optics provides ΦA(s)In ellipsoid coordinates ΦB(s) has simply form because B2 = 0
ΦB(s) = esB = I + sB

Define Φave(s) = _[ ΦA(s)ΦB(s) + ΦB(s)ΦA(s) ]λ By Taylor expanding Φ(s) = es(A+B) we findΦ(s) = Φave(s) + O(s3)That is, Φave(s) is a second-order accurate approximation of Φ(s)

4. Space Charge Algorithm (cont.)ϒ

We now have a stepping procedure which is second-orderaccurate in step length h.
Now consider the effects of “step doubling”
λλLet τ1(s+2h) denote the result of taking one step of length 2hLet τ2(s+2h) denote the result of taking twos steps of length hΔ(h) ≡ τ1(s+2h) - τ2(s+2h) = 6ch3
where the constant c = dτ(s’)/ds for some s’∈[s,s+2h]

ϒConsider ratio of |Δ| for steps of differing lengths h0 and h1
h1= h0 [|Δ1|/ |Δ0|]1/3

4. Space Charge Algorithm (cont.)
We use the formula h1= h0 [|Δ1|/ |Δ0|]1/3 as the basis for adaptive step sizing
Given |Δ1| = ε, a prescribed solution residual error we can tolerate
For each iteration kλλλ

Let |Δ0| = | τ1(sk+2h) - τ2(sk+2h) |, the residual error
Let h0 = hk be the step size at iteration kLet h1 = hk+1 be the step size at iteration k+1

hk+1 = hk [ε/| τ1(sk+2h) - τ2(sk+2h) |]1/3
where if hk+1 < hk, we must re-compute the kth step using the new steps size hk+1
to maintain the same solution accuracy