# CUDA.jl for particle tracking simulation

**URL:** <https://discourse.julialang.org/t/cuda-jl-for-particle-tracking-simulation/110906>\
**Category:** GPU\
**Tags:** speed-optimization, simulations\
**Created:** [February 28, 2024, 4:28pm UTC](https://discourse.julialang.org/t/cuda-jl-for-particle-tracking-simulation/110906 "2024-02-28T16:28:59Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![VitorSouzaLNLS](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vitorsouzalnls/32/207341_2.png) [@VitorSouzaLNLS](https://discourse.julialang.org/u/VitorSouzaLNLS)\
**Post date:** [February 28, 2024, 4:28pm UTC](https://discourse.julialang.org/t/cuda-jl-for-particle-tracking-simulation/110906/1 "2024-02-28T16:28:59Z")

</div>

Im trying to implement some consolidated CPU functions for particle tracking simulations (models of synchrotron accelerators) in GPU using CUDA.jl. The goal is to be able to paralelize tracking routines.

I had tried to start with the simplest function: a drift.

In CPU, the drift function is define as:

```julia
function drift(pos::Pos{T}, length::Float64)
    pnorm::T = 1 / (1 + pos.de)
    norml::T = pnorm * length
    pos.rx += norml * pos.px
    pos.ry += norml * pos.py
    pos.dl += 0.5 * norml * pnorm * (pos.px^2 + pos.py^2)
end 

```

Where the object `Pos` (a single electron state) is basicaly:

```julia
mutable struct Pos{T}
    rx::T
    px::T
    ry::T
    py::T
    de::T
    dl::T
end 

```

My GPU first implementation of the drift function is:

```julia
function drift(pos, length)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    if i <= size(pos)[2]
        @inbounds pnorm = 1 / (1 + pos[5, i])
        norml = pnorm * length
        @inbounds pos[1, i] += norml * pos[2,i]
        @inbounds pos[3, i] += norml * pos[4,i]
        @inbounds pos[6, i] += 0.5 * norml * pnorm * (pos[2,i]*pos[2,i] + pos[4,i]*pos[4,i])
    end
    return
end

```

Where the `pos` argument is now a `CuArray` (`CuDeviceVector`) with the dimension = (6, number\_of\_particles).

Then, I wrote the following function (`multiple_drifts`) for tracking the particles allong a line of multiple drifts.

```julia
function multiple_drifts(pos, lengths)
    for l in lengths
        drift(pos, l)
    end
    return
end

```

To do the comparison between the CPU and GPU functions i did:

```julia
dim::Int = 1 # the number of electrons
gpulengths = CUDA.rand(Float64, 6000*20) # the lengths of the drifts to pass
cpulengths = Array(gpulengths) # the CPU array of the lengths

nthreads = 256;
nblocks = Int(ceil(dim/256));

# GPU drift pass
# creation of the particle (all rx, px, ry, py, de, dl = 1e-6) {Float64}
particles = CUDA.fill(1e-6, (6, dim)); 
# the computation
t_gpu = CUDA.@elapsed begin CUDA.@sync @cuda threads=nthreads blocks=nblocks multiple_drifts(particles, lengs) end; 

#CPU drift pass
# creation of the particle (all rx, px, ry, py, de, dl = 1e-6) {Float64}
particle = Pos(1e-6); # using my constructor
# creation of the accelerator model
model = Accelerator()
# adding the lattice
model.lattice = [Elements.drift(length = l) for l in cpulengths] 
# the computation
t_cpu = @elapsed begin p_final = line_pass(model, p) end;

```

The final state of the particles is equally right… but the time is much slower in the GPU (35 ~ 40 times slower). The benchmarks are:

```julia
# CPU benchmark
# code line: p = Pos(1e-6); @benchmark begin line_pass(acc, $p, "end") end
BenchmarkTools.Trial: 1428 samples with 1 evaluation.
 Range (min … max): 2.686 ms … 21.966 ms ┊ GC (min … max): 0.00% … 65.95%
 Time (median): 3.275 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 3.487 ms ± 2.084 ms ┊ GC (mean ± σ): 7.82% ± 10.94%

  █ █                                                         
  █▅█▄▄▃▂▂▂▂▁▁▂▁▂▁▂▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▁▂▂▁▁▂▁▂▂▁▂▂ ▂
  2.69 ms Histogram: frequency by time 17.9 ms <

 Memory estimate: 1.96 MiB, allocs estimate: 120009.

# GPU benchmark
# code line: particles = CUDA.fill(1e-6, (6, dim)); @benchmark CUDA.@sync @cuda threads=nthreads blocks=nblocks multiple_drifts($particles, $lengs)
BenchmarkTools.Trial: 41 samples with 1 evaluation.
 Range (min … max): 122.164 ms … 124.992 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 123.547 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 123.516 ms ± 818.585 μs ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▃ ▃ ▃▃ ▃ ▃ ▃▃ █                
  █▇▁▁▁█▁▁██▁█▁▁▁▇▇▁▁▁▇▇▁▇▇▁▇▁█▇▁██▇▁▇▁▁▇▁▁▁▇▁▁▇█▇▇▁▇▇▇▇▁▁▁▇▁▁▇ ▁
  122 ms Histogram: frequency by time 125 ms <

 Memory estimate: 4.03 KiB, allocs estimate: 72.

```

I am really newbie in CUDA.jl and GPU implementations. I researched other things like shared memory, type stabilization, atomic operations etc… but I couldnt solve my problem.

Any help or suggestion?

(If you want more details about the CPU implementations like the line\_pass function or the Accelerator and Element structures, please tell me! Ill be happy to clarify.

---

<div class="post-metadata">

**Author:** ![trahflow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trahflow/32/30585_2.png) [@trahflow](https://discourse.julialang.org/u/trahflow)\
**Post date:** [February 28, 2024, 5:56pm UTC](https://discourse.julialang.org/t/cuda-jl-for-particle-tracking-simulation/110906/2 "2024-02-28T17:56:47Z")

</div>

Hi @VitorSouzaLNLS and welcome to this forum! 🙂

> [@VitorSouzaLNLS](#):
>
> but the time is much slower in the GPU (35 ~ 40 times slower).

this is somewhat expected. You’re doing a serial computation (drift for a number of steps) for just a single particle. That’s not something a GPU is good at. However if you would benchmark a drift of N \>\> 1 independent particles, the benchmark will look more favorably for the GPU, since it can do the computations for all (or rather many, whether it’s really _all_ depends on N) of the particles in parallel.

---

<div class="post-metadata">

**Author:** ![VitorSouzaLNLS](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vitorsouzalnls/32/207341_2.png) [@VitorSouzaLNLS](https://discourse.julialang.org/u/VitorSouzaLNLS)\
**Post date:** [February 28, 2024, 10:23pm UTC](https://discourse.julialang.org/t/cuda-jl-for-particle-tracking-simulation/110906/3 "2024-02-28T22:23:31Z")

</div>

Hi! @trahflow ! Thanks for your reply!

After your reply, I ran a test with 256 \* 100 particles (with the GPU function).

Here is the benchmarks GPU x CPU:

`>>> CPU (1 Thread)`

```julia
BenchmarkTools.Trial: 2284 samples with 1 evaluation.
 Range (min … max): 3.222 ms … 30.346 ms ┊ GC (min … max): 0.00% … 82.53%
 Time (median): 4.050 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 4.361 ms ± 2.246 ms ┊ GC (mean ± σ): 6.67% ± 10.65%

  ▇▇▆█▅▄▃▂▁▁ ▁ ▁
  ████████████▆▆▄▁▄▄▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄▁▁▁▁▁▁▁▁▅▁▄▁▁▁▁▁▄█ █
  3.22 ms Histogram: log(frequency) by time 19.3 ms <

 Memory estimate: 1.96 MiB, allocs estimate: 120009.

```

`>>> GPU (256 Threads, 100 blocks)`

```julia
BenchmarkTools.Trial: 12 samples with 1 evaluation.
 Range (min … max): 3.033 s … 4.308 s ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 3.243 s ┊ GC (median): 0.00%
 Time (mean ± σ): 3.290 s ± 331.083 ms ┊ GC (mean ± σ): 0.00% ± 0.00%

           █                                                  
  ▆▁▁▁▄▁▁▁▁█▄▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▄ ▁
  3.03 s Histogram: frequency by time 4.31 s <

 Memory estimate: 4.16 KiB, allocs estimate: 76.

```

`>>> GPU (512 Threads, 50 blocks)`

```julia
BenchmarkTools.Trial: 12 samples with 1 evaluation.
 Range (min … max): 3.250 s … 3.259 s ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 3.256 s ┊ GC (median): 0.00%
 Time (mean ± σ): 3.256 s ± 2.844 ms ┊ GC (mean ± σ): 0.00% ± 0.00%

  █ █ █ █ ██ █ █ █ ███  
  █▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁█▁██▁▁█▁▁▁▁▁▁▁▁█▁▁▁█▁▁███ ▁
  3.25 s Histogram: frequency by time 3.26 s <

 Memory estimate: 4.14 KiB, allocs estimate: 76.

```

`>>> GPU (768 Threads, 33 blocks)`

```julia
BenchmarkTools.Trial: 12 samples with 1 evaluation.
 Range (min … max): 2.475 s … 2.523 s ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 2.493 s ┊ GC (median): 0.00%
 Time (mean ± σ): 2.495 s ± 18.593 ms ┊ GC (mean ± σ): 0.00% ± 0.00%

  ▃ █  
  █▇▇▁▁▁▁▁▁▁▁▁▁▁▁▇▁▁▁▇▁▁▁▇▇▇▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  2.48 s Histogram: frequency by time 2.52 s <

 Memory estimate: 4.11 KiB, allocs estimate: 75.

```

So, in fact, if we take the mean time of the GPU runs and simply divide by the number of particles the time is really speed-up:

_(Per particle)_

- mean CPU = `4.361 ms`
- mean GPU (256 threads) = `0.128 ms`
- mean GPU (512 threads) = `0.127 ms`
- mean GPU (768 threads) = `0.097 ms`

So, indeed its faster!

I still dont have implemented a parallel evaluation of the CPU functions (line\_pass or drift). The benchmarks comparison should be more fair, right? (Instead of comparing the GPU evaluations with 256\*100 particles with CPU eval with 1 particle).

_One thing I still don’t understand is that the benchmark with 512 threads didn’t show any gain compared to the one with 256, but the one with 768 did (as I expected)._

---

<div class="post-metadata">

**Author:** ![trahflow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/trahflow/32/30585_2.png) [@trahflow](https://discourse.julialang.org/u/trahflow)\
**Post date:** [February 29, 2024, 8:42am UTC](https://discourse.julialang.org/t/cuda-jl-for-particle-tracking-simulation/110906/4 "2024-02-29T08:42:02Z")

</div>

> [@VitorSouzaLNLS](#):
>
> One thing I still don’t understand is that the benchmark with 512 threads didn’t show any gain compared to the one with 256, but the one with 768 did (as I expected).

The launch configuration ( 512 Threads, 50 blocks) means that 50 blocks with 512 threads each are being launched.  
This does not necessarily mean however, that only 512 threads run in parallel. Afaik the block scheduler can decide to run multiple blocks at the same time if enough resources for a full block are free. The [CUDA programming guide](https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#multiprocessor-level) might have some more details.

To finetune the launch configuration you would have to do some profiling, as described [here](https://cuda.juliagpu.org/stable/development/profiling/).  
A good starting point however would be to use the occupancy api. The CUDA.jl docs mention this in the introduction at the [bottom of this section](https://cuda.juliagpu.org/stable/tutorials/introduction/#Writing-a-parallel-GPU-kernel)

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [February 29, 2024, 10:40am UTC](https://discourse.julialang.org/t/cuda-jl-for-particle-tracking-simulation/110906/5 "2024-02-29T10:40:35Z")

</div>

You could also try KernelAbstractions.jl which allows you to re-use your code for CPU and GPU.

See [this example](https://juliagpu.github.io/KernelAbstractions.jl/stable/examples/matmul/).