# 2D stationary diffusion problem & (?) MethodOfLines

**URL:** <https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, methodoflines\
**Created:** [September 2, 2024, 9:01pm UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963 "2024-09-02T21:01:05Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [September 2, 2024, 9:01pm UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/1 "2024-09-02T21:01:05Z")

</div>

I’d like to model a mass exchanger which can be represented as diffusion in a fluid flow with given stationary flow field and given coordinate-dependent diffusion coefficient. 2D, rectangular domain, rectangular grid, steady state, all linear. Dirichlet bc on part of borders (inlets), Neumann otherwise.

MethodOfLines.jl apeared to be what I need, however the [first attempt](https://discourse.julialang.org/t/stability-of-steady-state-heat-equation-example-for-methodoflines-jl/98651/8) was unsuccessful, the example from the manual don’t work. I’ve reformulated the same toy problem as `SteadyStateProblem`, it appears to work, and now I have several questions.

> **Code**
>
> ```julia
> using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets, Plots
> 
> @parameters t x y
> @variables u(..)
> Dt = Differential(t)
> Dxx = Differential(x)^2
> Dyy = Differential(y)^2
> 
> eq = Dt(u(t, x, y)) ~ Dxx(u(t, x, y)) + Dyy(u(t, x, y))
> 
> bcs = [u(0, x, y) ~ x * y,
> u(t, 0, y) ~ x * y,
> u(t, 1, y) ~ x * y,
> u(t, x, 0) ~ x * y,
> u(t, x, 1) ~ x * y]
> 
> domains = [t ∈ Interval(0.0, 1.0),
> x ∈ Interval(0.0, 1.0),
> y ∈ Interval(0.0, 1.0)]
> 
> @named pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)])
> 
> g = 20 
> dx = dy = 1/g
> 
> order = 2
> discretization = MOLFiniteDifference([x => dx, y => dy], t)
> 
> t0 = @elapsed prob = discretize(pdesys,discretization)
> 
> steadystateprob = SteadyStateProblem(prob)
> 
> t1 = @elapsed sol = solve(steadystateprob, DynamicSS(Tsit5()))
> t2 = @elapsed sol = solve(steadystateprob, DynamicSS(Tsit5()))
> 
> s = reshape(copy(sol.u), g-1, g-1)
> heatmap(s)
> 
> ```

The solution is returned as 1D array. It is not a problem to reshape it, but why?

The time to compute increases quite steep with grid density. Most of the time is spent in the `discretize` step, and `solve` appears to spend most time by compiling, not computing:  
For the grid 100x100 the timings are

```julia
t0 = 925 # discretize
t1 = 181 # 1st solve
t2 = 0.358 # 2nd solve

```

The total time about 6 min for 75x75 and 18.5 min for 100x100: it looks like the time is proportional to the 4th power of grid side size (i.e. 2nd power of number of nodes).

My system is a bit more complex than the toy example above - should I expect much longer compute time?

I’d like to test the model with different parameter, e.g. different diffusion coefficients. Would it be possible to reuse discretization for some situations? My feeling, the grid like 100x100 or 50x200 should be (just) enough for realistic modelling in my case, but I have little experience, it could do with less.

One alternative would be to write out the finite differences by hand (a bit tedious). Still another possibility is some finite elements package (a lot to learn). Any more options? My purpose is to solve a specific technical problem with spending minimum time.

---

<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 3, 2024, 10:14am UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/2 "2024-09-03T10:14:54Z")

</div>

> [@Eben60](#):
>
> The solution is returned as 1D array. It is not a problem to reshape it, but why?

Just index the solution via the symbolic, i.e. `sol[u]`, to get the shaped version. The reason is because if you have multiple dependent variables, this symbolic form will be able to generate each one and know the interpretation.

> [@Eben60](#):
>
> Most of the time is spent in the `discretize` step, and `solve` appears to spend most time by compiling, not computing:

This is something being worked on by JuliaSimCompiler and symbolic array tooling.

> [@Eben60](#):
>
> I’d like to test the model with different parameter, e.g. different diffusion coefficients. Would it be possible to reuse discretization for some situations? My feeling, the grid like 100x100 or 50x200 should be (just) enough for realistic modelling in my case, but I have little experience, it could do with less.

Just do `remake(prob, p = [myparameter => 1.0))` etc.

> [@Eben60](#):
>
> One alternative would be to write out the finite differences by hand (a bit tedious). Still another possibility is some finite elements package (a lot to learn). Any more options? My purpose is to solve a specific technical problem with spending minimum time.

Right now manual ODE definition code for PDEs would give you faster compile time, though we’re working to remove that gap.

Finite element does not make that much sense if you have a square domains though.

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [September 3, 2024, 9:06pm UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/3 "2024-09-03T21:06:30Z")

</div>

Chris, first thanks a lot for your support!

Now, following issues:

> [@ChrisRackauckas](#):
>
> Just index the solution via the symbolic, i.e. `sol[u]`

```julia
@parameters t x y
@variables u(..)
###
julia> sol[u]
ERROR: ArgumentError: invalid index: u⋆ of type Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}

```

This one is actually just a minor inconvenience in my case.

> [@ChrisRackauckas](#):
>
> Just do `remake(prob, p = [myparameter => 1.0))` etc.

I’ve first tried to run the [Tutorial example](https://docs.sciml.ai/MethodOfLines/stable/tutorials/params/#Adding-parameters), which errored - I’ve opened an [issue](https://github.com/SciML/MethodOfLines.jl/issues/411).

Then tried to apply parameter to my code, with the same error.

> **My code**
>
> using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets, Plots
> 
> @parameters t x y  
> @parameters Px, Py  
> @variables u(…)  
> Dt = Differential(t)  
> Dxx = Differential(x)^2  
> Dyy = Differential(y)^2
> 
> eq = Dt(u(t, x, y)) ~ Dxx(u(t, x, y)) \* (1 - Px + Px \* y^2) + Dyy(u(t, x, y)) \* (1 - Py + Py \* x^2)
> 
> bcs = [u(0, x, y) ~ 1,  
> u(t, 0, y) ~ 0,  
> u(t, 1, y) ~ 0,  
> u(t, x, 0) ~ x, # y = 0  
> u(t, x, 1) ~ 1-x, # y = 1  
> ]
> 
> domains = [t ∈ Interval(0.0, 1.0),  
> x ∈ Interval(0.0, 1.0),  
> y ∈ Interval(0.0, 1.0)]
> 
> @named pdesys = PDESystem(eq, bcs, domains, [t, x, y], [u(t, x, y)], [Px=\>0, Py=\>0])

I hope it was some change in the syntax, not yet reflected in the manual, and not an actual bug.

---

<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 3, 2024, 9:41pm UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/4 "2024-09-03T21:41:27Z")

</div>

> [@Eben60](#):
>
> `sol[u]`

`sol[u(t,x,y)]`

Something happened with parameters in the recent change. I need to handle that.

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [September 4, 2024, 2:11pm UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/5 "2024-09-04T14:11:28Z")

</div>

For parameters, the Tutorial example is working with `MethodOfLines@0.11.0` and erroring with `MethodOfLines@0.11.1`. I’ve commented in the [issue](https://github.com/SciML/MethodOfLines.jl/issues/411).

> [@ChrisRackauckas](#):
>
> Just do `remake(prob, p = [myparameter => 1.0))` etc.

So now I’ve downgraded to 11.0, and it works. I hope now I can start to work on the actual model 🙂

> [@ChrisRackauckas](#):
>
> sol[u(t,x,y)]

```julia
julia> sol[u(t,x,y)]
ERROR: ArgumentError: u(t, x, y) is neither an observed nor an unknown variable.

```

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [September 7, 2024, 3:55pm UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/6 "2024-09-07T15:55:11Z")

</div>

The next problem: [Interfaces](https://docs.sciml.ai/MethodOfLines/dev/boundary_conditions/#Interfaces) (“You may want to connect regions with differing dynamics”)

Tried to run the Tutorial example. As the example contains a problem definition only, added the discretization and solving.

> **code**
>
> ```julia
> using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets
> 
> @parameters t x1 x2
> @variables c1(..)
> @variables c2(..)
> Dt = Differential(t)
> 
> Dx1 = Differential(x1)
> Dxx1 = Dx1^2
> 
> Dx2 = Differential(x2)
> Dxx2 = Dx2^2
> 
> D1(c) = 1 + c / 10
> D2(c) = 1 / 10 + c / 10
> 
> eqs = [Dt(c1(t, x1)) ~ Dx1(D1(c1(t, x1)) * Dx1(c1(t, x1))),
> Dt(c2(t, x2)) ~ Dx2(D2(c2(t, x2)) * Dx2(c2(t, x2)))]
> 
> bcs = [c1(0, x1) ~ 1 + cospi(2 * x1),
> c2(0, x2) ~ 1 + cospi(2 * x2),
> Dx1(c1(t, 0)) ~ 0,
> c1(t, 0.5) ~ c2(t, 0.5), # Relevant interface boundary condition
> -D1(c1(t, 0.5)) * Dx1(c1(t, 0.5)) ~ -D2(c2(t, 0.5)) * Dx2(c2(t, 0.5)), # Higher order interface condition
> Dx2(c2(t, 1)) ~ 0]
> 
> domains = [t ∈ Interval(0.0, 0.15),
> x1 ∈ Interval(0.0, 0.5),
> x2 ∈ Interval(0.5, 1.0)]
> 
> @named pdesys = PDESystem(eqs, bcs, domains,
> [t, x1, x2], [c1(t, x1), c2(t, x2)])
> 
> # from this point my code
> order = 2
> discretization = MOLFiniteDifference([x1 => 0.05, x2 => 0.05], t)
> ;
> prob = MethodOfLines.discretize(pdesys, discretization);
> sol = solve(prob, NewtonRaphson(), saveat = 0.1);
> 
> ```

Getting a warning

```julia
prob = MethodOfLines.discretize(pdesys, discretization);
 Warning: The system contains interface boundaries, which are not compatible with system transformation. The system will not be transformed. Please post an issue if you need this feature.

```

and then an error on calling `solve`.

```julia
ERROR: MethodError: anyeltypedual(::ArrayPartition{Union{}, Tuple{}}) is ambiguous.

```

`MethodOfLines@0.11.0`  
Anything wrong with my code, or with the Tutorial example?

---

<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 8, 2024, 5:07am UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/7 "2024-09-08T05:07:44Z")

</div>

We need to update some tutorials. See the discussion in [MethodOfLines complains about my model... · Issue #249 · SciML/MethodOfLines.jl · GitHub](https://github.com/SciML/MethodOfLines.jl/issues/249)

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [September 9, 2024, 10:30am UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/8 "2024-09-09T10:30:01Z")

</div>

Now updated to MethodOfLines@0.11.3 (i.e. yesterday’s release), adapted my code according to updated manual - it works now ☺

Also got the Tutorial [example with an interface](https://docs.sciml.ai/MethodOfLines/dev/boundary_conditions/#Interfaces) working (just had to change the algorithm in my code to `Rosenbrock23`).

Now trying to adapt that code with interface to my 2D problem. The code below describes a simplified case of a heat/mass exchanger with a counterflow, and semi-permeable wall in the middle.

> **Code counterflow exchanger**
>
> ```julia
> using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets 
> 
> @parameters t x y1 y2
> @variables c1(..)
> @variables c2(..)
> 
> Dt = Differential(t)
> 
> Dx = Differential(x)
> Dxx = Dx^2
> 
> Dy1 = Differential(y1)
> Dyy1 = Dy1^2
> 
> Dy2 = Differential(y2)
> Dyy2 = Dy2^2
> 
> li_vel(y) = 2 * (0.5 - y)
> 
> eq1 = Dt(c1(t, x, y1)) ~ Dxx(c1(t, x, y1)) + Dyy1(c1(t, x, y1)) - Dx(c1(t, x, y1))*li_vel(y1)
> eq2 = Dt(c2(t, x, y2)) ~ Dxx(c2(t, x, y2)) + Dyy2(c2(t, x, y2)) - Dx(c2(t, x, y2))*li_vel(y2)
> 
> eqs = [eq1, eq2]
> 
> domains = [t ∈ Interval(0.0, 0.15),
> x ∈ Interval(0.0, 1.0),
> y1 ∈ Interval(0.0, 0.5),
> y2 ∈ Interval(0.5, 1.0)]
> 
> bcs = [
> c1(0, x, y1) ~ 0.5, # init
> c2(0, x, y2) ~ 0.5,
> c1(t, 0, y1) ~ 0, # flow in
> c2(t, 1, y2) ~ 1,
> Dx(c1(t, 1, y1)) ~ 0, # flow out
> Dx(c2(t, 0, y1)) ~ 0, 
> Dy1(c1(t, x, 0)) ~ 0, # bottom & top walls
> Dy2(c2(t, x, 1)) ~ 0,
> Dy1(c1(t, x, 0.5)) + c1(t, x, 0.5) ~ c2(t, x, 0.5), # membrane
> Dy2(c2(t, x, 0.5)) + c2(t, x, 0.5) ~ c1(t, x, 0.5), 
> ]
> 
> domains = [t ∈ Interval(0.0, 0.15),
> y1 ∈ Interval(0.0, 0.5),
> y2 ∈ Interval(0.5, 1.0)]
> 
> @named pdesys = PDESystem(eqs, bcs, domains,
> [t, x, y1, y2], [c1(t, x, y1), c2(t, x, y2)])
> 
> discretization = MOLFiniteDifference([x => 0.1, y1 => 0.1, y2 => 0.1], t)
> 
> prob = MethodOfLines.discretize(pdesys, discretization) 
> # error on previous line
> alg = Rosenbrock23() 
> sol = solve(prob, alg, saveat = 0.1);
> 
> ```

Unfortunately I get an error on `MethodOfLines.discretize(pdesys, discretization)`:

```julia
ERROR: ArgumentError: invalid index: nothing of type Nothing

```

> **Stack trace**
>
> ```julia
> ERROR: ArgumentError: invalid index: nothing of type Nothing
> Stacktrace:
> [1] to_index(i::Nothing)
> @ Base ./indices.jl:300
> [2] to_index(A::Vector{Symbolics.VarDomainPairing}, i::Nothing)
> @ Base ./indices.jl:277
> [3] _to_indices1(A::Vector{Symbolics.VarDomainPairing}, inds::Tuple{Base.OneTo{Int64}}, I1::Nothing)
> @ Base ./indices.jl:359
> [4] to_indices
> @ ./indices.jl:354 [inlined]
> [5] to_indices
> @ ./indices.jl:345 [inlined]
> [6] getindex
> @ ./abstractarray.jl:1291 [inlined]
> [7] (::PDEBase.var"#54#62"{Vector{Symbolics.VarDomainPairing}})(x::SymbolicUtils.BasicSymbolic{Real})
> @ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/variable_map.jl:34
> [8] iterate
> @ ./generator.jl:47 [inlined]
> [9] collect_to!(dest::Vector{Pair{…}}, itr::Base.Generator{Vector{…}, PDEBase.var"#54#62"{…}}, offs::Int64, st::Int64)
> @ Base ./array.jl:892
> [10] collect_to_with_first!(dest::Vector{…}, v1::Pair{…}, itr::Base.Generator{…}, st::Int64)
> @ Base ./array.jl:870
> [11] _collect(c::Vector{…}, itr::Base.Generator{…}, ::Base.EltypeUnknown, isz::Base.HasShape{…})
> @ Base ./array.jl:864
> [12] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, PDEBase.var"#54#62"{Vector{…}}})
> @ Base ./array.jl:763
> [13] map(f::Function, A::Vector{Any})
> @ Base ./abstractarray.jl:3285
> [14] PDEBase.VariableMap(pdesys::PDESystem, disc::MOLFiniteDifference{…}; replaced_vars::Dict{…})
> @ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/variable_map.jl:33
> [15] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
> @ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/symbolic_discretize.jl:17
> [16] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@Kwargs{})
> @ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/discretization_state.jl:57
> [17] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
> @ PDEBase ~/.julia/packages/PDEBase/o7LGc/src/discretization_state.jl:54
> [18] top-level scope
> @ ~/Julia/MembH2Otransport.jl/src/PDE/m20_2D_middle-interf_4discourse.jl:52
> Some type information was truncated. Use `show(err)` to see complete types.
> 
> ```

Any ideas?  
The manual says

> Note that if you want to use a higher order interface condition, this may not work if you have no simple condition of the form `c1(t, 0.5) ~ c2(t, 0.5)`

I do not have that kind of conditions. However if I add a simple condition,

> **bcs with a simple condition**
>
> ```julia
> bcs = [
> c1(0, x, y1) ~ 0.5, # init
> c2(0, x, y2) ~ 0.5,
> c1(t, 0, y1) ~ 0, # flow in
> c2(t, 1, y2) ~ 1,
> Dx(c1(t, 1, y1)) ~ 0, # flow out
> Dx(c2(t, 0, y1)) ~ 0, 
> Dy1(c1(t, x, 0)) ~ 0, # bottom & top walls
> Dy2(c2(t, x, 1)) ~ 0,
> Dy1(c1(t, x, 0.5)) + c1(t, x, 0.5) ~ c2(t, x, 0.5), # membrane
> # Dy2(c2(t, x, 0.5)) + c2(t, x, 0.5) ~ c1(t, x, 0.5), 
> c2(t, x, 0.5) ~ c1(t, x, 0.5), # simple condition
> ]
> 
> ```

I still get the same error, so it is probably not the error cause, but the next problem.

---

<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 9, 2024, 11:27am UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/9 "2024-09-09T11:27:00Z")

</div>

Okay that’s an odd error 😅. Open an issue on that.

---

<div class="post-metadata">

**Author:** ![Eben60](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eben60/32/13475_2.png) [@Eben60](https://discourse.julialang.org/u/Eben60)\
**Post date:** [September 9, 2024, 1:59pm UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/10 "2024-09-09T13:59:36Z")

</div>

Issue [opened](https://github.com/SciML/MethodOfLines.jl/issues/413)

---

<div class="post-metadata">

**Author:** ![bpeters](https://avatars.discourse-cdn.com/v4/letter/b/ec9cab/32.png) [@bpeters](https://discourse.julialang.org/u/bpeters)\
**Post date:** [December 23, 2025, 10:23am UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/11 "2025-12-23T10:23:20Z")

</div>

@ChrisRackauckas

Dear all,

I am new in the community and use this entry point due to very harsh restrictions as a newcomer (Th at-symbol is replaced by _at_). I see a huge benefit using DiffEqGPU for my applications (luxdem.uni.lu), in particular to use the power of GPUs of different architectures. My application comes from thermodynamics and requires to solve the one-dimensional and transient heat equation (diffusion type differential equation). I have already seen MoL, however, I am not able to use it due to conservation issues and working in different frames of references. However, that is not the issue. Thus, my ode system comes from the discretisation of the one-dimensional heat equation. My application requires to solve app. n=1e6 small and similar ode-systems coming from individual particles. Each small ode system representing a particle consisting of app 10 - 100 unknowns with individual boundary conditions. I am not clear if I should use  
EnsembleGPUKernel for one big ODE system or EnsembleGPUArray to handle each small ODE system. From my understanding, EnsembleGPUKernel should be the method to be applied because

- the function f to evaluate the system i.e. u and would be a loop over all small ode system for which transport properties are dependent on u

- the boundary conditions for each small ode system come from the parameters defined in advance before solving

- is able to account for interactions between the small ode systems

That would then require to run EnsembleGPUKernel for one trajectory. Trying this with the following example taken from the homepage and modified the trajectory number:

```julia-auto
using DiffEqGPU, OrdinaryDiffEq, StaticArrays, CUDA

function lorenz(u, p, t)
σ = p[1]
ρ = p[2]
β = p[3]
du1 = σ * (u[2] - u[1])
du2 = u[1] * (ρ - u[3]) - u[2]
du3 = u[1] * u[2] - β * u[3]
return SVector{3}(du1, du2, du3)
end

u0 = _at_SVector [1.0f0; 0.0f0; 0.0f0]
tspan = (0.0f0, 10.0f0)
p = _at_SVector [10.0f0, 28.0f0, 8 / 3.0f0]
prob = ODEProblem{false}(lorenz, u0, tspan, p)
prob_func = (prob, i, repeat) → remake(prob, p = (_at_SVector rand(Float32, 3)) .* p)
monteprob = EnsembleProblem(prob, prob_func = prob_func, safetycopy = false)
sol = solve(monteprob, GPUTsit5(), EnsembleGPUKernel(CUDA.CUDABackend()), trajectories = 1)

```

gives me (very surprising):

julia\> include(“runODETestsFloat32GPUEnsembleGPUKernel.jl”)  
ERROR: LoadError: MethodError: Cannot `convert` an object of type GPUTsit5 to an object of type Vector{SVector{3, Float32}}  
The function `convert` exists, but no method is defined for this combination of argument types.

Closest candidates are:  
convert(::Type{Array{T, N}}, ::SizedArray{S, T, N, N, Array{T, N}}) where {S, T, N}  
\_at\_ StaticArrays ~/.julia/packages/StaticArrays/DsPgf/src/SizedArray.jl:88  
convert(::Type{Array{T, N}}, ::SizedArray{S, T, N, M, TData} where {M, TData\<:AbstractArray{T, M}}) where {T, S, N}  
\_at\_ StaticArrays ~/.julia/packages/StaticArrays/DsPgf/src/SizedArray.jl:82  
convert(::Type{Array{S, N}}, ::PooledArrays.PooledArray{T, R, N}) where {S, T, R, N}  
\_at\_ PooledArrays ~/.julia/packages/PooledArrays/Vy2X0/src/PooledArrays.jl:499  
…

Stacktrace:  
[1] \_\_init(prob::ODEProblem{…}, alg::CompositeAlgorithm{…}, timeseries\_init::GPUTsit5, ts\_init::Tuple{}, ks\_init::Tuple{}; saveat::Tuple{}, tstops::Tuple{}, d\_discontinuities::Tuple{}, save\_idxs::Nothing, save\_everystep::Bool, save\_on::Bool, save\_discretes::Bool, save\_start::Bool, save\_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float32, dtmin::Float32, dtmax::Float32, force\_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady\_min::Int64, qsteady\_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE\_DEFAULT\_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE\_DEFAULT\_ISOUTOFDOMAIN), unstable\_check::typeof(DiffEqBase.ODE\_DEFAULT\_UNSTABLE\_CHECK), verbose::Bool, timeseries\_errors::Bool, dense\_errors::Bool, advance\_to\_tstop::Bool, stop\_at\_next\_tstop::Bool, initialize\_save::Bool, progress::Bool, progress\_steps::Int64, progress\_name::String, progress\_message::typeof(DiffEqBase.ODE\_DEFAULT\_PROG\_MESSAGE), progress\_id::Symbol, userdata::Nothing, allow\_extrapolation::Bool, initialize\_integrator::Bool, alias::ODEAliasSpecifier, initializealg::DiffEqBase.DefaultInit, kwargs::\_at\_KwargsKwargsKwargsKwargsKwargsKwargsKwargsKwargs{…})  
\_at\_ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/GMkz9/src/solve.jl:334  
[2] \_\_solve(prob::ODEProblem{…}, alg::CompositeAlgorithm{…}, args::GPUTsi\_at\_Kwargs5; kw\_at\_KwargsKwargsrgs::\_at\_Kwargs{…})  
\_at\_ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/GMkz9/src/solve.jl:6  
[3] \_\_solve(prob::ODEProblem{…}, ::Nothing, args::\_at\_KwargsPUTsi\_at\_KwargsKwargs5; kwargs::\_at\_Kwargs{…})  
\_at\_ OrdinaryDiffEqDefault ~/.julia/packages/OrdinaryDiffEqDefault/MdlB6/src/default\_alg.jl:48  
[4] \_\_solve(prob::ODEProblem{…}, args::GPUTsit5; default\_set::Bool, sec\_at\_Kwargsnd\_ti\_at\_KwargsKwargse::Bool, kwargs::\_at\_Kwargs{})  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:758  
[5] \_\_solve(prob::ODEProblem{…}, args::GPUTsit5)  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:749  
[6] solve\_call(\_prob::ODEProblem{…}, args::GPUTsit5; merge\_callbacks::Bool, k\_at\_Kwargsargsh\_at\_KwargsKwargsndle::Nothing, kwargs::\_at\_Kwargs{})  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:142  
[7] solve\_call(\_prob::ODEProblem{…}, args::GPUTsit5)  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:109  
[8] solve\_up(prob::ODEProblem{…}, sensealg::Nothing, u0::SVector{…}, p::SVector{…}, ar\_at\_KwargshainRulesOri\_at\_KwargshainRulesOriginatorinators::G\_at\_KwargsUTsit\_at\_Kwargs; originator::SciMLBase.\_at\_KwargshainRulesOriginator, kwargs::\_at\_Kwargs{})  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:578  
[9] solve(prob::ODEProblem{…}, args::GPUTsit5;\_at\_Kwargssense\_at\_Kwargslg::Nothing, u0::Nothing\_at\_Kwargs p::Nothing, wrap::Val{…}, kwargs::\_at\_Kwargs{})  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:545  
[10] solve(prob::ODEProblem{…}, args::GPUTsit5)  
\_at\_ DiffEqBase ~/.julia/packages/DiffEqBase/6sydR/src/solve.jl:535  
[1\_at\_Kwargs] bat\_at\_Kwargsh\_func(i::Int64, prob::E\_at\_KwargssembleProblem{…}, alg::GPUTsit5; kwargs::\_at\_Kwargs{})  
\_at\_ SciMLBase ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:235  
[12] batch\_func  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:222 [inlined]  
[13] #811  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:294 [inlined]  
[14] responsible\_map  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:287 [inlined]  
[15] solve\_batch(prob::EnsembleProblem{…}, alg::GPUTsit5, ::Ensem\_at\_KwargsKwargsleSerial, II::Unit\_at\_Kwargsange{…}, pmap\_batch\_size::Int64; kwarg\_at\_Kwarg\_at\_Kwarg\_at\_Kwarg\_at\_Kwargs::\_at\_Kwargs{})  
\_at\_ SciMLBase ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:293  
[16] solve\_batch  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_e\_at\_Kwarg\_at\_Kwargssemble\_solve.jl:292 [inlined]  
[17] macro expansion  
\_at\_ ./timi\_at\_Kwarg\_at\_Kwargsg.jl:461 [inlined]  
[1\_at\_Kwarg\_at\_Kwargs] \_at\_Kwargs::SciMLBase.var"#800#801"{Int64, Int64, I\_at\_Kwarg\_at\_Kwarg\_at\_Kwarg\_at\_Kwargst64, \_at\_Kwargs{}, EnsembleProblem{…}, GPUTsit5, EnsembleSerial})()  
\_at\_ SciMLBase ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:188  
[19] with\_logstate(f::SciMLBase.var"#800#801"{…}, logstate::Base.CoreLogging.LogState)  
\_at\_ Base.CoreLogging ./logging/logging.jl:540  
[20] with\_logger  
\_at\_ ./logging/logging.jl:651 [inlined]  
[21] \_\_solve(prob::EnsemblePro\_at\_KwargsInt64lem{…}, alg::GPUTsit5, ens\_at\_KwargsInt64mblealg::EnsembleSerial; traje\_at\_Kwargsnt64t\_at\_Kwargskwargsri\_at\_Kwargss::Int64, batch\_si\_at\_Kw\_at\_Kwargsnt64r\_at\_Kwargskwargs\_at\_K\_at\_KwargsargsInt6\_at\_KwargsInt64kwargse:\_at\_KwargsInt64, progress\_aggregate::Bool, pmap\_batch\_siz\_at\_Kwarg\_at\_Kwargsnt64k\_at\_Kwargskwargsar\_at\_\_at\_Kwargsnt64w\_at\_Kwargskwargsrg\_at\_Kwargss::\_at\_Kwargsnt64,\_at\_Kwargskwargs::\_at\_Kwargs{})  
\_at\_ SciMLBase ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:172  
[22] \_\_solve  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:164 [inlined]  
[23] \_\_solve(ensembleprob::EnsembleProblem{…}, alg::GPUTsit5, ense\_at\_KwargsBoolb\_at\_Kwargskwargsea\_at\_Kwargsg::Ensembl\_at\_KwargsBoolG\_at\_KwargskwargsUK\_at\_Kwargsrnel{…}; tra\_at\_Kwargskwargsec\_at\_Kwargsories::Int64, batch\_size::Int64, unstable\_check::Function, adapti\_at\_Kwarg\_at\_KwargsBoolk\_at\_Kwargskwargsar\_at\_\_at\_KwargsBoolw\_at\_Kwargskwargsrg\_at\_Kwargsse:\_at\_KwargsBool,\_at\_Kwargskwargs::\_at\_Kwargs{})  
\_at\_ DiffEqGPU ~/.julia/packages/DiffEqGPU/JsBHq/src/solve.jl:10  
[24] \_\_solve  
\_at\_ ~/.julia/packages/DiffEqGPU/JsBHq/src/solve.jl:1 [inlined]  
[25] #solve#825  
\_at\_ ~/.julia/packages/SciMLBase/hHrRi/src/ensemble/basic\_ensemble\_solve.jl:359 [inlined]  
[26] top-level scope  
\_at\_ ~/xdemExaScale/xdem/odeSolver/test/runODETestsFloat32GPUEnsembleGPUKernel.jl:19  
[27] include(mapexpr::Function, mod::Module, \_path::String)  
\_at\_ Base ./Base.jl:307  
[28] top-level scope  
\_at\_ REPL[3]:1  
in expression starting at /home/bernhard/xdemExaScale/xdem/odeSolver/test/runODETestsFloat32GPUEnsembleGPUKernel.jl:19  
Some type information was truncated. Use `show(err)` to see complete types.

I hope I made myself clear and appreciate very much your advice.

Best

Bernhard

---

<div class="post-metadata">

**Author:** ![raman\_kumar](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/raman_kumar/32/26782_2.png) [@raman\_kumar](https://discourse.julialang.org/u/raman_kumar)\
**Post date:** [December 23, 2025, 11:02am UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/12 "2025-12-23T11:02:29Z")

</div>

Hi @bpeters , welcome to Julia discourse. I think this guide can help you make better discussion format [Please read: make it easier to help you](https://discourse.julialang.org/t/please-read-make-it-easier-to-help-you/14757).

---

<div class="post-metadata">

**Author:** ![bpeters](https://avatars.discourse-cdn.com/v4/letter/b/ec9cab/32.png) [@bpeters](https://discourse.julialang.org/u/bpeters)\
**Post date:** [December 24, 2025, 6:23am UTC](https://discourse.julialang.org/t/2d-stationary-diffusion-problem-methodoflines/118963/13 "2025-12-24T06:23:02Z")

</div>

Thanks, will do.