# Brusselator example from DiffEqGPU won't run or performed badly after simple fix

**URL:** <https://discourse.julialang.org/t/brusselator-example-from-diffeqgpu-wont-run-or-performed-badly-after-simple-fix/119410>\
**Category:** New to Julia\
**Tags:** diffeq, cuda\
**Created:** [September 14, 2024, 8:46pm UTC](https://discourse.julialang.org/t/brusselator-example-from-diffeqgpu-wont-run-or-performed-badly-after-simple-fix/119410 "2024-09-14T20:46:23Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![gluon-free](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gluon-free/32/211942_2.png) [@gluon-free](https://discourse.julialang.org/u/gluon-free)\
**Post date:** [September 14, 2024, 8:46pm UTC](https://discourse.julialang.org/t/brusselator-example-from-diffeqgpu-wont-run-or-performed-badly-after-simple-fix/119410/1 "2024-09-14T20:46:23Z")

</div>

Doing my first steps in Julia, I often find that good examples from “official” sources don’t work because things change so quickly. This [stiff PDE example](https://docs.sciml.ai/DiffEqGPU/stable/examples/bruss/) returning error  
"Scalar indexing is disallowed "

 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/9/f9deb96fef95e28282a5464754d9ef536c095b11.png)  
I googled this type of error and wrapped finite differences with exclusion condition:

`CUDA.allowscalar() do du[II[i, j, 1]] = α * (u[II[im1, j, 1]] + u[II[ip1, j, 1]] + u[II[i, jp1, 1]] + u[II[i, jm1, 1]] - 4u[II[i, j, 1]]) + B + u[II[i, j, 1]]^2 * u[II[i, j, 2]] - (A + 1) * u[II[i, j, 1]] + brusselator_f(x, y, t) end`

Now it works but performance is abysmal. Utilization of my 4090 is ~20-30% and this example is FP32! Computational time ~35 min:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/a/7af9f26a5c4c3b9e64b17e4858e698dc68e30546.png)

What else need to be changed to perform well?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 15, 2024, 7:09am UTC](https://discourse.julialang.org/t/brusselator-example-from-diffeqgpu-wont-run-or-performed-badly-after-simple-fix/119410/2 "2024-09-15T07:09:29Z")

</div>

Yeah if you are removing the allowscalar to get around an error then you are guaranteed bad performance.

@maleadt this seems like a regression in broadcast support in CUDA.jl?

---

<div class="post-metadata">

**Author:** ![maleadt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maleadt/32/10097_2.png) [@maleadt](https://discourse.julialang.org/u/maleadt)\
**Post date:** [September 15, 2024, 11:10am UTC](https://discourse.julialang.org/t/brusselator-example-from-diffeqgpu-wont-run-or-performed-badly-after-simple-fix/119410/3 "2024-09-15T11:10:03Z")

</div>

> [@ChrisRackauckas](#):
>
> @maleadt this seems like a regression in broadcast support in [CUDA.jl](https://juliahub.com/ui/Packages/General/CUDA)?

Hard to tell, e.g., it’s not clear what array types are being used here. If a CPU array is accidentally being used in that expression, it would fall back to a CPU broadcast doing scalar iteration, and that would be expected. If this did work fine on a previous version of CUDA.jl, please file an issue with a minimal example.

---

<div class="post-metadata">

**Author:** ![gluon-free](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gluon-free/32/211942_2.png) [@gluon-free](https://discourse.julialang.org/u/gluon-free)\
**Post date:** [September 15, 2024, 11:11am UTC](https://discourse.julialang.org/t/brusselator-example-from-diffeqgpu-wont-run-or-performed-badly-after-simple-fix/119410/4 "2024-09-15T11:11:40Z")

</div>

> [@maleadt](#):
>
> If this did work fine on a previous version of [CUDA.jl](https://juliahub.com/ui/Packages/General/CUDA), please file an issue with a minimal example.

I dont know, i am newbie and just tried this “seems to be outdated” code.  
I believe that possible fix is turn everything into CUDA arrays. But after that OrdinaryDiffEq.jl and Rosenbrock23() may not work because this is CPU things.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 15, 2024, 12:05pm UTC](https://discourse.julialang.org/t/brusselator-example-from-diffeqgpu-wont-run-or-performed-badly-after-simple-fix/119410/5 "2024-09-15T12:05:24Z")

</div>

It’s just a broadcast over CuArrays

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 15, 2024, 1:08pm UTC](https://discourse.julialang.org/t/brusselator-example-from-diffeqgpu-wont-run-or-performed-badly-after-simple-fix/119410/6 "2024-09-15T13:08:52Z")

</div>

The MWE should be:

```julia
N = 1000
du = cu(rand(Float32,N,N))
u = cu(rand(Float32,N,N))
A = 1f0
B = 1f0
α = 1f0
t = 1f0
II = LinearIndices((N, N, 2))
kernel_u!.(Ref(du), Ref(u), A, B, α, Ref(II), CartesianIndices((N, N)), t)

```

though I don’t have a GPU on me right now to run it, that’s what is scalar indexing.

---

<div class="post-metadata">

**Author:** ![gluon-free](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gluon-free/32/211942_2.png) [@gluon-free](https://discourse.julialang.org/u/gluon-free)\
**Post date:** [September 15, 2024, 5:38pm UTC](https://discourse.julialang.org/t/brusselator-example-from-diffeqgpu-wont-run-or-performed-badly-after-simple-fix/119410/7 "2024-09-15T17:38:25Z")

</div>

Can you please specify how to define that `du` and `u` now became CuArray inside this kernel?

```julia
kernel_u! = let N = N, xyd = xyd_brusselator, dx = step(xyd_brusselator)
    @inline function (du, u, A, B, α, II, I, t)
        i, j = Tuple(I)
        x = xyd[I[1]]
        y = xyd[I[2]]
        ip1 = limit(i + 1, N)
        im1 = limit(i - 1, N)
        jp1 = limit(j + 1, N)
        jm1 = limit(j - 1, N)
        du[II[i, j, 1]] = α * (u[II[im1, j, 1]] + u[II[ip1, j, 1]] + u[II[i, jp1, 1]] +
                           u[II[i, jm1, 1]] - 4u[II[i, j, 1]]) +
                          B + u[II[i, j, 1]]^2 * u[II[i, j, 2]] - (A + 1) * u[II[i, j, 1]] +
                          brusselator_f(x, y, t)
    end
end

```

I did in the head of program:

```julia
u = cu(zeros(Float32,N,N))
du = cu(zeros(Float32,N,N))

```

but this won’t translate inside this kernels.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 16, 2024, 6:03am UTC](https://discourse.julialang.org/t/brusselator-example-from-diffeqgpu-wont-run-or-performed-badly-after-simple-fix/119410/8 "2024-09-16T06:03:37Z")

</div>

What types are you getting in there?