Skip to content

Commit

Permalink
Beam -> bunch
Browse files Browse the repository at this point in the history
  • Loading branch information
mattsignorelli committed Nov 14, 2024
1 parent 51f1946 commit 5ee7500
Show file tree
Hide file tree
Showing 10 changed files with 74 additions and 74 deletions.
20 changes: 10 additions & 10 deletions docs/src/devel.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,12 +5,12 @@ The design of `BeamTracking.jl` is built with modularity, high performance, diff
The entire package is centered around one single `track!` function, which has the following format:

```julia
track!(beam::Beam, element, work=work)
track!(bunch::Bunch, element, work=work)
```

Here, `beam` is a `Beam` struct, which is described in detail [below](@ref beam). `element` is some element which to track the beam through, and `work` is an optional tuple of the minimal number of temporaries needed for use inside the tracking function (for some elements, it is an empty tuple).
Here, `bunch` is a `Bunch` struct, which is described in detail [below](@ref bunch). `element` is some element which to track the bunch through, and `work` is an optional tuple of the minimal number of temporaries needed for use inside the tracking function (for some elements, it is an empty tuple).

After calling `track!`, the `beam` struct is mutated to contain the particle phase space coordinates (and possiblly spin transport quaternions) after propagation through the `element`. With this `track!` function, all one needs to do is define their own custom `element` type, and then when looping through a vector of elements, Julia's multiple dispatch will take care of each particular element.
After calling `track!`, the `bunch` struct is mutated to contain the particle phase space coordinates (and possiblly spin transport quaternions) after propagation through the `element`. With this `track!` function, all one needs to do is define their own custom `element` type, and then when looping through a vector of elements, Julia's multiple dispatch will take care of each particular element.

For example, using the `Linear.Drift` and `Linear.Quadrupole` elements:

Expand All @@ -21,17 +21,17 @@ qd = Linear.Quadrupole(Bn1=18.5, L=0.5)

fodo = (qf, d, qd, d)

beam = Beam(x=1e-3, px=1e-4, pz=1e-4, beta_gamma_ref=35000.) # Creates a Beam with one particle
bunch = Bunch(x=1e-3, px=1e-4, pz=1e-4, beta_gamma_ref=35000.) # Creates a Bunch with one particle
for ele in fodo
track!(beam, ele)
track!(bunch, ele)
end
```

## [The `Beam` Struct](@id beam)
## [The `Bunch` Struct](@id bunch)

The `Beam` struct contains the following fields:
The `Bunch` struct contains the following fields:

- **`species::Species`** -- A `Species` type storing the beam species (e.g. electron)
- **`species::Species`** -- A `Species` type storing the bunch species (e.g. electron)
- **`beta_gamma_ref::Float64`** -- A reference Lorentz $\beta\gamma$ for normalizing the transverse momenta to
- **`v::T`** -- All particles' 6D phase space coordinates stored in a structure-of-arrays (SoA) memory layout
- **`q::U`** -- If spin tracking, then all particles' spin transport quaternions (using the quaternion type defined in [`ReferenceFrameRotations.jl`](https://github.com/JuliaSpace/ReferenceFrameRotations.jl)) stored in a structure-of-arrays layout. Else, `nothing`
Expand Down Expand Up @@ -70,13 +70,13 @@ v.x # accesses x array
v.px # accesses px array
v[1] # This goes from SoA to a single Coord struct representing the first Coord!
```
In the above example, `v` is what is stored in the `Beam` struct. The choice of SoA was made after careful speed benchmarks and the desire to have a mutable interface (`track!` instead of `track.` for a beam).
In the above example, `v` is what is stored in the `Bunch` struct. The choice of SoA was made after careful speed benchmarks and the desire to have a mutable interface (`track!` instead of `track.` for a bunch).

When there spin tracking, `q` is a `StructArray{<:Quaternion}`

## Rules for `track!` Implementations

