# Correct utilisation of CUDA kernel for simulations

**URL:** <https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882>\
**Category:** GPU\
**Created:** [February 7, 2024, 4:07pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882 "2024-02-07T16:07:34Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [February 7, 2024, 4:07pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/1 "2024-02-07T16:07:34Z")

</div>

Hello,

I am new to julia and GPU computing, before today I didn’t care about optimisation (it was fast enough for me).

I use this kind of CUDA kernel (example for 2D diffusion using finite difference with euler forward):

```julia
using CUDA

function kernel_diff!(ρ_new, ρ, D, Nx, Nz)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    if i <= Nx && j <= Nz
        i_ = mod(i-1,1:Nx); ip = mod(i+1,1:Nx)
        jm = mod(j-1,1:Nz); jp = mod(j+1,1:Nz)
        @inbounds ρ_new[i,j] = ρ[i,j] + D*(ρ[i_,j]+ρ[ip,j]+ρ[i,jm]+ρ[i,jp]-4*ρ[i,j])
    end
    return nothing
end

function diffusion!(ρ, D, Nt)
    Nx, Nz = size(ρ)
    WrapsT = 16
    Bx = ceil(Int, Nx/WrapsT)
    Bz = ceil(Int, Nz/WrapsT)
    block_dim = (WrapsT, WrapsT)
    grid_dim = (Bx, Bz)
    ρ_new = similar(ρ)
    for i=1:Nt
        @cuda threads = block_dim blocks = grid_dim kernel_diff!(ρ_new, ρ, D, Nx, Nz)
        ρ_new = ρ
    end
end

ρ = CUDA.rand(Float64, 1000, 1000)
D = 1e-3
Nt = 2e7

@CUDA.time CUDA.@sync diffusion!(ρ, D, Nt)

```

When I see the number of CPU allocations, i am wondering if i am doing things correctly.

What would I need to change to make it correct ?

Thank you for your help,

Best

PS : I know that there better way to solve diffusion problem, but I use way more complicated kernel to solve more complicated problem with complex boundary conditions. Diffusion is just an example here.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [February 7, 2024, 5:15pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/2 "2024-02-07T17:15:55Z")

</div>

> [@Ludovic\_Dumoulin](#):
>
> `ρ_new = ρ`

Maybe you want `copyto!(ρ, ρ_new)`? And maybe just call `diffusion!(ρ, D, Nt)`

Is this what you want?:

```julia
using CUDA
function kernel_diff!(ρ_new, ρ, D, Nx, Nz)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    if i <= Nx && j <= Nz
        i_ = mod(i-1,1:Nx); ip = mod(i+1,1:Nx)
        jm = mod(j-1,1:Nz); jp = mod(j+1,1:Nz)
        @inbounds ρ_new[i,j] = ρ[i,j] + D*(ρ[i_,j]+ρ[ip,j]+ρ[i,jm]+ρ[i,jp]-4*ρ[i,j])
    end
    return nothing
end
function diffusion!(ρ_new, ρ, D, Nt)
    Nx, Nz = size(ρ)
    gpukernel = @cuda launch=false kernel_diff!(ρ_new, ρ, D, Nx, Nz)
    config = launch_configuration(gpukernel.fun)
    maxThreads = config.threads
    Tx = min(maxThreads, Nx)
    Ty = min(fld(maxThreads, Tx), Nz)
    Bx, By = cld(Nx, Tx), cld(Nz, Ty) # Blocks in grid.
    threads = (Tx, Ty)
    blocks = Bx, By
    for i=1:Nt
        CUDA.@sync gpukernel(ρ_new, ρ, D, Nx, Nz; threads = threads, blocks = blocks)
        copyto!(ρ, ρ_new)
    end
    ρ_new
end

ρ = CUDA.rand(Float64, 1000, 1000)
ρ_new = similar(ρ)
D = 1e-3
Nt = 10000

using BenchmarkTools

@benchmark diffusion!($ρ_new, $ρ, $D, $Nt)

```

```julia
BenchmarkTools.Trial: 2 samples with 1 evaluation.
 Range (min … max): 3.158 s … 3.464 s ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 3.311 s ┊ GC (median): 0.00%
 Time (mean ± σ): 3.311 s ± 216.491 ms ┊ GC (mean ± σ): 0.00% ± 0.00%

  █ █
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  3.16 s Histogram: frequency by time 3.46 s <

 Memory estimate: 44.20 MiB, allocs estimate: 720281.

```

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [February 8, 2024, 9:41am UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/3 "2024-02-08T09:41:23Z")

</div>

Thank you for your reply,

If I use `copyto!()` it takes longer

> 4.956804 seconds (892.41 k CPU allocations: 47.940 MiB, 0.13% gc time) (1 GPU allocation: 7.629 MiB, 0.00% memmgmt time)

(with `ρ_new = ρ`:

> 3.071523 seconds (890.81 k CPU allocations: 47.846 MiB, 0.21% gc time) (1 GPU allocation: 7.629 MiB, 0.00% memmgmt time)

)

In fact I am a bit lost when I see the number of ways to do `ρ_new = ρ` that take very different times…

If I write a kernel to do `ρ_new = ρ`:

```julia
function kernel_update!(ρ, ρ_new, Nx, Nz)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    if i <= Nx && j <= Nz
        @inbounds ρ[i,j] = ρ_new[i,j]
    end
    return nothing
end

```

it takes less time than with `copyto!()` but I have more CPU allocations:

> 4.751879 seconds (1.73 M CPU allocations: 91.440 MiB, 0.36% gc time) (1 GPU allocation: 7.629 MiB, 0.00% memmgmt time)

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [February 8, 2024, 9:57am UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/4 "2024-02-08T09:57:50Z")

</div>

> [@Ludovic\_Dumoulin](#):
>
> > 4.956804 seconds (892.41 k CPU alloca

Hi! Of course, because ‘rho\_new = rho’ really do nothing. If you want to make changes in array - you don’t need rho\_new, if you want to get another array you should use copy or copyto! You can decrease allocations if put sim cycle in kernel.

---

<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 8, 2024, 10:39am UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/5 "2024-02-08T10:39:54Z")

</div>

probably

```julia
ρ, ρ_new = ρ_new, ρ

```

is the most efficient way?

You’re calling your kernel lots of times within a for loop.  
Maybe compare that to another kernel that _contains_ that loop and is called just a _single_ time.  
Alternatively, you can maybe have a look at the [graph execution API](https://cuda.juliagpu.org/stable/lib/driver/#CUDA.@captured) to speed calling the kernel in a loop. I don’t really have experience with this though, so not sure if your case applies.

---

<div class="post-metadata">

**Author:** ![PharmCat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pharmcat/32/6953_2.png) [@PharmCat](https://discourse.julialang.org/u/PharmCat)\
**Post date:** [February 8, 2024, 11:31am UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/6 "2024-02-08T11:31:18Z")

</div>

> [@trahflow](#):
>
> ρ, ρ\_new = ρ\_new, ρ

this is the best solution 🙂  
so, I think that allocation not a problem - GC takes 0.36% of time

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [February 8, 2024, 3:24pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/7 "2024-02-08T15:24:36Z")

</div>

> [@PharmCat](#):
>
> Of course, because ‘rho\_new = rho’ really do nothing.

It’s true ^^, thank you 🙂

> [@trahflow](#):
>
> You’re calling your kernel lots of times within a for loop.  
> Maybe compare that to another kernel that _contains_ that loop and is called just a _single_ time.

The problem is that I need to have the entire array updated before going to the next step.

> [@PharmCat](#):
>
> this is the best solution

Thank you, I’ll try soon

---

<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 8, 2024, 3:49pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/8 "2024-02-08T15:49:23Z")

</div>

> [@Ludovic\_Dumoulin](#):
>
> The problem is that I need to have the entire array updated before going to the next step.

sure, but that doesn’t prevent you from moving the time-step-loop into the kernel. You just need to make sure to place proper synchronization barriers in each iteration.  
That should at least save you the kernel call overhead.  
Not sure whether that will result in any speedup without a proper profiling, but it should be a quick thing to try…

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [February 8, 2024, 4:21pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/9 "2024-02-08T16:21:47Z")

</div>

You could try Stencils.jl for CUDA diffusions… all this is worked out for you and pretty fast.

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [February 8, 2024, 4:30pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/10 "2024-02-08T16:30:53Z")

</div>

Thank you, but in reality I don’t do diffusion.  
I already checked Stencils.jl but it is not suitable for what I am doing

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [February 8, 2024, 4:33pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/11 "2024-02-08T16:33:28Z")

</div>

But… your running some kind of function of a stencil on GPU?

I’m not sure why it wouldn’t help?

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [February 8, 2024, 4:34pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/12 "2024-02-08T16:34:20Z")

</div>

I’ll check, I never tried to synchronize inside a kernel. I usualy call kernel step by step.  
Typically, I perform a lot of operations within a single kernel.

---

<div class="post-metadata">

**Author:** ![Raf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raf/32/3383_2.png) [@Raf](https://discourse.julialang.org/u/Raf)\
**Post date:** [February 8, 2024, 4:37pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/13 "2024-02-08T16:37:57Z")

</div>

You can just use stencils for the window part, you don’t need to use the whole `mapstencil`…

I’m suggesting it because reading a static stencil from an array will be an order of magnitude faster than your loop/ranges with runtime valued size

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [February 9, 2024, 12:47pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/14 "2024-02-09T12:47:50Z")

</div>

Sorry, I made a mistake, I confused it with another package ([ParallelStencil.jl](https://github.com/omlins/ParallelStencil.jl/tree/main?tab=readme-ov-file)).  
There is no real documentation for Stencils.jl, but it seems that it is not suitable when I have multiple field to update in parallel ?  
In this example I have only the density but usually in one kernel I update many fields (and not always the same).

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [February 9, 2024, 2:54pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/15 "2024-02-09T14:54:39Z")

</div>

If I want to synchronize inside a kernel I can synchronize only the threads of one block, then do have I to use only one block ? I was thinking it is impossible for large array

---

<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 9, 2024, 3:19pm UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/16 "2024-02-09T15:19:57Z")

</div>

If the code is not independent across blocks (as is the case for your diffusion example), then this might indeed become difficult.  
In that case I think you should first figure out where the bottleneck actually is and do some profiling.  
There’s a great tutorial on how to do that by @maleadt here: [GitHub - maleadt/cscs2023](https://github.com/maleadt/cscs2023/) (especially 1-3 and 1-4)

---

<div class="post-metadata">

**Author:** ![Ludovic\_Dumoulin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ludovic_dumoulin/32/22414_2.png) [@Ludovic\_Dumoulin](https://discourse.julialang.org/u/Ludovic_Dumoulin)\
**Post date:** [February 13, 2025, 9:53am UTC](https://discourse.julialang.org/t/correct-utilisation-of-cuda-kernel-for-simulations/109882/17 "2025-02-13T09:53:24Z")

</div>

Hello,  
I’m sorry to revive this old conversation. Since my issue is entirely related to this thread and follows its continuity, I assumed this was the right approach.

> [@trahflow](#):
>
> probably
> 
> ```julia
> ρ, ρ_new = ρ_new, ρ
> 
> ```

In `Float64`, `ρ, ρ_new = ρ_new, ρ` works very well, but as soon as I switch to `Float32`, I get this error:

```julia
ERROR: LoadError: MethodError: Cannot `convert` an object of type 
  CuDeviceArray{Float32,2,1} to an object of type 
  CuDeviceArray{Float64,2,1}

```

This is really surprising to me because `copyto!(ρ, ρ_new)` works fine regardless of the float type.

Thank you for your help!  
Best