Vectorized operations should be used. The function should be type-stable and, when pre-allocating the necessary `work`, have zero allocations included when tracking a single _non-parametric_ GTPSA map or tracking a beam of particles as immutable numbers (`Float64` or `Dual` numbers for example). For parametric GTPSA maps (e.g. when a quadrupole strength is included as a parameter in the GTPSA), the allocation restriction is loosened in order to maintain readable code. Tests for this are included in the `test_matrix` function in [`runtests.jl`](https://github.com/bmad-sim/BeamTracking.jl/blob/main/test/runtests.jl)
Vectorized operations should be used. The function should be type-stable and, when pre-allocating the necessary `work`, have zero allocations included when tracking a single _non-parametric_ GTPSA map or tracking a bunch of particles as immutable numbers (`Float64` or `Dual` numbers for example). For parametric GTPSA maps (e.g. when a quadrupole strength is included as a parameter in the GTPSA), the allocation restriction is loosened in order to maintain readable code. Tests for this are included in the `test_matrix` function in [`runtests.jl`](https://github.com/bmad-sim/BeamTracking.jl/blob/main/test/runtests.jl)

## Compatibility with `GTPSA.jl`

Expand Down
4 changes: 2 additions & 2 deletions src/BeamTracking.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ using GTPSA,

include("aapc.jl")

export Beam,
export Bunch,
Coord,
Quaternion,
Particle,
Expand Down Expand Up @@ -41,7 +41,7 @@ export Beam,
include("types.jl")

# Empty tracking method to be imported by submodules
track!(beam::Beam, ::Nothing; work=nothing) = beam
track!(bunch::Bunch, ::Nothing; work=nothing) = bunch

# --------------------------------------------------

Expand Down
22 changes: 11 additions & 11 deletions src/Linear/Linear.jl
Original file line number Diff line number Diff line change
Expand Up @@ -51,11 +51,11 @@ Base.@kwdef struct Solenoid{T}
end


function track!(beam::Beam, ele::Linear.Drift; work=get_work(beam, Val{0}()))
function track!(bunch::Bunch, ele::Linear.Drift; work=get_work(bunch, Val{0}()))
L = ele.L
v = beam.v
v = bunch.v

gamma_ref = sr_gamma(beam.beta_gamma_ref)
gamma_ref = sr_gamma(bunch.beta_gamma_ref)

@FastGTPSA! begin
@. v.x = v.x + L*v.px
Expand All @@ -65,16 +65,16 @@ function track!(beam::Beam, ele::Linear.Drift; work=get_work(beam, Val{0}()))

# Spin unchanged

return beam
return bunch
end


function track!(beam::Beam, ele::Linear.Quadrupole; work=get_work(beam, Val{1}()))
v = beam.v
function track!(bunch::Bunch, ele::Linear.Quadrupole; work=get_work(bunch, Val{1}()))
v = bunch.v
L = ele.L

Kn1 = ele.Bn1 / brho(massof(beam.species), beam.beta_gamma_ref, chargeof(beam.species))
gamma_ref = sr_gamma(beam.beta_gamma_ref)
Kn1 = ele.Bn1 / brho(massof(bunch.species), bunch.beta_gamma_ref, chargeof(bunch.species))
gamma_ref = sr_gamma(bunch.beta_gamma_ref)

if Kn1 >= 0
k = sqrt(Kn1)
Expand Down Expand Up @@ -112,13 +112,13 @@ function track!(beam::Beam, ele::Linear.Quadrupole; work=get_work(beam, Val{1}()
@. v.z = v.z + L/gamma_ref^2*v.pz
end

return beam
return bunch
end

"""
Routine to linearly tracking through a Sector Magnet
"""
function track!(beamf::Beam, ele::Linear.SBend, beami::Beam)
function track!(beamf::Bunch, ele::Linear.SBend, beami::Bunch)
@assert !(beamf === beami) "Aliasing beamf === beami not allowed!"
zi = beami.vec
zf = beamf.vec
Expand All @@ -144,7 +144,7 @@ end
"""
Routine to linearly tracking through a Combined Magnet
"""
function track!(beamf::Beam, ele::Linear.Combined, beami::Beam)
function track!(beamf::Bunch, ele::Linear.Combined, beami::Bunch)
@assert !(beamf === beami) "Aliasing beamf === beami not allowed!"
zi = beami.vec
zf = beamf.vec
Expand Down
22 changes: 11 additions & 11 deletions src/MatrixKick/MatrixKick.jl
Original file line number Diff line number Diff line change
Expand Up @@ -15,12 +15,12 @@ Base.@kwdef struct Quadrupole{T}
end


function track!(beam::Beam, ele::MatrixKick.Drift; work=get_work(beam, Val{1}()))
function track!(bunch::Bunch, ele::MatrixKick.Drift; work=get_work(bunch, Val{1}()))
L = ele.L
v = beam.v
v = bunch.v

tilde_m = 1 / beam.beta_gamma_ref
beta_ref = sr_beta(beam.beta_gamma_ref)
tilde_m = 1 / bunch.beta_gamma_ref
beta_ref = sr_beta(bunch.beta_gamma_ref)

@FastGTPSA! begin
@. work[1] = 1 / sqrt((1.0 + v.pz)^2 - (v.px^2 + v.py^2))
Expand All @@ -31,22 +31,22 @@ function track!(beam::Beam, ele::MatrixKick.Drift; work=get_work(beam, Val{1}())

# Spin unchanged

return beam
end # function track!(::Beam, ::Drift)
return bunch
end # function track!(::Bunch, ::Drift)



# This integrator uses the so-called Matrix-Kick-Matrix method to implement
# an integrator accurate though second-order in the integration step-size.
function track!(beam::Beam, ele::MatrixKick.Quadrupole; work=get_work(beam, Val{6}()))
function track!(bunch::Bunch, ele::MatrixKick.Quadrupole; work=get_work(bunch, Val{6}()))
@assert !(beamf === beami) "Aliasing beamf === beami not allowed!"
L = ele.L

# κ^2 (kappa-squared) := (q g / P0) / (1 + δ)
# numerator of κ^2
k2_num = Bn1 / brho(massof(beam.species), beam.beta_gamma_ref, chargeof(beam.species))
k2_num = Bn1 / brho(massof(bunch.species), bunch.beta_gamma_ref, chargeof(bunch.species))

v = beam.v
v = bunch.v
v_work = StructArray{Coord{eltype(work[1])}}((work[1], work[2], work[3], work[4], work[5], work[6]))

trackQuadMx!(v_work, v, k2_num, L / 2)
Expand All @@ -56,12 +56,12 @@ function track!(beam::Beam, ele::MatrixKick.Quadrupole; work=get_work(beam, Val{

v .= v_work

return beam
return bunch
end # function track!(::Quadrupole)


"""
trackQuadMx!(beamf::Beam, beami::Beam, k2_num::Float64, s::Float64)
trackQuadMx!(beamf::Bunch, beami::Bunch, k2_num::Float64, s::Float64)
track "matrix part" of quadrupole
"""
Expand Down
2 changes: 1 addition & 1 deletion src/Misc/Misc.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ Base.@kwdef struct Taylor{T<:TPS,U<:Union{Nothing,Quaternion{<:T}}}
q::U
end

function track!(beam::Beam, ele::Misc.Taylor)
function track!(bunch::Bunch, ele::Misc.Taylor)

end

Expand Down
32 changes: 16 additions & 16 deletions src/types.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
#=
Beam, Coord, and Particle type definitions. StructArrays is
Bunch, Coord, and Particle type definitions. StructArrays is
used to handle an internal SoA layout of memory which also
allows us to mutate, so both the phase space Coord struct
(defined here) and Quaternion struct (defined in
Expand All @@ -21,30 +21,30 @@ Base.@kwdef struct Coord{T} <: FieldVector{6, T}
pz::T = 0.0
end

# Static quaternion type defined by ReferenceFrameRotations
# Static quaternion type defined by ReferenceFrameRotatio

struct Beam{T<:StructVector{<:Coord}, U<:Union{Nothing, StructVector{<:Quaternion}}}
struct Bunch{T<:StructVector{<:Coord}, U<:Union{Nothing, StructVector{<:Quaternion}}}
species::Species
beta_gamma_ref::Float64
v::T
q::U
end

# Initialize a Beam with either a single particle (scalars)
# Initialize a Bunch with either a single particle (scalars)
"""
Beam(; species::Species=Species("electron"), beta_gamma_ref=1.0,
Bunch(; species::Species=Species("electron"), beta_gamma_ref=1.0,
spin::Union{Bool,Nothing}=nothing, gtpsa_map::Union{Bool,Nothing}=nothing,
x::Union{Number,AbstractVector}=0.0, px::Union{Number,AbstractVector}=0.0,
y::Union{Number,AbstractVector}=0.0, py::Union{Number,AbstractVector}=0.0,
z::Union{Number,AbstractVector}=0.0, pz::Union{Number,AbstractVector}=0.0 )
Initializes a `Beam`. Any of the specified phase space coordinates may be scalar `Number`s or
Initializes a `Bunch`. Any of the specified phase space coordinates may be scalar `Number`s or
`Vector`s to store as a structure-of-arrays. Internally, all phase space coordinates are stored
as `Vector`s. If all phase space coordinates are scalar `Number`s, then a `Beam` is created with
as `Vector`s. If all phase space coordinates are scalar `Number`s, then a `Bunch` is created with
a single particle. If any of the coordinates are specified as `Vector`s, then all other scalar-
specified quantities are `fill`-ed as `Vector`s. For example, `Beam(x=1.0, y=[2,3])` creates a
beam with two particles having the phase space coordinates `[1.0, 0.0, 2.0, 0.0, 0.0, 0.0]`
specified quantities are `fill`-ed as `Vector`s. For example, `Bunch(x=1.0, y=[2,3])` creates a
bunch with two particles having the phase space coordinates `[1.0, 0.0, 2.0, 0.0, 0.0, 0.0]`
and `[1.0, 0.0, 3.0, 0.0, 0.0, 0.0]`.
## Other keyword arguments
Expand All @@ -53,12 +53,12 @@ and `[1.0, 0.0, 3.0, 0.0, 0.0, 0.0]`.
• `spin` -- If true, spin tracking is turned on and a quaternion for each particle is tracked
• `gtpsa_map` -- If true, GTPSA map tracking is used for each particle using the Descriptor defined in GTPSA.desc_current
"""
function Beam(; species::Species=Species("electron"), beta_gamma_ref=1.0,
function Bunch(; species::Species=Species("electron"), beta_gamma_ref=1.0,
spin::Union{Bool,Nothing}=nothing, gtpsa_map::Union{Bool,Nothing}=nothing,
x::Union{Number,AbstractVector}=0.0, px::Union{Number,AbstractVector}=0.0,
y::Union{Number,AbstractVector}=0.0, py::Union{Number,AbstractVector}=0.0,
z::Union{Number,AbstractVector}=0.0, pz::Union{Number,AbstractVector}=0.0 )

idx_vector = findfirst(t->t isa AbstractVector, (x, px, y, py, z, pz))
if isnothing(idx_vector)
N_particle = 1
Expand Down Expand Up @@ -131,7 +131,7 @@ function Beam(; species::Species=Species("electron"), beta_gamma_ref=1.0,
end


return Beam(species, Float64(beta_gamma_ref), v, q)
return Bunch(species, Float64(beta_gamma_ref), v, q)
end

struct Particle{T,U<:Union{Nothing,Quaternion{T}}}
Expand All @@ -141,8 +141,8 @@ struct Particle{T,U<:Union{Nothing,Quaternion{T}}}
q::U
end

function Particle(beam::Beam, idx::Integer=1)
v = beam.v[idx] # StructArrays handles this!
q = isnothing(beam.q) ? nothing : beam.q[idx]
return Particle(beam.species, beam.beta_gamma_ref, v, q)
function Particle(bunch::Bunch, idx::Integer=1)
v = bunch.v[idx] # StructArrays handles this!
q = isnothing(bunch.q) ? nothing : bunch.q[idx]
return Particle(bunch.species, bunch.beta_gamma_ref, v, q)
end
14 changes: 7 additions & 7 deletions src/work.jl
Original file line number Diff line number Diff line change
@@ -1,17 +1,17 @@
"""
get_work(beam::Beam, ::Val{N}) where N -> work
get_work(bunch::Bunch, ::Val{N}) where N -> work
Returns a tuple of `N` arrays of type `eltype(Beam.v.x)` and
length `length(Beam.v.x)` which may be used as temporaries.
Returns a tuple of `N` arrays of type `eltype(Bunch.v.x)` and
length `length(Bunch.v.x)` which may be used as temporaries.
### Arguments
- `beam` -- Beam to extract types and number of particles from
- `bunch` -- Bunch to extract types and number of particles from
- `::Val{N}` -- Number of `N` temporary arrays desired
"""
function get_work(beam::Beam, ::Val{N}) where {N}
sample = first(beam.v.x)
function get_work(bunch::Bunch, ::Val{N}) where {N}
sample = first(bunch.v.x)
T = typeof(sample)
N_particle = length(beam.v.x)
N_particle = length(bunch.v.x)

# Instead of using zeros, we do this to ensure
# same GTPSA descriptor if T isa TPS.
Expand Down
4 changes: 2 additions & 2 deletions test/element_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ zf1.y .= [0., 0., 0., 1.e-3,
zf1.py .= [0., 0., 0., 0., 3.e-4, -3.e-4, 3.e-4, 3.e-4 ]
zf1.z .= [1.4999999141696263-4, 0., -1.499999914426978-4, 1.4999999141696263-4, 1.4995115162148818-4, 1.4995115162148818-4, 2.9990230324297636-4, 2.9990230324297636-4]
zf1.pz .= [1.e-3, 0., -1.e-3, 1.e-3, 1.e-3, 1.e-3, 1.e-3, 1.e-3 ]
beami1 = Beam(e_minus, bg1, zi1)
beamf1 = Beam(e_minus, bg1, zf1)
beami1 = Bunch(e_minus, bg1, zi1)
beamf1 = Bunch(e_minus, bg1, zf1)


# test individual elements
Expand Down
16 changes: 8 additions & 8 deletions test/linear.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,15 +17,15 @@ using GTPSA
#Sbend
#sb_ele1 = Linear.SBend(L = 1.5, B0 = 0.000446) # positive field

#single particle beam
#single particle bunch
d = Descriptor(6,1)
bi = Beam(d)
bf_d = Beam(d)
bf_q1 = Beam(d)
bf_q2 = Beam(d)
bf_q3 = Beam(d)
bf_q4 = Beam(d)
bf_s1 = Beam(d)
bi = Bunch(d)
bf_d = Bunch(d)
bf_q1 = Bunch(d)
bf_q2 = Bunch(d)
bf_q3 = Bunch(d)
bf_q4 = Bunch(d)
bf_s1 = Bunch(d)

#Drift
track!(bf_d, drift_ele, bi)
Expand Down
12 changes: 6 additions & 6 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,21 +13,21 @@ function test_matrix(ele, n_work, M_expected; type_stable=true, no_allocs=true,
beta_gamma_ref=1.0, species=Species("electron") )

GTPSA.desc_current = Descriptor(6, 1) # 6 variables, order 1
beam = Beam(beta_gamma_ref=beta_gamma_ref, species=species, gtpsa_map=true)
work = BeamTracking.get_work(beam, Val(n_work))
bunch = Bunch(beta_gamma_ref=beta_gamma_ref, species=species, gtpsa_map=true)
work = BeamTracking.get_work(bunch, Val(n_work))

track!(beam, ele, work=work)
p = Particle(beam, 1)
track!(bunch, ele, work=work)
p = Particle(bunch, 1)

# 1) Correctness
@test norm(GTPSA.jacobian(p.v) - M_expected) < tol
# 2) Type stability
if type_stable
@test_opt track!(beam, ele, work=work)
@test_opt track!(bunch, ele, work=work)
end
# 3) No Allocations
if no_allocs
@test @ballocated(track!($beam, $ele, work=$work)) == 0
@test @ballocated(track!($bunch, $ele, work=$work)) == 0
end
end

Expand Down

0 comments on commit 5ee7500

Please sign in to comment.